Annotation of rpl/lapack/lapack/zgesvj.f, revision 1.9
1.4 bertrand 1: *> \brief <b> ZGESVJ </b>
1.1 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.4 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
8: *> \htmlonly
1.4 bertrand 9: *> Download ZGESVJ + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvj.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvj.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvj.f">
1.1 bertrand 15: *> [TXT]</a>
1.4 bertrand 16: *> \endhtmlonly
1.1 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
22: * LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
1.4 bertrand 23: *
1.1 bertrand 24: * .. Scalar Arguments ..
25: * INTEGER INFO, LDA, LDV, LWORK, LRWORK, M, MV, N
26: * CHARACTER*1 JOBA, JOBU, JOBV
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 A( LDA, * ), V( LDV, * ), CWORK( LWORK )
30: * DOUBLE PRECISION RWORK( LRWORK ), SVA( N )
31: * ..
1.4 bertrand 32: *
1.1 bertrand 33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZGESVJ computes the singular value decomposition (SVD) of a complex
40: *> M-by-N matrix A, where M >= N. The SVD of A is written as
41: *> [++] [xx] [x0] [xx]
42: *> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx]
43: *> [++] [xx]
44: *> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
45: *> matrix, and V is an N-by-N unitary matrix. The diagonal elements
46: *> of SIGMA are the singular values of A. The columns of U and V are the
47: *> left and the right singular vectors of A, respectively.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] JOBA
54: *> \verbatim
1.6 bertrand 55: *> JOBA is CHARACTER*1
1.1 bertrand 56: *> Specifies the structure of A.
57: *> = 'L': The input matrix A is lower triangular;
58: *> = 'U': The input matrix A is upper triangular;
59: *> = 'G': The input matrix A is general M-by-N matrix, M >= N.
60: *> \endverbatim
61: *>
62: *> \param[in] JOBU
63: *> \verbatim
64: *> JOBU is CHARACTER*1
65: *> Specifies whether to compute the left singular vectors
66: *> (columns of U):
1.4 bertrand 67: *> = 'U' or 'F': The left singular vectors corresponding to the nonzero
1.1 bertrand 68: *> singular values are computed and returned in the leading
69: *> columns of A. See more details in the description of A.
70: *> The default numerical orthogonality threshold is set to
1.4 bertrand 71: *> approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=DLAMCH('E').
1.1 bertrand 72: *> = 'C': Analogous to JOBU='U', except that user can control the
73: *> level of numerical orthogonality of the computed left
74: *> singular vectors. TOL can be set to TOL = CTOL*EPS, where
75: *> CTOL is given on input in the array WORK.
76: *> No CTOL smaller than ONE is allowed. CTOL greater
77: *> than 1 / EPS is meaningless. The option 'C'
78: *> can be used if M*EPS is satisfactory orthogonality
79: *> of the computed left singular vectors, so CTOL=M could
80: *> save few sweeps of Jacobi rotations.
81: *> See the descriptions of A and WORK(1).
82: *> = 'N': The matrix U is not computed. However, see the
83: *> description of A.
84: *> \endverbatim
85: *>
86: *> \param[in] JOBV
87: *> \verbatim
88: *> JOBV is CHARACTER*1
89: *> Specifies whether to compute the right singular vectors, that
90: *> is, the matrix V:
1.4 bertrand 91: *> = 'V' or 'J': the matrix V is computed and returned in the array V
1.8 bertrand 92: *> = 'A': the Jacobi rotations are applied to the MV-by-N
1.1 bertrand 93: *> array V. In other words, the right singular vector
1.4 bertrand 94: *> matrix V is not computed explicitly; instead it is
1.1 bertrand 95: *> applied to an MV-by-N matrix initially stored in the
96: *> first MV rows of V.
1.8 bertrand 97: *> = 'N': the matrix V is not computed and the array V is not
1.1 bertrand 98: *> referenced
99: *> \endverbatim
100: *>
101: *> \param[in] M
102: *> \verbatim
103: *> M is INTEGER
1.4 bertrand 104: *> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.
1.1 bertrand 105: *> \endverbatim
106: *>
107: *> \param[in] N
108: *> \verbatim
109: *> N is INTEGER
110: *> The number of columns of the input matrix A.
111: *> M >= N >= 0.
112: *> \endverbatim
113: *>
114: *> \param[in,out] A
115: *> \verbatim
116: *> A is COMPLEX*16 array, dimension (LDA,N)
117: *> On entry, the M-by-N matrix A.
118: *> On exit,
1.8 bertrand 119: *> If JOBU = 'U' .OR. JOBU = 'C':
120: *> If INFO = 0 :
1.1 bertrand 121: *> RANKA orthonormal columns of U are returned in the
122: *> leading RANKA columns of the array A. Here RANKA <= N
123: *> is the number of computed singular values of A that are
124: *> above the underflow threshold DLAMCH('S'). The singular
125: *> vectors corresponding to underflowed or zero singular
126: *> values are not computed. The value of RANKA is returned
127: *> in the array RWORK as RANKA=NINT(RWORK(2)). Also see the
128: *> descriptions of SVA and RWORK. The computed columns of U
129: *> are mutually numerically orthogonal up to approximately
1.8 bertrand 130: *> TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'),
1.1 bertrand 131: *> see the description of JOBU.
1.8 bertrand 132: *> If INFO > 0,
1.1 bertrand 133: *> the procedure ZGESVJ did not converge in the given number
134: *> of iterations (sweeps). In that case, the computed
135: *> columns of U may not be orthogonal up to TOL. The output
136: *> U (stored in A), SIGMA (given by the computed singular
137: *> values in SVA(1:N)) and V is still a decomposition of the
138: *> input matrix A in the sense that the residual
139: *> || A - SCALE * U * SIGMA * V^* ||_2 / ||A||_2 is small.
1.8 bertrand 140: *> If JOBU = 'N':
141: *> If INFO = 0 :
1.1 bertrand 142: *> Note that the left singular vectors are 'for free' in the
143: *> one-sided Jacobi SVD algorithm. However, if only the
144: *> singular values are needed, the level of numerical
145: *> orthogonality of U is not an issue and iterations are
146: *> stopped when the columns of the iterated matrix are
147: *> numerically orthogonal up to approximately M*EPS. Thus,
148: *> on exit, A contains the columns of U scaled with the
149: *> corresponding singular values.
1.8 bertrand 150: *> If INFO > 0:
1.1 bertrand 151: *> the procedure ZGESVJ did not converge in the given number
152: *> of iterations (sweeps).
153: *> \endverbatim
154: *>
155: *> \param[in] LDA
156: *> \verbatim
157: *> LDA is INTEGER
158: *> The leading dimension of the array A. LDA >= max(1,M).
159: *> \endverbatim
160: *>
161: *> \param[out] SVA
162: *> \verbatim
163: *> SVA is DOUBLE PRECISION array, dimension (N)
164: *> On exit,
1.8 bertrand 165: *> If INFO = 0 :
1.1 bertrand 166: *> depending on the value SCALE = RWORK(1), we have:
1.8 bertrand 167: *> If SCALE = ONE:
1.1 bertrand 168: *> SVA(1:N) contains the computed singular values of A.
169: *> During the computation SVA contains the Euclidean column
170: *> norms of the iterated matrices in the array A.
171: *> If SCALE .NE. ONE:
172: *> The singular values of A are SCALE*SVA(1:N), and this
173: *> factored representation is due to the fact that some of the
174: *> singular values of A might underflow or overflow.
175: *>
1.8 bertrand 176: *> If INFO > 0:
1.1 bertrand 177: *> the procedure ZGESVJ did not converge in the given number of
178: *> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
179: *> \endverbatim
180: *>
181: *> \param[in] MV
182: *> \verbatim
183: *> MV is INTEGER
1.8 bertrand 184: *> If JOBV = 'A', then the product of Jacobi rotations in ZGESVJ
1.1 bertrand 185: *> is applied to the first MV rows of V. See the description of JOBV.
186: *> \endverbatim
187: *>
188: *> \param[in,out] V
189: *> \verbatim
190: *> V is COMPLEX*16 array, dimension (LDV,N)
191: *> If JOBV = 'V', then V contains on exit the N-by-N matrix of
192: *> the right singular vectors;
193: *> If JOBV = 'A', then V contains the product of the computed right
194: *> singular vector matrix and the initial matrix in
195: *> the array V.
196: *> If JOBV = 'N', then V is not referenced.
197: *> \endverbatim
198: *>
199: *> \param[in] LDV
200: *> \verbatim
201: *> LDV is INTEGER
1.8 bertrand 202: *> The leading dimension of the array V, LDV >= 1.
203: *> If JOBV = 'V', then LDV >= max(1,N).
204: *> If JOBV = 'A', then LDV >= max(1,MV) .
1.1 bertrand 205: *> \endverbatim
206: *>
207: *> \param[in,out] CWORK
208: *> \verbatim
1.6 bertrand 209: *> CWORK is COMPLEX*16 array, dimension (max(1,LWORK))
1.4 bertrand 210: *> Used as workspace.
1.8 bertrand 211: *> If on entry LWORK = -1, then a workspace query is assumed and
1.4 bertrand 212: *> no computation is done; CWORK(1) is set to the minial (and optimal)
213: *> length of CWORK.
1.1 bertrand 214: *> \endverbatim
215: *>
216: *> \param[in] LWORK
217: *> \verbatim
218: *> LWORK is INTEGER.
219: *> Length of CWORK, LWORK >= M+N.
220: *> \endverbatim
221: *>
222: *> \param[in,out] RWORK
223: *> \verbatim
1.6 bertrand 224: *> RWORK is DOUBLE PRECISION array, dimension (max(6,LRWORK))
1.1 bertrand 225: *> On entry,
1.8 bertrand 226: *> If JOBU = 'C' :
1.1 bertrand 227: *> RWORK(1) = CTOL, where CTOL defines the threshold for convergence.
228: *> The process stops if all columns of A are mutually
229: *> orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
230: *> It is required that CTOL >= ONE, i.e. it is not
231: *> allowed to force the routine to obtain orthogonality
232: *> below EPSILON.
233: *> On exit,
234: *> RWORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
235: *> are the computed singular values of A.
236: *> (See description of SVA().)
237: *> RWORK(2) = NINT(RWORK(2)) is the number of the computed nonzero
238: *> singular values.
239: *> RWORK(3) = NINT(RWORK(3)) is the number of the computed singular
240: *> values that are larger than the underflow threshold.
241: *> RWORK(4) = NINT(RWORK(4)) is the number of sweeps of Jacobi
242: *> rotations needed for numerical convergence.
243: *> RWORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
244: *> This is useful information in cases when ZGESVJ did
245: *> not converge, as it can be used to estimate whether
1.8 bertrand 246: *> the output is still useful and for post festum analysis.
1.1 bertrand 247: *> RWORK(6) = the largest absolute value over all sines of the
248: *> Jacobi rotation angles in the last sweep. It can be
249: *> useful for a post festum analysis.
1.8 bertrand 250: *> If on entry LRWORK = -1, then a workspace query is assumed and
1.4 bertrand 251: *> no computation is done; RWORK(1) is set to the minial (and optimal)
252: *> length of RWORK.
1.1 bertrand 253: *> \endverbatim
254: *>
255: *> \param[in] LRWORK
256: *> \verbatim
1.4 bertrand 257: *> LRWORK is INTEGER
1.1 bertrand 258: *> Length of RWORK, LRWORK >= MAX(6,N).
259: *> \endverbatim
260: *>
261: *> \param[out] INFO
262: *> \verbatim
263: *> INFO is INTEGER
1.8 bertrand 264: *> = 0: successful exit.
265: *> < 0: if INFO = -i, then the i-th argument had an illegal value
266: *> > 0: ZGESVJ did not converge in the maximal allowed number
1.4 bertrand 267: *> (NSWEEP=30) of sweeps. The output may still be useful.
1.1 bertrand 268: *> See the description of RWORK.
269: *> \endverbatim
270: *>
271: * Authors:
272: * ========
273: *
1.4 bertrand 274: *> \author Univ. of Tennessee
275: *> \author Univ. of California Berkeley
276: *> \author Univ. of Colorado Denver
277: *> \author NAG Ltd.
1.1 bertrand 278: *
1.4 bertrand 279: *> \ingroup complex16GEcomputational
1.1 bertrand 280: *
281: *> \par Further Details:
282: * =====================
283: *>
284: *> \verbatim
285: *>
286: *> The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
287: *> rotations. In the case of underflow of the tangent of the Jacobi angle, a
288: *> modified Jacobi transformation of Drmac [3] is used. Pivot strategy uses
289: *> column interchanges of de Rijk [1]. The relative accuracy of the computed
290: *> singular values and the accuracy of the computed singular vectors (in
291: *> angle metric) is as guaranteed by the theory of Demmel and Veselic [2].
292: *> The condition number that determines the accuracy in the full rank case
293: *> is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
294: *> spectral condition number. The best performance of this Jacobi SVD
295: *> procedure is achieved if used in an accelerated version of Drmac and
296: *> Veselic [4,5], and it is the kernel routine in the SIGMA library [6].
1.9 ! bertrand 297: *> Some tuning parameters (marked with [TP]) are available for the
1.4 bertrand 298: *> implementer.
1.1 bertrand 299: *> The computational range for the nonzero singular values is the machine
300: *> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
301: *> denormalized singular values can be computed with the corresponding
302: *> gradual loss of accurate digits.
303: *> \endverbatim
304: *
1.4 bertrand 305: *> \par Contributor:
1.1 bertrand 306: * ==================
307: *>
308: *> \verbatim
309: *>
310: *> ============
311: *>
1.4 bertrand 312: *> Zlatko Drmac (Zagreb, Croatia)
313: *>
1.1 bertrand 314: *> \endverbatim
315: *
316: *> \par References:
317: * ================
318: *>
1.6 bertrand 319: *> \verbatim
320: *>
1.1 bertrand 321: *> [1] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
1.4 bertrand 322: *> singular value decomposition on a vector computer.
323: *> SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
1.1 bertrand 324: *> [2] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
325: *> [3] Z. Drmac: Implementation of Jacobi rotations for accurate singular
326: *> value computation in floating point arithmetic.
327: *> SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
328: *> [4] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
329: *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
330: *> LAPACK Working note 169.
331: *> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
332: *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
333: *> LAPACK Working note 170.
334: *> [6] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
335: *> QSVD, (H,K)-SVD computations.
336: *> Department of Mathematics, University of Zagreb, 2008, 2015.
337: *> \endverbatim
338: *
1.4 bertrand 339: *> \par Bugs, examples and comments:
340: * =================================
1.1 bertrand 341: *>
342: *> \verbatim
343: *> ===========================
344: *> Please report all bugs and send interesting test examples and comments to
345: *> drmac@math.hr. Thank you.
346: *> \endverbatim
347: *>
348: * =====================================================================
1.4 bertrand 349: SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
1.1 bertrand 350: $ LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
351: *
1.9 ! bertrand 352: * -- LAPACK computational routine --
1.1 bertrand 353: * -- LAPACK is a software package provided by Univ. of Tennessee, --
354: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
355: *
1.4 bertrand 356: IMPLICIT NONE
1.1 bertrand 357: * .. Scalar Arguments ..
358: INTEGER INFO, LDA, LDV, LWORK, LRWORK, M, MV, N
359: CHARACTER*1 JOBA, JOBU, JOBV
360: * ..
361: * .. Array Arguments ..
362: COMPLEX*16 A( LDA, * ), V( LDV, * ), CWORK( LWORK )
363: DOUBLE PRECISION RWORK( LRWORK ), SVA( N )
364: * ..
365: *
366: * =====================================================================
367: *
368: * .. Local Parameters ..
369: DOUBLE PRECISION ZERO, HALF, ONE
370: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0)
371: COMPLEX*16 CZERO, CONE
372: PARAMETER ( CZERO = (0.0D0, 0.0D0), CONE = (1.0D0, 0.0D0) )
373: INTEGER NSWEEP
374: PARAMETER ( NSWEEP = 30 )
375: * ..
376: * .. Local Scalars ..
377: COMPLEX*16 AAPQ, OMPQ
1.4 bertrand 378: DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
379: $ BIGTHETA, CS, CTOL, EPSLN, MXAAPQ,
380: $ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
381: $ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, THSIGN, TOL
1.1 bertrand 382: INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
1.4 bertrand 383: $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
1.1 bertrand 384: $ N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
1.4 bertrand 385: LOGICAL APPLV, GOSCALE, LOWER, LQUERY, LSVEC, NOSCALE, ROTOK,
1.1 bertrand 386: $ RSVEC, UCTOL, UPPER
387: * ..
388: * ..
389: * .. Intrinsic Functions ..
1.4 bertrand 390: INTRINSIC ABS, MAX, MIN, CONJG, DBLE, SIGN, SQRT
1.1 bertrand 391: * ..
392: * .. External Functions ..
393: * ..
394: * from BLAS
395: DOUBLE PRECISION DZNRM2
396: COMPLEX*16 ZDOTC
397: EXTERNAL ZDOTC, DZNRM2
398: INTEGER IDAMAX
399: EXTERNAL IDAMAX
400: * from LAPACK
401: DOUBLE PRECISION DLAMCH
402: EXTERNAL DLAMCH
403: LOGICAL LSAME
404: EXTERNAL LSAME
405: * ..
406: * .. External Subroutines ..
407: * ..
408: * from BLAS
1.6 bertrand 409: EXTERNAL ZCOPY, ZROT, ZDSCAL, ZSWAP, ZAXPY
1.1 bertrand 410: * from LAPACK
1.2 bertrand 411: EXTERNAL DLASCL, ZLASCL, ZLASET, ZLASSQ, XERBLA
1.1 bertrand 412: EXTERNAL ZGSVJ0, ZGSVJ1
413: * ..
414: * .. Executable Statements ..
415: *
416: * Test the input arguments
417: *
1.4 bertrand 418: LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
1.1 bertrand 419: UCTOL = LSAME( JOBU, 'C' )
1.4 bertrand 420: RSVEC = LSAME( JOBV, 'V' ) .OR. LSAME( JOBV, 'J' )
1.1 bertrand 421: APPLV = LSAME( JOBV, 'A' )
422: UPPER = LSAME( JOBA, 'U' )
423: LOWER = LSAME( JOBA, 'L' )
424: *
1.4 bertrand 425: LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
1.1 bertrand 426: IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
427: INFO = -1
428: ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
429: INFO = -2
430: ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
431: INFO = -3
432: ELSE IF( M.LT.0 ) THEN
433: INFO = -4
434: ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
435: INFO = -5
436: ELSE IF( LDA.LT.M ) THEN
437: INFO = -7
438: ELSE IF( MV.LT.0 ) THEN
439: INFO = -9
440: ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
441: $ ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
442: INFO = -11
443: ELSE IF( UCTOL .AND. ( RWORK( 1 ).LE.ONE ) ) THEN
444: INFO = -12
1.4 bertrand 445: ELSE IF( ( LWORK.LT.( M+N ) ) .AND. ( .NOT.LQUERY ) ) THEN
1.1 bertrand 446: INFO = -13
1.4 bertrand 447: ELSE IF( ( LRWORK.LT.MAX( N, 6 ) ) .AND. ( .NOT.LQUERY ) ) THEN
448: INFO = -15
1.1 bertrand 449: ELSE
450: INFO = 0
451: END IF
452: *
453: * #:(
454: IF( INFO.NE.0 ) THEN
455: CALL XERBLA( 'ZGESVJ', -INFO )
456: RETURN
1.4 bertrand 457: ELSE IF ( LQUERY ) THEN
458: CWORK(1) = M + N
459: RWORK(1) = MAX( N, 6 )
460: RETURN
1.1 bertrand 461: END IF
462: *
463: * #:) Quick return for void matrix
464: *
465: IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
466: *
467: * Set numerical parameters
468: * The stopping criterion for Jacobi rotations is
469: *
470: * max_{i<>j}|A(:,i)^* * A(:,j)| / (||A(:,i)||*||A(:,j)||) < CTOL*EPS
471: *
472: * where EPS is the round-off and CTOL is defined as follows:
473: *
474: IF( UCTOL ) THEN
475: * ... user controlled
476: CTOL = RWORK( 1 )
477: ELSE
478: * ... default
479: IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
1.4 bertrand 480: CTOL = SQRT( DBLE( M ) )
1.1 bertrand 481: ELSE
1.2 bertrand 482: CTOL = DBLE( M )
1.1 bertrand 483: END IF
484: END IF
485: * ... and the machine dependent parameters are
1.4 bertrand 486: *[!] (Make sure that SLAMCH() works properly on the target machine.)
1.1 bertrand 487: *
488: EPSLN = DLAMCH( 'Epsilon' )
1.4 bertrand 489: ROOTEPS = SQRT( EPSLN )
1.1 bertrand 490: SFMIN = DLAMCH( 'SafeMinimum' )
1.4 bertrand 491: ROOTSFMIN = SQRT( SFMIN )
1.1 bertrand 492: SMALL = SFMIN / EPSLN
493: BIG = DLAMCH( 'Overflow' )
494: * BIG = ONE / SFMIN
495: ROOTBIG = ONE / ROOTSFMIN
1.4 bertrand 496: * LARGE = BIG / SQRT( DBLE( M*N ) )
1.1 bertrand 497: BIGTHETA = ONE / ROOTEPS
498: *
499: TOL = CTOL*EPSLN
1.4 bertrand 500: ROOTTOL = SQRT( TOL )
1.1 bertrand 501: *
1.2 bertrand 502: IF( DBLE( M )*EPSLN.GE.ONE ) THEN
1.1 bertrand 503: INFO = -4
504: CALL XERBLA( 'ZGESVJ', -INFO )
505: RETURN
506: END IF
507: *
508: * Initialize the right singular vector matrix.
509: *
510: IF( RSVEC ) THEN
511: MVL = N
512: CALL ZLASET( 'A', MVL, N, CZERO, CONE, V, LDV )
513: ELSE IF( APPLV ) THEN
514: MVL = MV
515: END IF
516: RSVEC = RSVEC .OR. APPLV
517: *
518: * Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
519: *(!) If necessary, scale A to protect the largest singular value
520: * from overflow. It is possible that saving the largest singular
521: * value destroys the information about the small ones.
522: * This initial scaling is almost minimal in the sense that the
523: * goal is to make sure that no column norm overflows, and that
524: * SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
525: * in A are detected, the procedure returns with INFO=-6.
526: *
1.4 bertrand 527: SKL = ONE / SQRT( DBLE( M )*DBLE( N ) )
1.1 bertrand 528: NOSCALE = .TRUE.
529: GOSCALE = .TRUE.
530: *
531: IF( LOWER ) THEN
532: * the input matrix is M-by-N lower triangular (trapezoidal)
533: DO 1874 p = 1, N
534: AAPP = ZERO
535: AAQQ = ONE
536: CALL ZLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
537: IF( AAPP.GT.BIG ) THEN
538: INFO = -6
539: CALL XERBLA( 'ZGESVJ', -INFO )
540: RETURN
541: END IF
1.4 bertrand 542: AAQQ = SQRT( AAQQ )
1.1 bertrand 543: IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
544: SVA( p ) = AAPP*AAQQ
545: ELSE
546: NOSCALE = .FALSE.
547: SVA( p ) = AAPP*( AAQQ*SKL )
548: IF( GOSCALE ) THEN
549: GOSCALE = .FALSE.
550: DO 1873 q = 1, p - 1
551: SVA( q ) = SVA( q )*SKL
552: 1873 CONTINUE
553: END IF
554: END IF
555: 1874 CONTINUE
556: ELSE IF( UPPER ) THEN
557: * the input matrix is M-by-N upper triangular (trapezoidal)
558: DO 2874 p = 1, N
559: AAPP = ZERO
560: AAQQ = ONE
561: CALL ZLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
562: IF( AAPP.GT.BIG ) THEN
563: INFO = -6
564: CALL XERBLA( 'ZGESVJ', -INFO )
565: RETURN
566: END IF
1.4 bertrand 567: AAQQ = SQRT( AAQQ )
1.1 bertrand 568: IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
569: SVA( p ) = AAPP*AAQQ
570: ELSE
571: NOSCALE = .FALSE.
572: SVA( p ) = AAPP*( AAQQ*SKL )
573: IF( GOSCALE ) THEN
574: GOSCALE = .FALSE.
575: DO 2873 q = 1, p - 1
576: SVA( q ) = SVA( q )*SKL
577: 2873 CONTINUE
578: END IF
579: END IF
580: 2874 CONTINUE
581: ELSE
582: * the input matrix is M-by-N general dense
583: DO 3874 p = 1, N
584: AAPP = ZERO
585: AAQQ = ONE
586: CALL ZLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
587: IF( AAPP.GT.BIG ) THEN
588: INFO = -6
589: CALL XERBLA( 'ZGESVJ', -INFO )
590: RETURN
591: END IF
1.4 bertrand 592: AAQQ = SQRT( AAQQ )
1.1 bertrand 593: IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
594: SVA( p ) = AAPP*AAQQ
595: ELSE
596: NOSCALE = .FALSE.
597: SVA( p ) = AAPP*( AAQQ*SKL )
598: IF( GOSCALE ) THEN
599: GOSCALE = .FALSE.
600: DO 3873 q = 1, p - 1
601: SVA( q ) = SVA( q )*SKL
602: 3873 CONTINUE
603: END IF
604: END IF
605: 3874 CONTINUE
606: END IF
607: *
608: IF( NOSCALE )SKL = ONE
609: *
610: * Move the smaller part of the spectrum from the underflow threshold
611: *(!) Start by determining the position of the nonzero entries of the
612: * array SVA() relative to ( SFMIN, BIG ).
613: *
614: AAPP = ZERO
615: AAQQ = BIG
616: DO 4781 p = 1, N
1.4 bertrand 617: IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) )
618: AAPP = MAX( AAPP, SVA( p ) )
1.1 bertrand 619: 4781 CONTINUE
620: *
621: * #:) Quick return for zero matrix
622: *
623: IF( AAPP.EQ.ZERO ) THEN
624: IF( LSVEC )CALL ZLASET( 'G', M, N, CZERO, CONE, A, LDA )
625: RWORK( 1 ) = ONE
626: RWORK( 2 ) = ZERO
627: RWORK( 3 ) = ZERO
628: RWORK( 4 ) = ZERO
629: RWORK( 5 ) = ZERO
630: RWORK( 6 ) = ZERO
631: RETURN
632: END IF
633: *
634: * #:) Quick return for one-column matrix
635: *
636: IF( N.EQ.1 ) THEN
637: IF( LSVEC )CALL ZLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
638: $ A( 1, 1 ), LDA, IERR )
639: RWORK( 1 ) = ONE / SKL
640: IF( SVA( 1 ).GE.SFMIN ) THEN
641: RWORK( 2 ) = ONE
642: ELSE
643: RWORK( 2 ) = ZERO
644: END IF
645: RWORK( 3 ) = ZERO
646: RWORK( 4 ) = ZERO
647: RWORK( 5 ) = ZERO
648: RWORK( 6 ) = ZERO
649: RETURN
650: END IF
651: *
652: * Protect small singular values from underflow, and try to
653: * avoid underflows/overflows in computing Jacobi rotations.
654: *
1.4 bertrand 655: SN = SQRT( SFMIN / EPSLN )
656: TEMP1 = SQRT( BIG / DBLE( N ) )
657: IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
1.1 bertrand 658: $ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
1.4 bertrand 659: TEMP1 = MIN( BIG, TEMP1 / AAPP )
1.1 bertrand 660: * AAQQ = AAQQ*TEMP1
661: * AAPP = AAPP*TEMP1
662: ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
1.4 bertrand 663: TEMP1 = MIN( SN / AAQQ, BIG / (AAPP*SQRT( DBLE(N)) ) )
1.1 bertrand 664: * AAQQ = AAQQ*TEMP1
665: * AAPP = AAPP*TEMP1
666: ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
1.4 bertrand 667: TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP )
1.1 bertrand 668: * AAQQ = AAQQ*TEMP1
669: * AAPP = AAPP*TEMP1
670: ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
1.4 bertrand 671: TEMP1 = MIN( SN / AAQQ, BIG / ( SQRT( DBLE( N ) )*AAPP ) )
1.1 bertrand 672: * AAQQ = AAQQ*TEMP1
673: * AAPP = AAPP*TEMP1
674: ELSE
675: TEMP1 = ONE
676: END IF
677: *
678: * Scale, if necessary
679: *
680: IF( TEMP1.NE.ONE ) THEN
1.2 bertrand 681: CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
1.1 bertrand 682: END IF
683: SKL = TEMP1*SKL
684: IF( SKL.NE.ONE ) THEN
685: CALL ZLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
686: SKL = ONE / SKL
687: END IF
688: *
689: * Row-cyclic Jacobi SVD algorithm with column pivoting
690: *
691: EMPTSW = ( N*( N-1 ) ) / 2
692: NOTROT = 0
1.4 bertrand 693:
1.1 bertrand 694: DO 1868 q = 1, N
695: CWORK( q ) = CONE
1.4 bertrand 696: 1868 CONTINUE
1.1 bertrand 697: *
698: *
699: *
700: SWBAND = 3
701: *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
702: * if ZGESVJ is used as a computational routine in the preconditioned
703: * Jacobi SVD algorithm ZGEJSV. For sweeps i=1:SWBAND the procedure
704: * works on pivots inside a band-like region around the diagonal.
705: * The boundaries are determined dynamically, based on the number of
706: * pivots above a threshold.
707: *
1.4 bertrand 708: KBL = MIN( 8, N )
1.1 bertrand 709: *[TP] KBL is a tuning parameter that defines the tile size in the
710: * tiling of the p-q loops of pivot pairs. In general, an optimal
711: * value of KBL depends on the matrix dimensions and on the
712: * parameters of the computer's memory.
713: *
714: NBL = N / KBL
715: IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
716: *
717: BLSKIP = KBL**2
718: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
719: *
1.4 bertrand 720: ROWSKIP = MIN( 5, KBL )
1.1 bertrand 721: *[TP] ROWSKIP is a tuning parameter.
722: *
723: LKAHEAD = 1
724: *[TP] LKAHEAD is a tuning parameter.
725: *
726: * Quasi block transformations, using the lower (upper) triangular
727: * structure of the input matrix. The quasi-block-cycling usually
728: * invokes cubic convergence. Big part of this cycle is done inside
729: * canonical subspaces of dimensions less than M.
730: *
1.4 bertrand 731: IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN
1.1 bertrand 732: *[TP] The number of partition levels and the actual partition are
733: * tuning parameters.
734: N4 = N / 4
735: N2 = N / 2
736: N34 = 3*N4
737: IF( APPLV ) THEN
738: q = 0
739: ELSE
740: q = 1
741: END IF
742: *
743: IF( LOWER ) THEN
744: *
745: * This works very well on lower triangular matrices, in particular
746: * in the framework of the preconditioned Jacobi SVD (xGEJSV).
747: * The idea is simple:
748: * [+ 0 0 0] Note that Jacobi transformations of [0 0]
749: * [+ + 0 0] [0 0]
750: * [+ + x 0] actually work on [x 0] [x 0]
751: * [+ + x x] [x x]. [x x]
752: *
753: CALL ZGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
754: $ CWORK( N34+1 ), SVA( N34+1 ), MVL,
755: $ V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
756: $ 2, CWORK( N+1 ), LWORK-N, IERR )
757:
758: CALL ZGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
759: $ CWORK( N2+1 ), SVA( N2+1 ), MVL,
760: $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
761: $ CWORK( N+1 ), LWORK-N, IERR )
762:
763: CALL ZGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
764: $ CWORK( N2+1 ), SVA( N2+1 ), MVL,
765: $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
766: $ CWORK( N+1 ), LWORK-N, IERR )
767:
768: CALL ZGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
769: $ CWORK( N4+1 ), SVA( N4+1 ), MVL,
770: $ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
771: $ CWORK( N+1 ), LWORK-N, IERR )
772: *
773: CALL ZGSVJ0( JOBV, M, N4, A, LDA, CWORK, SVA, MVL, V, LDV,
774: $ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N,
775: $ IERR )
776: *
777: CALL ZGSVJ1( JOBV, M, N2, N4, A, LDA, CWORK, SVA, MVL, V,
778: $ LDV, EPSLN, SFMIN, TOL, 1, CWORK( N+1 ),
779: $ LWORK-N, IERR )
780: *
781: *
782: ELSE IF( UPPER ) THEN
783: *
784: *
785: CALL ZGSVJ0( JOBV, N4, N4, A, LDA, CWORK, SVA, MVL, V, LDV,
786: $ EPSLN, SFMIN, TOL, 2, CWORK( N+1 ), LWORK-N,
787: $ IERR )
788: *
789: CALL ZGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, CWORK( N4+1 ),
790: $ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
791: $ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N,
792: $ IERR )
793: *
794: CALL ZGSVJ1( JOBV, N2, N2, N4, A, LDA, CWORK, SVA, MVL, V,
795: $ LDV, EPSLN, SFMIN, TOL, 1, CWORK( N+1 ),
796: $ LWORK-N, IERR )
797: *
798: CALL ZGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
799: $ CWORK( N2+1 ), SVA( N2+1 ), MVL,
800: $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
801: $ CWORK( N+1 ), LWORK-N, IERR )
802:
803: END IF
804: *
805: END IF
806: *
807: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
808: *
809: DO 1993 i = 1, NSWEEP
810: *
811: * .. go go go ...
812: *
813: MXAAPQ = ZERO
814: MXSINJ = ZERO
815: ISWROT = 0
816: *
817: NOTROT = 0
818: PSKIPPED = 0
819: *
820: * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
821: * 1 <= p < q <= N. This is the first step toward a blocked implementation
822: * of the rotations. New implementation, based on block transformations,
823: * is under development.
824: *
825: DO 2000 ibr = 1, NBL
826: *
827: igl = ( ibr-1 )*KBL + 1
828: *
1.4 bertrand 829: DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr )
1.1 bertrand 830: *
831: igl = igl + ir1*KBL
832: *
1.4 bertrand 833: DO 2001 p = igl, MIN( igl+KBL-1, N-1 )
1.1 bertrand 834: *
835: * .. de Rijk's pivoting
836: *
837: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
838: IF( p.NE.q ) THEN
839: CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
1.4 bertrand 840: IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1,
1.1 bertrand 841: $ V( 1, q ), 1 )
842: TEMP1 = SVA( p )
843: SVA( p ) = SVA( q )
844: SVA( q ) = TEMP1
845: AAPQ = CWORK(p)
846: CWORK(p) = CWORK(q)
847: CWORK(q) = AAPQ
848: END IF
849: *
850: IF( ir1.EQ.0 ) THEN
851: *
852: * Column norms are periodically updated by explicit
853: * norm computation.
854: *[!] Caveat:
855: * Unfortunately, some BLAS implementations compute DZNRM2(M,A(1,p),1)
856: * as SQRT(S=CDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to
857: * overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
858: * underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
859: * Hence, DZNRM2 cannot be trusted, not even in the case when
860: * the true norm is far from the under(over)flow boundaries.
861: * If properly implemented SCNRM2 is available, the IF-THEN-ELSE-END IF
862: * below should be replaced with "AAPP = DZNRM2( M, A(1,p), 1 )".
863: *
1.4 bertrand 864: IF( ( SVA( p ).LT.ROOTBIG ) .AND.
1.1 bertrand 865: $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN
866: SVA( p ) = DZNRM2( M, A( 1, p ), 1 )
867: ELSE
868: TEMP1 = ZERO
869: AAPP = ONE
870: CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
1.4 bertrand 871: SVA( p ) = TEMP1*SQRT( AAPP )
1.1 bertrand 872: END IF
873: AAPP = SVA( p )
874: ELSE
875: AAPP = SVA( p )
876: END IF
877: *
878: IF( AAPP.GT.ZERO ) THEN
879: *
880: PSKIPPED = 0
881: *
1.4 bertrand 882: DO 2002 q = p + 1, MIN( igl+KBL-1, N )
1.1 bertrand 883: *
884: AAQQ = SVA( q )
885: *
886: IF( AAQQ.GT.ZERO ) THEN
887: *
888: AAPP0 = AAPP
889: IF( AAQQ.GE.ONE ) THEN
890: ROTOK = ( SMALL*AAPP ).LE.AAQQ
891: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
1.4 bertrand 892: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
1.1 bertrand 893: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
894: ELSE
1.4 bertrand 895: CALL ZCOPY( M, A( 1, p ), 1,
1.1 bertrand 896: $ CWORK(N+1), 1 )
1.4 bertrand 897: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
1.1 bertrand 898: $ M, 1, CWORK(N+1), LDA, IERR )
899: AAPQ = ZDOTC( M, CWORK(N+1), 1,
900: $ A( 1, q ), 1 ) / AAQQ
901: END IF
902: ELSE
903: ROTOK = AAPP.LE.( AAQQ / SMALL )
904: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
1.4 bertrand 905: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
906: $ A( 1, q ), 1 ) / AAPP ) / AAQQ
1.1 bertrand 907: ELSE
1.4 bertrand 908: CALL ZCOPY( M, A( 1, q ), 1,
1.1 bertrand 909: $ CWORK(N+1), 1 )
910: CALL ZLASCL( 'G', 0, 0, AAQQ,
911: $ ONE, M, 1,
912: $ CWORK(N+1), LDA, IERR )
913: AAPQ = ZDOTC( M, A(1, p ), 1,
914: $ CWORK(N+1), 1 ) / AAPP
915: END IF
916: END IF
917: *
1.4 bertrand 918:
919: * AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
920: AAPQ1 = -ABS(AAPQ)
921: MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
1.1 bertrand 922: *
923: * TO rotate or NOT to rotate, THAT is the question ...
924: *
925: IF( ABS( AAPQ1 ).GT.TOL ) THEN
1.4 bertrand 926: OMPQ = AAPQ / ABS(AAPQ)
1.1 bertrand 927: *
928: * .. rotate
929: *[RTD] ROTATED = ROTATED + ONE
930: *
931: IF( ir1.EQ.0 ) THEN
932: NOTROT = 0
933: PSKIPPED = 0
934: ISWROT = ISWROT + 1
935: END IF
936: *
937: IF( ROTOK ) THEN
938: *
1.4 bertrand 939: AQOAP = AAQQ / AAPP
1.1 bertrand 940: APOAQ = AAPP / AAQQ
941: THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1
942: *
943: IF( ABS( THETA ).GT.BIGTHETA ) THEN
1.4 bertrand 944: *
1.1 bertrand 945: T = HALF / THETA
946: CS = ONE
947:
948: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1.4 bertrand 949: $ CS, CONJG(OMPQ)*T )
1.1 bertrand 950: IF ( RSVEC ) THEN
1.4 bertrand 951: CALL ZROT( MVL, V(1,p), 1,
952: $ V(1,q), 1, CS, CONJG(OMPQ)*T )
1.1 bertrand 953: END IF
1.4 bertrand 954:
955: SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1 bertrand 956: $ ONE+T*APOAQ*AAPQ1 ) )
1.4 bertrand 957: AAPP = AAPP*SQRT( MAX( ZERO,
1.1 bertrand 958: $ ONE-T*AQOAP*AAPQ1 ) )
1.4 bertrand 959: MXSINJ = MAX( MXSINJ, ABS( T ) )
1.1 bertrand 960: *
961: ELSE
962: *
963: * .. choose correct signum for THETA and rotate
964: *
1.4 bertrand 965: THSIGN = -SIGN( ONE, AAPQ1 )
966: T = ONE / ( THETA+THSIGN*
967: $ SQRT( ONE+THETA*THETA ) )
968: CS = SQRT( ONE / ( ONE+T*T ) )
1.1 bertrand 969: SN = T*CS
970: *
1.4 bertrand 971: MXSINJ = MAX( MXSINJ, ABS( SN ) )
972: SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1 bertrand 973: $ ONE+T*APOAQ*AAPQ1 ) )
1.4 bertrand 974: AAPP = AAPP*SQRT( MAX( ZERO,
1.1 bertrand 975: $ ONE-T*AQOAP*AAPQ1 ) )
976: *
977: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1.4 bertrand 978: $ CS, CONJG(OMPQ)*SN )
1.1 bertrand 979: IF ( RSVEC ) THEN
1.4 bertrand 980: CALL ZROT( MVL, V(1,p), 1,
981: $ V(1,q), 1, CS, CONJG(OMPQ)*SN )
982: END IF
983: END IF
984: CWORK(p) = -CWORK(q) * OMPQ
1.1 bertrand 985: *
986: ELSE
987: * .. have to use modified Gram-Schmidt like transformation
988: CALL ZCOPY( M, A( 1, p ), 1,
989: $ CWORK(N+1), 1 )
990: CALL ZLASCL( 'G', 0, 0, AAPP, ONE, M,
991: $ 1, CWORK(N+1), LDA,
992: $ IERR )
993: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, M,
994: $ 1, A( 1, q ), LDA, IERR )
995: CALL ZAXPY( M, -AAPQ, CWORK(N+1), 1,
996: $ A( 1, q ), 1 )
997: CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M,
998: $ 1, A( 1, q ), LDA, IERR )
1.4 bertrand 999: SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1 bertrand 1000: $ ONE-AAPQ1*AAPQ1 ) )
1.4 bertrand 1001: MXSINJ = MAX( MXSINJ, SFMIN )
1.1 bertrand 1002: END IF
1003: * END IF ROTOK THEN ... ELSE
1004: *
1005: * In the case of cancellation in updating SVA(q), SVA(p)
1006: * recompute SVA(q), SVA(p).
1007: *
1008: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1009: $ THEN
1010: IF( ( AAQQ.LT.ROOTBIG ) .AND.
1011: $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
1012: SVA( q ) = DZNRM2( M, A( 1, q ), 1 )
1013: ELSE
1014: T = ZERO
1015: AAQQ = ONE
1016: CALL ZLASSQ( M, A( 1, q ), 1, T,
1017: $ AAQQ )
1.4 bertrand 1018: SVA( q ) = T*SQRT( AAQQ )
1.1 bertrand 1019: END IF
1020: END IF
1021: IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
1022: IF( ( AAPP.LT.ROOTBIG ) .AND.
1023: $ ( AAPP.GT.ROOTSFMIN ) ) THEN
1024: AAPP = DZNRM2( M, A( 1, p ), 1 )
1025: ELSE
1026: T = ZERO
1027: AAPP = ONE
1028: CALL ZLASSQ( M, A( 1, p ), 1, T,
1029: $ AAPP )
1.4 bertrand 1030: AAPP = T*SQRT( AAPP )
1.1 bertrand 1031: END IF
1032: SVA( p ) = AAPP
1033: END IF
1034: *
1035: ELSE
1036: * A(:,p) and A(:,q) already numerically orthogonal
1037: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1038: *[RTD] SKIPPED = SKIPPED + 1
1039: PSKIPPED = PSKIPPED + 1
1040: END IF
1041: ELSE
1042: * A(:,q) is zero column
1043: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1044: PSKIPPED = PSKIPPED + 1
1045: END IF
1046: *
1047: IF( ( i.LE.SWBAND ) .AND.
1048: $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
1049: IF( ir1.EQ.0 )AAPP = -AAPP
1050: NOTROT = 0
1051: GO TO 2103
1052: END IF
1053: *
1054: 2002 CONTINUE
1055: * END q-LOOP
1056: *
1057: 2103 CONTINUE
1058: * bailed out of q-loop
1059: *
1060: SVA( p ) = AAPP
1061: *
1062: ELSE
1063: SVA( p ) = AAPP
1064: IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
1.4 bertrand 1065: $ NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p
1.1 bertrand 1066: END IF
1067: *
1068: 2001 CONTINUE
1069: * end of the p-loop
1070: * end of doing the block ( ibr, ibr )
1071: 1002 CONTINUE
1072: * end of ir1-loop
1073: *
1074: * ... go to the off diagonal blocks
1075: *
1076: igl = ( ibr-1 )*KBL + 1
1077: *
1078: DO 2010 jbc = ibr + 1, NBL
1079: *
1080: jgl = ( jbc-1 )*KBL + 1
1081: *
1082: * doing the block at ( ibr, jbc )
1083: *
1084: IJBLSK = 0
1.4 bertrand 1085: DO 2100 p = igl, MIN( igl+KBL-1, N )
1.1 bertrand 1086: *
1087: AAPP = SVA( p )
1088: IF( AAPP.GT.ZERO ) THEN
1089: *
1090: PSKIPPED = 0
1091: *
1.4 bertrand 1092: DO 2200 q = jgl, MIN( jgl+KBL-1, N )
1.1 bertrand 1093: *
1094: AAQQ = SVA( q )
1095: IF( AAQQ.GT.ZERO ) THEN
1096: AAPP0 = AAPP
1097: *
1098: * .. M x 2 Jacobi SVD ..
1099: *
1100: * Safe Gram matrix computation
1101: *
1102: IF( AAQQ.GE.ONE ) THEN
1103: IF( AAPP.GE.AAQQ ) THEN
1104: ROTOK = ( SMALL*AAPP ).LE.AAQQ
1105: ELSE
1106: ROTOK = ( SMALL*AAQQ ).LE.AAPP
1107: END IF
1108: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
1.4 bertrand 1109: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
1.1 bertrand 1110: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
1111: ELSE
1112: CALL ZCOPY( M, A( 1, p ), 1,
1113: $ CWORK(N+1), 1 )
1114: CALL ZLASCL( 'G', 0, 0, AAPP,
1115: $ ONE, M, 1,
1116: $ CWORK(N+1), LDA, IERR )
1117: AAPQ = ZDOTC( M, CWORK(N+1), 1,
1118: $ A( 1, q ), 1 ) / AAQQ
1119: END IF
1120: ELSE
1121: IF( AAPP.GE.AAQQ ) THEN
1122: ROTOK = AAPP.LE.( AAQQ / SMALL )
1123: ELSE
1124: ROTOK = AAQQ.LE.( AAPP / SMALL )
1125: END IF
1126: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
1.4 bertrand 1127: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
1128: $ A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
1129: $ / MIN(AAQQ,AAPP)
1.1 bertrand 1130: ELSE
1131: CALL ZCOPY( M, A( 1, q ), 1,
1132: $ CWORK(N+1), 1 )
1133: CALL ZLASCL( 'G', 0, 0, AAQQ,
1134: $ ONE, M, 1,
1135: $ CWORK(N+1), LDA, IERR )
1136: AAPQ = ZDOTC( M, A( 1, p ), 1,
1137: $ CWORK(N+1), 1 ) / AAPP
1138: END IF
1139: END IF
1140: *
1.4 bertrand 1141:
1142: * AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
1.1 bertrand 1143: AAPQ1 = -ABS(AAPQ)
1.4 bertrand 1144: MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
1.1 bertrand 1145: *
1146: * TO rotate or NOT to rotate, THAT is the question ...
1147: *
1148: IF( ABS( AAPQ1 ).GT.TOL ) THEN
1.4 bertrand 1149: OMPQ = AAPQ / ABS(AAPQ)
1.1 bertrand 1150: NOTROT = 0
1151: *[RTD] ROTATED = ROTATED + 1
1152: PSKIPPED = 0
1153: ISWROT = ISWROT + 1
1154: *
1155: IF( ROTOK ) THEN
1156: *
1157: AQOAP = AAQQ / AAPP
1158: APOAQ = AAPP / AAQQ
1159: THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1
1160: IF( AAQQ.GT.AAPP0 )THETA = -THETA
1161: *
1162: IF( ABS( THETA ).GT.BIGTHETA ) THEN
1163: T = HALF / THETA
1.4 bertrand 1164: CS = ONE
1.1 bertrand 1165: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1.4 bertrand 1166: $ CS, CONJG(OMPQ)*T )
1.1 bertrand 1167: IF( RSVEC ) THEN
1.4 bertrand 1168: CALL ZROT( MVL, V(1,p), 1,
1169: $ V(1,q), 1, CS, CONJG(OMPQ)*T )
1.1 bertrand 1170: END IF
1.4 bertrand 1171: SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1 bertrand 1172: $ ONE+T*APOAQ*AAPQ1 ) )
1.4 bertrand 1173: AAPP = AAPP*SQRT( MAX( ZERO,
1.1 bertrand 1174: $ ONE-T*AQOAP*AAPQ1 ) )
1.4 bertrand 1175: MXSINJ = MAX( MXSINJ, ABS( T ) )
1.1 bertrand 1176: ELSE
1177: *
1178: * .. choose correct signum for THETA and rotate
1179: *
1.4 bertrand 1180: THSIGN = -SIGN( ONE, AAPQ1 )
1.1 bertrand 1181: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
1182: T = ONE / ( THETA+THSIGN*
1.4 bertrand 1183: $ SQRT( ONE+THETA*THETA ) )
1184: CS = SQRT( ONE / ( ONE+T*T ) )
1.1 bertrand 1185: SN = T*CS
1.4 bertrand 1186: MXSINJ = MAX( MXSINJ, ABS( SN ) )
1187: SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1 bertrand 1188: $ ONE+T*APOAQ*AAPQ1 ) )
1.4 bertrand 1189: AAPP = AAPP*SQRT( MAX( ZERO,
1.1 bertrand 1190: $ ONE-T*AQOAP*AAPQ1 ) )
1191: *
1192: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1.4 bertrand 1193: $ CS, CONJG(OMPQ)*SN )
1.1 bertrand 1194: IF( RSVEC ) THEN
1.4 bertrand 1195: CALL ZROT( MVL, V(1,p), 1,
1196: $ V(1,q), 1, CS, CONJG(OMPQ)*SN )
1.1 bertrand 1197: END IF
1198: END IF
1.4 bertrand 1199: CWORK(p) = -CWORK(q) * OMPQ
1.1 bertrand 1200: *
1201: ELSE
1202: * .. have to use modified Gram-Schmidt like transformation
1203: IF( AAPP.GT.AAQQ ) THEN
1204: CALL ZCOPY( M, A( 1, p ), 1,
1205: $ CWORK(N+1), 1 )
1206: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
1207: $ M, 1, CWORK(N+1),LDA,
1208: $ IERR )
1209: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
1210: $ M, 1, A( 1, q ), LDA,
1211: $ IERR )
1212: CALL ZAXPY( M, -AAPQ, CWORK(N+1),
1213: $ 1, A( 1, q ), 1 )
1214: CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
1215: $ M, 1, A( 1, q ), LDA,
1216: $ IERR )
1.4 bertrand 1217: SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1 bertrand 1218: $ ONE-AAPQ1*AAPQ1 ) )
1.4 bertrand 1219: MXSINJ = MAX( MXSINJ, SFMIN )
1.1 bertrand 1220: ELSE
1221: CALL ZCOPY( M, A( 1, q ), 1,
1222: $ CWORK(N+1), 1 )
1223: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
1224: $ M, 1, CWORK(N+1),LDA,
1225: $ IERR )
1226: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
1227: $ M, 1, A( 1, p ), LDA,
1228: $ IERR )
1.4 bertrand 1229: CALL ZAXPY( M, -CONJG(AAPQ),
1.1 bertrand 1230: $ CWORK(N+1), 1, A( 1, p ), 1 )
1231: CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
1232: $ M, 1, A( 1, p ), LDA,
1233: $ IERR )
1.4 bertrand 1234: SVA( p ) = AAPP*SQRT( MAX( ZERO,
1.1 bertrand 1235: $ ONE-AAPQ1*AAPQ1 ) )
1.4 bertrand 1236: MXSINJ = MAX( MXSINJ, SFMIN )
1.1 bertrand 1237: END IF
1238: END IF
1239: * END IF ROTOK THEN ... ELSE
1240: *
1241: * In the case of cancellation in updating SVA(q), SVA(p)
1242: * .. recompute SVA(q), SVA(p)
1243: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1244: $ THEN
1245: IF( ( AAQQ.LT.ROOTBIG ) .AND.
1246: $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
1247: SVA( q ) = DZNRM2( M, A( 1, q ), 1)
1248: ELSE
1249: T = ZERO
1250: AAQQ = ONE
1251: CALL ZLASSQ( M, A( 1, q ), 1, T,
1252: $ AAQQ )
1.4 bertrand 1253: SVA( q ) = T*SQRT( AAQQ )
1.1 bertrand 1254: END IF
1255: END IF
1256: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
1257: IF( ( AAPP.LT.ROOTBIG ) .AND.
1258: $ ( AAPP.GT.ROOTSFMIN ) ) THEN
1259: AAPP = DZNRM2( M, A( 1, p ), 1 )
1260: ELSE
1261: T = ZERO
1262: AAPP = ONE
1263: CALL ZLASSQ( M, A( 1, p ), 1, T,
1264: $ AAPP )
1.4 bertrand 1265: AAPP = T*SQRT( AAPP )
1.1 bertrand 1266: END IF
1267: SVA( p ) = AAPP
1268: END IF
1269: * end of OK rotation
1270: ELSE
1271: NOTROT = NOTROT + 1
1272: *[RTD] SKIPPED = SKIPPED + 1
1273: PSKIPPED = PSKIPPED + 1
1274: IJBLSK = IJBLSK + 1
1275: END IF
1276: ELSE
1277: NOTROT = NOTROT + 1
1278: PSKIPPED = PSKIPPED + 1
1279: IJBLSK = IJBLSK + 1
1280: END IF
1281: *
1282: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
1283: $ THEN
1284: SVA( p ) = AAPP
1285: NOTROT = 0
1286: GO TO 2011
1287: END IF
1288: IF( ( i.LE.SWBAND ) .AND.
1289: $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
1290: AAPP = -AAPP
1291: NOTROT = 0
1292: GO TO 2203
1293: END IF
1294: *
1295: 2200 CONTINUE
1296: * end of the q-loop
1297: 2203 CONTINUE
1298: *
1299: SVA( p ) = AAPP
1300: *
1301: ELSE
1302: *
1303: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
1.4 bertrand 1304: $ MIN( jgl+KBL-1, N ) - jgl + 1
1.1 bertrand 1305: IF( AAPP.LT.ZERO )NOTROT = 0
1306: *
1307: END IF
1308: *
1309: 2100 CONTINUE
1310: * end of the p-loop
1311: 2010 CONTINUE
1312: * end of the jbc-loop
1313: 2011 CONTINUE
1314: *2011 bailed out of the jbc-loop
1.4 bertrand 1315: DO 2012 p = igl, MIN( igl+KBL-1, N )
1.1 bertrand 1316: SVA( p ) = ABS( SVA( p ) )
1317: 2012 CONTINUE
1318: ***
1319: 2000 CONTINUE
1320: *2000 :: end of the ibr-loop
1321: *
1322: * .. update SVA(N)
1323: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
1324: $ THEN
1325: SVA( N ) = DZNRM2( M, A( 1, N ), 1 )
1326: ELSE
1327: T = ZERO
1328: AAPP = ONE
1329: CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
1.4 bertrand 1330: SVA( N ) = T*SQRT( AAPP )
1.1 bertrand 1331: END IF
1332: *
1333: * Additional steering devices
1334: *
1335: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
1336: $ ( ISWROT.LE.N ) ) )SWBAND = i
1337: *
1.4 bertrand 1338: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( DBLE( N ) )*
1.2 bertrand 1339: $ TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1.1 bertrand 1340: GO TO 1994
1341: END IF
1342: *
1343: IF( NOTROT.GE.EMPTSW )GO TO 1994
1344: *
1345: 1993 CONTINUE
1346: * end i=1:NSWEEP loop
1347: *
1348: * #:( Reaching this point means that the procedure has not converged.
1349: INFO = NSWEEP - 1
1350: GO TO 1995
1351: *
1352: 1994 CONTINUE
1353: * #:) Reaching this point means numerical convergence after the i-th
1354: * sweep.
1355: *
1356: INFO = 0
1357: * #:) INFO = 0 confirms successful iterations.
1358: 1995 CONTINUE
1359: *
1360: * Sort the singular values and find how many are above
1361: * the underflow threshold.
1362: *
1363: N2 = 0
1364: N4 = 0
1365: DO 5991 p = 1, N - 1
1366: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
1367: IF( p.NE.q ) THEN
1368: TEMP1 = SVA( p )
1369: SVA( p ) = SVA( q )
1370: SVA( q ) = TEMP1
1371: CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
1372: IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
1373: END IF
1374: IF( SVA( p ).NE.ZERO ) THEN
1375: N4 = N4 + 1
1376: IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
1377: END IF
1378: 5991 CONTINUE
1379: IF( SVA( N ).NE.ZERO ) THEN
1380: N4 = N4 + 1
1381: IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
1382: END IF
1383: *
1384: * Normalize the left singular vectors.
1385: *
1386: IF( LSVEC .OR. UCTOL ) THEN
1.4 bertrand 1387: DO 1998 p = 1, N4
1388: * CALL ZDSCAL( M, ONE / SVA( p ), A( 1, p ), 1 )
1389: CALL ZLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M, IERR )
1.1 bertrand 1390: 1998 CONTINUE
1391: END IF
1392: *
1393: * Scale the product of Jacobi rotations.
1394: *
1395: IF( RSVEC ) THEN
1396: DO 2399 p = 1, N
1397: TEMP1 = ONE / DZNRM2( MVL, V( 1, p ), 1 )
1398: CALL ZDSCAL( MVL, TEMP1, V( 1, p ), 1 )
1399: 2399 CONTINUE
1400: END IF
1401: *
1402: * Undo scaling, if necessary (and possible).
1.4 bertrand 1403: IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL ) ) )
1.1 bertrand 1404: $ .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT.
1405: $ ( SFMIN / SKL ) ) ) ) THEN
1406: DO 2400 p = 1, N
1.4 bertrand 1407: SVA( p ) = SKL*SVA( p )
1.1 bertrand 1408: 2400 CONTINUE
1409: SKL = ONE
1410: END IF
1411: *
1412: RWORK( 1 ) = SKL
1413: * The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
1414: * then some of the singular values may overflow or underflow and
1415: * the spectrum is given in this factored representation.
1416: *
1.2 bertrand 1417: RWORK( 2 ) = DBLE( N4 )
1.1 bertrand 1418: * N4 is the number of computed nonzero singular values of A.
1419: *
1.2 bertrand 1420: RWORK( 3 ) = DBLE( N2 )
1.1 bertrand 1421: * N2 is the number of singular values of A greater than SFMIN.
1422: * If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
1423: * that may carry some information.
1424: *
1.2 bertrand 1425: RWORK( 4 ) = DBLE( i )
1.1 bertrand 1426: * i is the index of the last sweep before declaring convergence.
1427: *
1428: RWORK( 5 ) = MXAAPQ
1429: * MXAAPQ is the largest absolute value of scaled pivots in the
1430: * last sweep
1431: *
1432: RWORK( 6 ) = MXSINJ
1433: * MXSINJ is the largest absolute value of the sines of Jacobi angles
1434: * in the last sweep
1435: *
1436: RETURN
1437: * ..
1438: * .. END OF ZGESVJ
1439: * ..
1440: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgesvj.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack