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