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