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