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