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