Annotation of rpl/lapack/lapack/zbdsqr.f, revision 1.1
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