Annotation of rpl/lapack/lapack/dbdsqr.f, revision 1.18
1.8 bertrand 1: *> \brief \b DBDSQR
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.14 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.14 bertrand 9: *> Download DBDSQR + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f">
1.8 bertrand 15: *> [TXT]</a>
1.14 bertrand 16: *> \endhtmlonly
1.8 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
22: * LDU, C, LDC, WORK, INFO )
1.14 bertrand 23: *
1.8 bertrand 24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
26: * INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
30: * $ VT( LDVT, * ), WORK( * )
31: * ..
1.14 bertrand 32: *
1.8 bertrand 33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DBDSQR computes the singular values and, optionally, the right and/or
40: *> left singular vectors from the singular value decomposition (SVD) of
41: *> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
42: *> zero-shift QR algorithm. The SVD of B has the form
1.14 bertrand 43: *>
1.8 bertrand 44: *> B = Q * S * P**T
1.14 bertrand 45: *>
1.8 bertrand 46: *> where S is the diagonal matrix of singular values, Q is an orthogonal
47: *> matrix of left singular vectors, and P is an orthogonal matrix of
48: *> right singular vectors. If left singular vectors are requested, this
49: *> subroutine actually returns U*Q instead of Q, and, if right singular
50: *> vectors are requested, this subroutine returns P**T*VT instead of
51: *> P**T, for given real input matrices U and VT. When U and VT are the
52: *> orthogonal matrices that reduce a general matrix A to bidiagonal
53: *> form: A = U*B*VT, as computed by DGEBRD, then
54: *>
55: *> A = (U*Q) * S * (P**T*VT)
56: *>
57: *> is the SVD of A. Optionally, the subroutine may also compute Q**T*C
58: *> for a given real input matrix C.
59: *>
60: *> See "Computing Small Singular Values of Bidiagonal Matrices With
61: *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
62: *> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
63: *> no. 5, pp. 873-912, Sept 1990) and
64: *> "Accurate singular values and differential qd algorithms," by
65: *> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
66: *> Department, University of California at Berkeley, July 1992
67: *> for a detailed description of the algorithm.
68: *> \endverbatim
69: *
70: * Arguments:
71: * ==========
72: *
73: *> \param[in] UPLO
74: *> \verbatim
75: *> UPLO is CHARACTER*1
76: *> = 'U': B is upper bidiagonal;
77: *> = 'L': B is lower bidiagonal.
78: *> \endverbatim
79: *>
80: *> \param[in] N
81: *> \verbatim
82: *> N is INTEGER
83: *> The order of the matrix B. N >= 0.
84: *> \endverbatim
85: *>
86: *> \param[in] NCVT
87: *> \verbatim
88: *> NCVT is INTEGER
89: *> The number of columns of the matrix VT. NCVT >= 0.
90: *> \endverbatim
91: *>
92: *> \param[in] NRU
93: *> \verbatim
94: *> NRU is INTEGER
95: *> The number of rows of the matrix U. NRU >= 0.
96: *> \endverbatim
97: *>
98: *> \param[in] NCC
99: *> \verbatim
100: *> NCC is INTEGER
101: *> The number of columns of the matrix C. NCC >= 0.
102: *> \endverbatim
103: *>
104: *> \param[in,out] D
105: *> \verbatim
106: *> D is DOUBLE PRECISION array, dimension (N)
107: *> On entry, the n diagonal elements of the bidiagonal matrix B.
108: *> On exit, if INFO=0, the singular values of B in decreasing
109: *> order.
110: *> \endverbatim
111: *>
112: *> \param[in,out] E
113: *> \verbatim
114: *> E is DOUBLE PRECISION array, dimension (N-1)
115: *> On entry, the N-1 offdiagonal elements of the bidiagonal
1.14 bertrand 116: *> matrix B.
1.8 bertrand 117: *> On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
118: *> will contain the diagonal and superdiagonal elements of a
119: *> bidiagonal matrix orthogonally equivalent to the one given
120: *> as input.
121: *> \endverbatim
122: *>
123: *> \param[in,out] VT
124: *> \verbatim
125: *> VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
126: *> On entry, an N-by-NCVT matrix VT.
127: *> On exit, VT is overwritten by P**T * VT.
128: *> Not referenced if NCVT = 0.
129: *> \endverbatim
130: *>
131: *> \param[in] LDVT
132: *> \verbatim
133: *> LDVT is INTEGER
134: *> The leading dimension of the array VT.
135: *> LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
136: *> \endverbatim
137: *>
138: *> \param[in,out] U
139: *> \verbatim
140: *> U is DOUBLE PRECISION array, dimension (LDU, N)
141: *> On entry, an NRU-by-N matrix U.
142: *> On exit, U is overwritten by U * Q.
143: *> Not referenced if NRU = 0.
144: *> \endverbatim
145: *>
146: *> \param[in] LDU
147: *> \verbatim
148: *> LDU is INTEGER
149: *> The leading dimension of the array U. LDU >= max(1,NRU).
150: *> \endverbatim
151: *>
152: *> \param[in,out] C
153: *> \verbatim
154: *> C is DOUBLE PRECISION array, dimension (LDC, NCC)
155: *> On entry, an N-by-NCC matrix C.
156: *> On exit, C is overwritten by Q**T * C.
157: *> Not referenced if NCC = 0.
158: *> \endverbatim
159: *>
160: *> \param[in] LDC
161: *> \verbatim
162: *> LDC is INTEGER
163: *> The leading dimension of the array C.
164: *> LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
165: *> \endverbatim
166: *>
167: *> \param[out] WORK
168: *> \verbatim
1.18 ! bertrand 169: *> WORK is DOUBLE PRECISION array, dimension (4*(N-1))
1.8 bertrand 170: *> \endverbatim
171: *>
172: *> \param[out] INFO
173: *> \verbatim
174: *> INFO is INTEGER
175: *> = 0: successful exit
176: *> < 0: If INFO = -i, the i-th argument had an illegal value
177: *> > 0:
178: *> if NCVT = NRU = NCC = 0,
179: *> = 1, a split was marked by a positive value in E
180: *> = 2, current block of Z not diagonalized after 30*N
181: *> iterations (in inner while loop)
1.14 bertrand 182: *> = 3, termination criterion of outer while loop not met
1.8 bertrand 183: *> (program created more than N unreduced blocks)
184: *> else NCVT = NRU = NCC = 0,
185: *> the algorithm did not converge; D and E contain the
186: *> elements of a bidiagonal matrix which is orthogonally
187: *> similar to the input matrix B; if INFO = i, i
188: *> elements of E have not converged to zero.
189: *> \endverbatim
190: *
191: *> \par Internal Parameters:
192: * =========================
193: *>
194: *> \verbatim
195: *> TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
196: *> TOLMUL controls the convergence criterion of the QR loop.
197: *> If it is positive, TOLMUL*EPS is the desired relative
198: *> precision in the computed singular values.
199: *> If it is negative, abs(TOLMUL*EPS*sigma_max) is the
200: *> desired absolute accuracy in the computed singular
201: *> values (corresponds to relative accuracy
202: *> abs(TOLMUL*EPS) in the largest singular value.
203: *> abs(TOLMUL) should be between 1 and 1/EPS, and preferably
204: *> between 10 (for fast convergence) and .1/EPS
205: *> (for there to be some accuracy in the results).
206: *> Default is to lose at either one eighth or 2 of the
207: *> available decimal digits in each computed singular value
208: *> (whichever is smaller).
209: *>
210: *> MAXITR INTEGER, default = 6
211: *> MAXITR controls the maximum number of passes of the
212: *> algorithm through its inner loop. The algorithms stops
213: *> (and so fails to converge) if the number of passes
214: *> through the inner loop exceeds MAXITR*N**2.
1.16 bertrand 215: *>
216: *> \endverbatim
217: *
218: *> \par Note:
219: * ===========
220: *>
221: *> \verbatim
222: *> Bug report from Cezary Dendek.
223: *> On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
224: *> removed since it can overflow pretty easily (for N larger or equal
225: *> than 18,919). We instead use MAXITDIVN = MAXITR*N.
1.8 bertrand 226: *> \endverbatim
227: *
228: * Authors:
229: * ========
230: *
1.14 bertrand 231: *> \author Univ. of Tennessee
232: *> \author Univ. of California Berkeley
233: *> \author Univ. of Colorado Denver
234: *> \author NAG Ltd.
1.8 bertrand 235: *
1.16 bertrand 236: *> \date June 2017
1.8 bertrand 237: *
238: *> \ingroup auxOTHERcomputational
239: *
240: * =====================================================================
1.1 bertrand 241: SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
242: $ LDU, C, LDC, WORK, INFO )
243: *
1.16 bertrand 244: * -- LAPACK computational routine (version 3.7.1) --
1.1 bertrand 245: * -- LAPACK is a software package provided by Univ. of Tennessee, --
246: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16 bertrand 247: * June 2017
1.1 bertrand 248: *
249: * .. Scalar Arguments ..
250: CHARACTER UPLO
251: INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
252: * ..
253: * .. Array Arguments ..
254: DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
255: $ VT( LDVT, * ), WORK( * )
256: * ..
257: *
258: * =====================================================================
259: *
260: * .. Parameters ..
261: DOUBLE PRECISION ZERO
262: PARAMETER ( ZERO = 0.0D0 )
263: DOUBLE PRECISION ONE
264: PARAMETER ( ONE = 1.0D0 )
265: DOUBLE PRECISION NEGONE
266: PARAMETER ( NEGONE = -1.0D0 )
267: DOUBLE PRECISION HNDRTH
268: PARAMETER ( HNDRTH = 0.01D0 )
269: DOUBLE PRECISION TEN
270: PARAMETER ( TEN = 10.0D0 )
271: DOUBLE PRECISION HNDRD
272: PARAMETER ( HNDRD = 100.0D0 )
273: DOUBLE PRECISION MEIGTH
274: PARAMETER ( MEIGTH = -0.125D0 )
275: INTEGER MAXITR
276: PARAMETER ( MAXITR = 6 )
277: * ..
278: * .. Local Scalars ..
279: LOGICAL LOWER, ROTATE
1.16 bertrand 280: INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
281: $ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM
1.1 bertrand 282: DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
283: $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
284: $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
285: $ SN, THRESH, TOL, TOLMUL, UNFL
286: * ..
287: * .. External Functions ..
288: LOGICAL LSAME
289: DOUBLE PRECISION DLAMCH
290: EXTERNAL LSAME, DLAMCH
291: * ..
292: * .. External Subroutines ..
293: EXTERNAL DLARTG, DLAS2, DLASQ1, DLASR, DLASV2, DROT,
294: $ DSCAL, DSWAP, XERBLA
295: * ..
296: * .. Intrinsic Functions ..
297: INTRINSIC ABS, DBLE, MAX, MIN, SIGN, SQRT
298: * ..
299: * .. Executable Statements ..
300: *
301: * Test the input parameters.
302: *
303: INFO = 0
304: LOWER = LSAME( UPLO, 'L' )
305: IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN
306: INFO = -1
307: ELSE IF( N.LT.0 ) THEN
308: INFO = -2
309: ELSE IF( NCVT.LT.0 ) THEN
310: INFO = -3
311: ELSE IF( NRU.LT.0 ) THEN
312: INFO = -4
313: ELSE IF( NCC.LT.0 ) THEN
314: INFO = -5
315: ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
316: $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
317: INFO = -9
318: ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
319: INFO = -11
320: ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
321: $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
322: INFO = -13
323: END IF
324: IF( INFO.NE.0 ) THEN
325: CALL XERBLA( 'DBDSQR', -INFO )
326: RETURN
327: END IF
328: IF( N.EQ.0 )
329: $ RETURN
330: IF( N.EQ.1 )
331: $ GO TO 160
332: *
333: * ROTATE is true if any singular vectors desired, false otherwise
334: *
335: ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
336: *
337: * If no singular vectors desired, use qd algorithm
338: *
339: IF( .NOT.ROTATE ) THEN
340: CALL DLASQ1( N, D, E, WORK, INFO )
1.8 bertrand 341: *
342: * If INFO equals 2, dqds didn't finish, try to finish
1.14 bertrand 343: *
1.8 bertrand 344: IF( INFO .NE. 2 ) RETURN
345: INFO = 0
1.1 bertrand 346: END IF
347: *
348: NM1 = N - 1
349: NM12 = NM1 + NM1
350: NM13 = NM12 + NM1
351: IDIR = 0
352: *
353: * Get machine constants
354: *
355: EPS = DLAMCH( 'Epsilon' )
356: UNFL = DLAMCH( 'Safe minimum' )
357: *
358: * If matrix lower bidiagonal, rotate to be upper bidiagonal
359: * by applying Givens rotations on the left
360: *
361: IF( LOWER ) THEN
362: DO 10 I = 1, N - 1
363: CALL DLARTG( D( I ), E( I ), CS, SN, R )
364: D( I ) = R
365: E( I ) = SN*D( I+1 )
366: D( I+1 ) = CS*D( I+1 )
367: WORK( I ) = CS
368: WORK( NM1+I ) = SN
369: 10 CONTINUE
370: *
371: * Update singular vectors if desired
372: *
373: IF( NRU.GT.0 )
374: $ CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), WORK( N ), U,
375: $ LDU )
376: IF( NCC.GT.0 )
377: $ CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), WORK( N ), C,
378: $ LDC )
379: END IF
380: *
381: * Compute singular values to relative accuracy TOL
382: * (By setting TOL to be negative, algorithm will compute
383: * singular values to absolute accuracy ABS(TOL)*norm(input matrix))
384: *
385: TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
386: TOL = TOLMUL*EPS
387: *
388: * Compute approximate maximum, minimum singular values
389: *
390: SMAX = ZERO
391: DO 20 I = 1, N
392: SMAX = MAX( SMAX, ABS( D( I ) ) )
393: 20 CONTINUE
394: DO 30 I = 1, N - 1
395: SMAX = MAX( SMAX, ABS( E( I ) ) )
396: 30 CONTINUE
397: SMINL = ZERO
398: IF( TOL.GE.ZERO ) THEN
399: *
400: * Relative accuracy desired
401: *
402: SMINOA = ABS( D( 1 ) )
403: IF( SMINOA.EQ.ZERO )
404: $ GO TO 50
405: MU = SMINOA
406: DO 40 I = 2, N
407: MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
408: SMINOA = MIN( SMINOA, MU )
409: IF( SMINOA.EQ.ZERO )
410: $ GO TO 50
411: 40 CONTINUE
412: 50 CONTINUE
413: SMINOA = SMINOA / SQRT( DBLE( N ) )
1.16 bertrand 414: THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) )
1.1 bertrand 415: ELSE
416: *
417: * Absolute accuracy desired
418: *
1.16 bertrand 419: THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) )
1.1 bertrand 420: END IF
421: *
422: * Prepare for main iteration loop for the singular values
423: * (MAXIT is the maximum number of passes through the inner
424: * loop permitted before nonconvergence signalled.)
425: *
1.16 bertrand 426: MAXITDIVN = MAXITR*N
427: ITERDIVN = 0
428: ITER = -1
1.1 bertrand 429: OLDLL = -1
430: OLDM = -1
431: *
432: * M points to last element of unconverged part of matrix
433: *
434: M = N
435: *
436: * Begin main iteration loop
437: *
438: 60 CONTINUE
439: *
440: * Check for convergence or exceeding iteration count
441: *
442: IF( M.LE.1 )
443: $ GO TO 160
1.16 bertrand 444: *
445: IF( ITER.GE.N ) THEN
446: ITER = ITER - N
447: ITERDIVN = ITERDIVN + 1
448: IF( ITERDIVN.GE.MAXITDIVN )
449: $ GO TO 200
450: END IF
1.1 bertrand 451: *
452: * Find diagonal block of matrix to work on
453: *
454: IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
455: $ D( M ) = ZERO
456: SMAX = ABS( D( M ) )
457: SMIN = SMAX
458: DO 70 LLL = 1, M - 1
459: LL = M - LLL
460: ABSS = ABS( D( LL ) )
461: ABSE = ABS( E( LL ) )
462: IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
463: $ D( LL ) = ZERO
464: IF( ABSE.LE.THRESH )
465: $ GO TO 80
466: SMIN = MIN( SMIN, ABSS )
467: SMAX = MAX( SMAX, ABSS, ABSE )
468: 70 CONTINUE
469: LL = 0
470: GO TO 90
471: 80 CONTINUE
472: E( LL ) = ZERO
473: *
474: * Matrix splits since E(LL) = 0
475: *
476: IF( LL.EQ.M-1 ) THEN
477: *
478: * Convergence of bottom singular value, return to top of loop
479: *
480: M = M - 1
481: GO TO 60
482: END IF
483: 90 CONTINUE
484: LL = LL + 1
485: *
486: * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
487: *
488: IF( LL.EQ.M-1 ) THEN
489: *
490: * 2 by 2 block, handle separately
491: *
492: CALL DLASV2( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
493: $ COSR, SINL, COSL )
494: D( M-1 ) = SIGMX
495: E( M-1 ) = ZERO
496: D( M ) = SIGMN
497: *
498: * Compute singular vectors, if desired
499: *
500: IF( NCVT.GT.0 )
501: $ CALL DROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,
502: $ SINR )
503: IF( NRU.GT.0 )
504: $ CALL DROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
505: IF( NCC.GT.0 )
506: $ CALL DROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
507: $ SINL )
508: M = M - 2
509: GO TO 60
510: END IF
511: *
512: * If working on new submatrix, choose shift direction
513: * (from larger end diagonal element towards smaller)
514: *
515: IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
516: IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
517: *
518: * Chase bulge from top (big end) to bottom (small end)
519: *
520: IDIR = 1
521: ELSE
522: *
523: * Chase bulge from bottom (big end) to top (small end)
524: *
525: IDIR = 2
526: END IF
527: END IF
528: *
529: * Apply convergence tests
530: *
531: IF( IDIR.EQ.1 ) THEN
532: *
533: * Run convergence test in forward direction
534: * First apply standard test to bottom of matrix
535: *
536: IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
537: $ ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
538: E( M-1 ) = ZERO
539: GO TO 60
540: END IF
541: *
542: IF( TOL.GE.ZERO ) THEN
543: *
544: * If relative accuracy desired,
545: * apply convergence criterion forward
546: *
547: MU = ABS( D( LL ) )
548: SMINL = MU
549: DO 100 LLL = LL, M - 1
550: IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
551: E( LLL ) = ZERO
552: GO TO 60
553: END IF
554: MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
555: SMINL = MIN( SMINL, MU )
556: 100 CONTINUE
557: END IF
558: *
559: ELSE
560: *
561: * Run convergence test in backward direction
562: * First apply standard test to top of matrix
563: *
564: IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
565: $ ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
566: E( LL ) = ZERO
567: GO TO 60
568: END IF
569: *
570: IF( TOL.GE.ZERO ) THEN
571: *
572: * If relative accuracy desired,
573: * apply convergence criterion backward
574: *
575: MU = ABS( D( M ) )
576: SMINL = MU
577: DO 110 LLL = M - 1, LL, -1
578: IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
579: E( LLL ) = ZERO
580: GO TO 60
581: END IF
582: MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
583: SMINL = MIN( SMINL, MU )
584: 110 CONTINUE
585: END IF
586: END IF
587: OLDLL = LL
588: OLDM = M
589: *
590: * Compute shift. First, test if shifting would ruin relative
591: * accuracy, and if so set the shift to zero.
592: *
593: IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
594: $ MAX( EPS, HNDRTH*TOL ) ) THEN
595: *
596: * Use a zero shift to avoid loss of relative accuracy
597: *
598: SHIFT = ZERO
599: ELSE
600: *
601: * Compute the shift from 2-by-2 block at end of matrix
602: *
603: IF( IDIR.EQ.1 ) THEN
604: SLL = ABS( D( LL ) )
605: CALL DLAS2( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
606: ELSE
607: SLL = ABS( D( M ) )
608: CALL DLAS2( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
609: END IF
610: *
611: * Test if shift negligible, and if so set to zero
612: *
613: IF( SLL.GT.ZERO ) THEN
614: IF( ( SHIFT / SLL )**2.LT.EPS )
615: $ SHIFT = ZERO
616: END IF
617: END IF
618: *
619: * Increment iteration count
620: *
621: ITER = ITER + M - LL
622: *
623: * If SHIFT = 0, do simplified QR iteration
624: *
625: IF( SHIFT.EQ.ZERO ) THEN
626: IF( IDIR.EQ.1 ) THEN
627: *
628: * Chase bulge from top to bottom
629: * Save cosines and sines for later singular vector updates
630: *
631: CS = ONE
632: OLDCS = ONE
633: DO 120 I = LL, M - 1
634: CALL DLARTG( D( I )*CS, E( I ), CS, SN, R )
635: IF( I.GT.LL )
636: $ E( I-1 ) = OLDSN*R
637: CALL DLARTG( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
638: WORK( I-LL+1 ) = CS
639: WORK( I-LL+1+NM1 ) = SN
640: WORK( I-LL+1+NM12 ) = OLDCS
641: WORK( I-LL+1+NM13 ) = OLDSN
642: 120 CONTINUE
643: H = D( M )*CS
644: D( M ) = H*OLDCS
645: E( M-1 ) = H*OLDSN
646: *
647: * Update singular vectors
648: *
649: IF( NCVT.GT.0 )
650: $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
651: $ WORK( N ), VT( LL, 1 ), LDVT )
652: IF( NRU.GT.0 )
653: $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
654: $ WORK( NM13+1 ), U( 1, LL ), LDU )
655: IF( NCC.GT.0 )
656: $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
657: $ WORK( NM13+1 ), C( LL, 1 ), LDC )
658: *
659: * Test convergence
660: *
661: IF( ABS( E( M-1 ) ).LE.THRESH )
662: $ E( M-1 ) = ZERO
663: *
664: ELSE
665: *
666: * Chase bulge from bottom to top
667: * Save cosines and sines for later singular vector updates
668: *
669: CS = ONE
670: OLDCS = ONE
671: DO 130 I = M, LL + 1, -1
672: CALL DLARTG( D( I )*CS, E( I-1 ), CS, SN, R )
673: IF( I.LT.M )
674: $ E( I ) = OLDSN*R
675: CALL DLARTG( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
676: WORK( I-LL ) = CS
677: WORK( I-LL+NM1 ) = -SN
678: WORK( I-LL+NM12 ) = OLDCS
679: WORK( I-LL+NM13 ) = -OLDSN
680: 130 CONTINUE
681: H = D( LL )*CS
682: D( LL ) = H*OLDCS
683: E( LL ) = H*OLDSN
684: *
685: * Update singular vectors
686: *
687: IF( NCVT.GT.0 )
688: $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
689: $ WORK( NM13+1 ), VT( LL, 1 ), LDVT )
690: IF( NRU.GT.0 )
691: $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
692: $ WORK( N ), U( 1, LL ), LDU )
693: IF( NCC.GT.0 )
694: $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
695: $ WORK( N ), C( LL, 1 ), LDC )
696: *
697: * Test convergence
698: *
699: IF( ABS( E( LL ) ).LE.THRESH )
700: $ E( LL ) = ZERO
701: END IF
702: ELSE
703: *
704: * Use nonzero shift
705: *
706: IF( IDIR.EQ.1 ) THEN
707: *
708: * Chase bulge from top to bottom
709: * Save cosines and sines for later singular vector updates
710: *
711: F = ( ABS( D( LL ) )-SHIFT )*
712: $ ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
713: G = E( LL )
714: DO 140 I = LL, M - 1
715: CALL DLARTG( F, G, COSR, SINR, R )
716: IF( I.GT.LL )
717: $ E( I-1 ) = R
718: F = COSR*D( I ) + SINR*E( I )
719: E( I ) = COSR*E( I ) - SINR*D( I )
720: G = SINR*D( I+1 )
721: D( I+1 ) = COSR*D( I+1 )
722: CALL DLARTG( F, G, COSL, SINL, R )
723: D( I ) = R
724: F = COSL*E( I ) + SINL*D( I+1 )
725: D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
726: IF( I.LT.M-1 ) THEN
727: G = SINL*E( I+1 )
728: E( I+1 ) = COSL*E( I+1 )
729: END IF
730: WORK( I-LL+1 ) = COSR
731: WORK( I-LL+1+NM1 ) = SINR
732: WORK( I-LL+1+NM12 ) = COSL
733: WORK( I-LL+1+NM13 ) = SINL
734: 140 CONTINUE
735: E( M-1 ) = F
736: *
737: * Update singular vectors
738: *
739: IF( NCVT.GT.0 )
740: $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCVT, WORK( 1 ),
741: $ WORK( N ), VT( LL, 1 ), LDVT )
742: IF( NRU.GT.0 )
743: $ CALL DLASR( 'R', 'V', 'F', NRU, M-LL+1, WORK( NM12+1 ),
744: $ WORK( NM13+1 ), U( 1, LL ), LDU )
745: IF( NCC.GT.0 )
746: $ CALL DLASR( 'L', 'V', 'F', M-LL+1, NCC, WORK( NM12+1 ),
747: $ WORK( NM13+1 ), C( LL, 1 ), LDC )
748: *
749: * Test convergence
750: *
751: IF( ABS( E( M-1 ) ).LE.THRESH )
752: $ E( M-1 ) = ZERO
753: *
754: ELSE
755: *
756: * Chase bulge from bottom to top
757: * Save cosines and sines for later singular vector updates
758: *
759: F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
760: $ D( M ) )
761: G = E( M-1 )
762: DO 150 I = M, LL + 1, -1
763: CALL DLARTG( F, G, COSR, SINR, R )
764: IF( I.LT.M )
765: $ E( I ) = R
766: F = COSR*D( I ) + SINR*E( I-1 )
767: E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
768: G = SINR*D( I-1 )
769: D( I-1 ) = COSR*D( I-1 )
770: CALL DLARTG( F, G, COSL, SINL, R )
771: D( I ) = R
772: F = COSL*E( I-1 ) + SINL*D( I-1 )
773: D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
774: IF( I.GT.LL+1 ) THEN
775: G = SINL*E( I-2 )
776: E( I-2 ) = COSL*E( I-2 )
777: END IF
778: WORK( I-LL ) = COSR
779: WORK( I-LL+NM1 ) = -SINR
780: WORK( I-LL+NM12 ) = COSL
781: WORK( I-LL+NM13 ) = -SINL
782: 150 CONTINUE
783: E( LL ) = F
784: *
785: * Test convergence
786: *
787: IF( ABS( E( LL ) ).LE.THRESH )
788: $ E( LL ) = ZERO
789: *
790: * Update singular vectors if desired
791: *
792: IF( NCVT.GT.0 )
793: $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCVT, WORK( NM12+1 ),
794: $ WORK( NM13+1 ), VT( LL, 1 ), LDVT )
795: IF( NRU.GT.0 )
796: $ CALL DLASR( 'R', 'V', 'B', NRU, M-LL+1, WORK( 1 ),
797: $ WORK( N ), U( 1, LL ), LDU )
798: IF( NCC.GT.0 )
799: $ CALL DLASR( 'L', 'V', 'B', M-LL+1, NCC, WORK( 1 ),
800: $ WORK( N ), C( LL, 1 ), LDC )
801: END IF
802: END IF
803: *
804: * QR iteration finished, go back and check convergence
805: *
806: GO TO 60
807: *
808: * All singular values converged, so make them positive
809: *
810: 160 CONTINUE
811: DO 170 I = 1, N
812: IF( D( I ).LT.ZERO ) THEN
813: D( I ) = -D( I )
814: *
815: * Change sign of singular vectors, if desired
816: *
817: IF( NCVT.GT.0 )
818: $ CALL DSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
819: END IF
820: 170 CONTINUE
821: *
822: * Sort the singular values into decreasing order (insertion sort on
823: * singular values, but only one transposition per singular vector)
824: *
825: DO 190 I = 1, N - 1
826: *
827: * Scan for smallest D(I)
828: *
829: ISUB = 1
830: SMIN = D( 1 )
831: DO 180 J = 2, N + 1 - I
832: IF( D( J ).LE.SMIN ) THEN
833: ISUB = J
834: SMIN = D( J )
835: END IF
836: 180 CONTINUE
837: IF( ISUB.NE.N+1-I ) THEN
838: *
839: * Swap singular values and vectors
840: *
841: D( ISUB ) = D( N+1-I )
842: D( N+1-I ) = SMIN
843: IF( NCVT.GT.0 )
844: $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
845: $ LDVT )
846: IF( NRU.GT.0 )
847: $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
848: IF( NCC.GT.0 )
849: $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
850: END IF
851: 190 CONTINUE
852: GO TO 220
853: *
854: * Maximum number of iterations exceeded, failure to converge
855: *
856: 200 CONTINUE
857: INFO = 0
858: DO 210 I = 1, N - 1
859: IF( E( I ).NE.ZERO )
860: $ INFO = INFO + 1
861: 210 CONTINUE
862: 220 CONTINUE
863: RETURN
864: *
865: * End of DBDSQR
866: *
867: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dbdsqr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack