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