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