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