Annotation of rpl/lapack/lapack/zbbcsd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
! 2: $ THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T,
! 3: $ V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E,
! 4: $ B22D, B22E, RWORK, LRWORK, INFO )
! 5: IMPLICIT NONE
! 6: *
! 7: * -- LAPACK routine (version 3.3.0) --
! 8: *
! 9: * -- Contributed by Brian Sutton of the Randolph-Macon College --
! 10: * -- November 2010
! 11: *
! 12: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 13: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 14: *
! 15: * .. Scalar Arguments ..
! 16: CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
! 17: INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
! 18: * ..
! 19: * .. Array Arguments ..
! 20: DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
! 21: $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
! 22: $ PHI( * ), THETA( * ), RWORK( * )
! 23: COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
! 24: $ V2T( LDV2T, * )
! 25: * ..
! 26: *
! 27: * Purpose
! 28: * =======
! 29: *
! 30: * ZBBCSD computes the CS decomposition of a unitary matrix in
! 31: * bidiagonal-block form,
! 32: *
! 33: *
! 34: * [ B11 | B12 0 0 ]
! 35: * [ 0 | 0 -I 0 ]
! 36: * X = [----------------]
! 37: * [ B21 | B22 0 0 ]
! 38: * [ 0 | 0 0 I ]
! 39: *
! 40: * [ C | -S 0 0 ]
! 41: * [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H
! 42: * = [---------] [---------------] [---------] .
! 43: * [ | U2 ] [ S | C 0 0 ] [ | V2 ]
! 44: * [ 0 | 0 0 I ]
! 45: *
! 46: * X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
! 47: * than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
! 48: * transposed and/or permuted. This can be done in constant time using
! 49: * the TRANS and SIGNS options. See ZUNCSD for details.)
! 50: *
! 51: * The bidiagonal matrices B11, B12, B21, and B22 are represented
! 52: * implicitly by angles THETA(1:Q) and PHI(1:Q-1).
! 53: *
! 54: * The unitary matrices U1, U2, V1T, and V2T are input/output.
! 55: * The input matrices are pre- or post-multiplied by the appropriate
! 56: * singular vector matrices.
! 57: *
! 58: * Arguments
! 59: * =========
! 60: *
! 61: * JOBU1 (input) CHARACTER
! 62: * = 'Y': U1 is updated;
! 63: * otherwise: U1 is not updated.
! 64: *
! 65: * JOBU2 (input) CHARACTER
! 66: * = 'Y': U2 is updated;
! 67: * otherwise: U2 is not updated.
! 68: *
! 69: * JOBV1T (input) CHARACTER
! 70: * = 'Y': V1T is updated;
! 71: * otherwise: V1T is not updated.
! 72: *
! 73: * JOBV2T (input) CHARACTER
! 74: * = 'Y': V2T is updated;
! 75: * otherwise: V2T is not updated.
! 76: *
! 77: * TRANS (input) CHARACTER
! 78: * = 'T': X, U1, U2, V1T, and V2T are stored in row-major
! 79: * order;
! 80: * otherwise: X, U1, U2, V1T, and V2T are stored in column-
! 81: * major order.
! 82: *
! 83: * M (input) INTEGER
! 84: * The number of rows and columns in X, the unitary matrix in
! 85: * bidiagonal-block form.
! 86: *
! 87: * P (input) INTEGER
! 88: * The number of rows in the top-left block of X. 0 <= P <= M.
! 89: *
! 90: * Q (input) INTEGER
! 91: * The number of columns in the top-left block of X.
! 92: * 0 <= Q <= MIN(P,M-P,M-Q).
! 93: *
! 94: * THETA (input/output) DOUBLE PRECISION array, dimension (Q)
! 95: * On entry, the angles THETA(1),...,THETA(Q) that, along with
! 96: * PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
! 97: * form. On exit, the angles whose cosines and sines define the
! 98: * diagonal blocks in the CS decomposition.
! 99: *
! 100: * PHI (input/workspace) DOUBLE PRECISION array, dimension (Q-1)
! 101: * The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
! 102: * THETA(Q), define the matrix in bidiagonal-block form.
! 103: *
! 104: * U1 (input/output) COMPLEX*16 array, dimension (LDU1,P)
! 105: * On entry, an LDU1-by-P matrix. On exit, U1 is postmultiplied
! 106: * by the left singular vector matrix common to [ B11 ; 0 ] and
! 107: * [ B12 0 0 ; 0 -I 0 0 ].
! 108: *
! 109: * LDU1 (input) INTEGER
! 110: * The leading dimension of the array U1.
! 111: *
! 112: * U2 (input/output) COMPLEX*16 array, dimension (LDU2,M-P)
! 113: * On entry, an LDU2-by-(M-P) matrix. On exit, U2 is
! 114: * postmultiplied by the left singular vector matrix common to
! 115: * [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
! 116: *
! 117: * LDU2 (input) INTEGER
! 118: * The leading dimension of the array U2.
! 119: *
! 120: * V1T (input/output) COMPLEX*16 array, dimension (LDV1T,Q)
! 121: * On entry, a LDV1T-by-Q matrix. On exit, V1T is premultiplied
! 122: * by the conjugate transpose of the right singular vector
! 123: * matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
! 124: *
! 125: * LDV1T (input) INTEGER
! 126: * The leading dimension of the array V1T.
! 127: *
! 128: * V2T (input/output) COMPLEX*16 array, dimenison (LDV2T,M-Q)
! 129: * On entry, a LDV2T-by-(M-Q) matrix. On exit, V2T is
! 130: * premultiplied by the conjugate transpose of the right
! 131: * singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
! 132: * [ B22 0 0 ; 0 0 I ].
! 133: *
! 134: * LDV2T (input) INTEGER
! 135: * The leading dimension of the array V2T.
! 136: *
! 137: * B11D (output) DOUBLE PRECISION array, dimension (Q)
! 138: * When ZBBCSD converges, B11D contains the cosines of THETA(1),
! 139: * ..., THETA(Q). If ZBBCSD fails to converge, then B11D
! 140: * contains the diagonal of the partially reduced top-left
! 141: * block.
! 142: *
! 143: * B11E (output) DOUBLE PRECISION array, dimension (Q-1)
! 144: * When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails
! 145: * to converge, then B11E contains the superdiagonal of the
! 146: * partially reduced top-left block.
! 147: *
! 148: * B12D (output) DOUBLE PRECISION array, dimension (Q)
! 149: * When ZBBCSD converges, B12D contains the negative sines of
! 150: * THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
! 151: * B12D contains the diagonal of the partially reduced top-right
! 152: * block.
! 153: *
! 154: * B12E (output) DOUBLE PRECISION array, dimension (Q-1)
! 155: * When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails
! 156: * to converge, then B12E contains the subdiagonal of the
! 157: * partially reduced top-right block.
! 158: *
! 159: * RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 160: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 161: *
! 162: * LRWORK (input) INTEGER
! 163: * The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).
! 164: *
! 165: * If LRWORK = -1, then a workspace query is assumed; the
! 166: * routine only calculates the optimal size of the RWORK array,
! 167: * returns this value as the first entry of the work array, and
! 168: * no error message related to LRWORK is issued by XERBLA.
! 169: *
! 170: * INFO (output) INTEGER
! 171: * = 0: successful exit.
! 172: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 173: * > 0: if ZBBCSD did not converge, INFO specifies the number
! 174: * of nonzero entries in PHI, and B11D, B11E, etc.,
! 175: * contain the partially reduced matrix.
! 176: *
! 177: * Reference
! 178: * =========
! 179: *
! 180: * [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
! 181: * Algorithms, 50(1):33-65, 2009.
! 182: *
! 183: * Internal Parameters
! 184: * ===================
! 185: *
! 186: * TOLMUL DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))
! 187: * TOLMUL controls the convergence criterion of the QR loop.
! 188: * Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
! 189: * are within TOLMUL*EPS of either bound.
! 190: *
! 191: * ===================================================================
! 192: *
! 193: * .. Parameters ..
! 194: INTEGER MAXITR
! 195: PARAMETER ( MAXITR = 6 )
! 196: DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO
! 197: PARAMETER ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0,
! 198: $ ONE = 1.0D0, PIOVER2 = 1.57079632679489662D0,
! 199: $ TEN = 10.0D0, ZERO = 0.0D0 )
! 200: COMPLEX*16 NEGONECOMPLEX
! 201: PARAMETER ( NEGONECOMPLEX = (-1.0D0,0.0D0) )
! 202: * ..
! 203: * .. Local Scalars ..
! 204: LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12,
! 205: $ RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
! 206: $ WANTV2T
! 207: INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
! 208: $ IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
! 209: $ LRWORKMIN, LRWORKOPT, MAXIT, MINI
! 210: DOUBLE PRECISION B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
! 211: $ EPS, MU, NU, R, SIGMA11, SIGMA21,
! 212: $ TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
! 213: $ UNFL, X1, X2, Y1, Y2
! 214: *
! 215: EXTERNAL DLARTGP, DLARTGS, DLAS2, XERBLA, ZLASR, ZSCAL,
! 216: $ ZSWAP
! 217: * ..
! 218: * .. External Functions ..
! 219: DOUBLE PRECISION DLAMCH
! 220: LOGICAL LSAME
! 221: EXTERNAL LSAME, DLAMCH
! 222: * ..
! 223: * .. Intrinsic Functions ..
! 224: INTRINSIC ABS, ATAN2, COS, MAX, MIN, SIN, SQRT
! 225: * ..
! 226: * .. Executable Statements ..
! 227: *
! 228: * Test input arguments
! 229: *
! 230: INFO = 0
! 231: LQUERY = LRWORK .EQ. -1
! 232: WANTU1 = LSAME( JOBU1, 'Y' )
! 233: WANTU2 = LSAME( JOBU2, 'Y' )
! 234: WANTV1T = LSAME( JOBV1T, 'Y' )
! 235: WANTV2T = LSAME( JOBV2T, 'Y' )
! 236: COLMAJOR = .NOT. LSAME( TRANS, 'T' )
! 237: *
! 238: IF( M .LT. 0 ) THEN
! 239: INFO = -6
! 240: ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
! 241: INFO = -7
! 242: ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
! 243: INFO = -8
! 244: ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN
! 245: INFO = -8
! 246: ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
! 247: INFO = -12
! 248: ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
! 249: INFO = -14
! 250: ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
! 251: INFO = -16
! 252: ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
! 253: INFO = -18
! 254: END IF
! 255: *
! 256: * Quick return if Q = 0
! 257: *
! 258: IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN
! 259: LRWORKMIN = 1
! 260: RWORK(1) = LRWORKMIN
! 261: RETURN
! 262: END IF
! 263: *
! 264: * Compute workspace
! 265: *
! 266: IF( INFO .EQ. 0 ) THEN
! 267: IU1CS = 1
! 268: IU1SN = IU1CS + Q
! 269: IU2CS = IU1SN + Q
! 270: IU2SN = IU2CS + Q
! 271: IV1TCS = IU2SN + Q
! 272: IV1TSN = IV1TCS + Q
! 273: IV2TCS = IV1TSN + Q
! 274: IV2TSN = IV2TCS + Q
! 275: LRWORKOPT = IV2TSN + Q - 1
! 276: LRWORKMIN = LRWORKOPT
! 277: RWORK(1) = LRWORKOPT
! 278: IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN
! 279: INFO = -28
! 280: END IF
! 281: END IF
! 282: *
! 283: IF( INFO .NE. 0 ) THEN
! 284: CALL XERBLA( 'ZBBCSD', -INFO )
! 285: RETURN
! 286: ELSE IF( LQUERY ) THEN
! 287: RETURN
! 288: END IF
! 289: *
! 290: * Get machine constants
! 291: *
! 292: EPS = DLAMCH( 'Epsilon' )
! 293: UNFL = DLAMCH( 'Safe minimum' )
! 294: TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) )
! 295: TOL = TOLMUL*EPS
! 296: THRESH = MAX( TOL, MAXITR*Q*Q*UNFL )
! 297: *
! 298: * Test for negligible sines or cosines
! 299: *
! 300: DO I = 1, Q
! 301: IF( THETA(I) .LT. THRESH ) THEN
! 302: THETA(I) = ZERO
! 303: ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
! 304: THETA(I) = PIOVER2
! 305: END IF
! 306: END DO
! 307: DO I = 1, Q-1
! 308: IF( PHI(I) .LT. THRESH ) THEN
! 309: PHI(I) = ZERO
! 310: ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
! 311: PHI(I) = PIOVER2
! 312: END IF
! 313: END DO
! 314: *
! 315: * Initial deflation
! 316: *
! 317: IMAX = Q
! 318: DO WHILE( ( IMAX .GT. 1 ) .AND. ( PHI(IMAX-1) .EQ. ZERO ) )
! 319: IMAX = IMAX - 1
! 320: END DO
! 321: IMIN = IMAX - 1
! 322: IF ( IMIN .GT. 1 ) THEN
! 323: DO WHILE( PHI(IMIN-1) .NE. ZERO )
! 324: IMIN = IMIN - 1
! 325: IF ( IMIN .LE. 1 ) EXIT
! 326: END DO
! 327: END IF
! 328: *
! 329: * Initialize iteration counter
! 330: *
! 331: MAXIT = MAXITR*Q*Q
! 332: ITER = 0
! 333: *
! 334: * Begin main iteration loop
! 335: *
! 336: DO WHILE( IMAX .GT. 1 )
! 337: *
! 338: * Compute the matrix entries
! 339: *
! 340: B11D(IMIN) = COS( THETA(IMIN) )
! 341: B21D(IMIN) = -SIN( THETA(IMIN) )
! 342: DO I = IMIN, IMAX - 1
! 343: B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) )
! 344: B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) )
! 345: B12D(I) = SIN( THETA(I) ) * COS( PHI(I) )
! 346: B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) )
! 347: B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) )
! 348: B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) )
! 349: B22D(I) = COS( THETA(I) ) * COS( PHI(I) )
! 350: B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) )
! 351: END DO
! 352: B12D(IMAX) = SIN( THETA(IMAX) )
! 353: B22D(IMAX) = COS( THETA(IMAX) )
! 354: *
! 355: * Abort if not converging; otherwise, increment ITER
! 356: *
! 357: IF( ITER .GT. MAXIT ) THEN
! 358: INFO = 0
! 359: DO I = 1, Q
! 360: IF( PHI(I) .NE. ZERO )
! 361: $ INFO = INFO + 1
! 362: END DO
! 363: RETURN
! 364: END IF
! 365: *
! 366: ITER = ITER + IMAX - IMIN
! 367: *
! 368: * Compute shifts
! 369: *
! 370: THETAMAX = THETA(IMIN)
! 371: THETAMIN = THETA(IMIN)
! 372: DO I = IMIN+1, IMAX
! 373: IF( THETA(I) > THETAMAX )
! 374: $ THETAMAX = THETA(I)
! 375: IF( THETA(I) < THETAMIN )
! 376: $ THETAMIN = THETA(I)
! 377: END DO
! 378: *
! 379: IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN
! 380: *
! 381: * Zero on diagonals of B11 and B22; induce deflation with a
! 382: * zero shift
! 383: *
! 384: MU = ZERO
! 385: NU = ONE
! 386: *
! 387: ELSE IF( THETAMIN .LT. THRESH ) THEN
! 388: *
! 389: * Zero on diagonals of B12 and B22; induce deflation with a
! 390: * zero shift
! 391: *
! 392: MU = ONE
! 393: NU = ZERO
! 394: *
! 395: ELSE
! 396: *
! 397: * Compute shifts for B11 and B21 and use the lesser
! 398: *
! 399: CALL DLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11,
! 400: $ DUMMY )
! 401: CALL DLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21,
! 402: $ DUMMY )
! 403: *
! 404: IF( SIGMA11 .LE. SIGMA21 ) THEN
! 405: MU = SIGMA11
! 406: NU = SQRT( ONE - MU**2 )
! 407: IF( MU .LT. THRESH ) THEN
! 408: MU = ZERO
! 409: NU = ONE
! 410: END IF
! 411: ELSE
! 412: NU = SIGMA21
! 413: MU = SQRT( 1.0 - NU**2 )
! 414: IF( NU .LT. THRESH ) THEN
! 415: MU = ONE
! 416: NU = ZERO
! 417: END IF
! 418: END IF
! 419: END IF
! 420: *
! 421: * Rotate to produce bulges in B11 and B21
! 422: *
! 423: IF( MU .LE. NU ) THEN
! 424: CALL DLARTGS( B11D(IMIN), B11E(IMIN), MU,
! 425: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
! 426: ELSE
! 427: CALL DLARTGS( B21D(IMIN), B21E(IMIN), NU,
! 428: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) )
! 429: END IF
! 430: *
! 431: TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) +
! 432: $ RWORK(IV1TSN+IMIN-1)*B11E(IMIN)
! 433: B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) -
! 434: $ RWORK(IV1TSN+IMIN-1)*B11D(IMIN)
! 435: B11D(IMIN) = TEMP
! 436: B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1)
! 437: B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1)
! 438: TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) +
! 439: $ RWORK(IV1TSN+IMIN-1)*B21E(IMIN)
! 440: B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) -
! 441: $ RWORK(IV1TSN+IMIN-1)*B21D(IMIN)
! 442: B21D(IMIN) = TEMP
! 443: B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1)
! 444: B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1)
! 445: *
! 446: * Compute THETA(IMIN)
! 447: *
! 448: THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ),
! 449: $ SQRT( B11D(IMIN)**2+B11BULGE**2 ) )
! 450: *
! 451: * Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
! 452: *
! 453: IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
! 454: CALL DLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1),
! 455: $ RWORK(IU1CS+IMIN-1), R )
! 456: ELSE IF( MU .LE. NU ) THEN
! 457: CALL DLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU,
! 458: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
! 459: ELSE
! 460: CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
! 461: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
! 462: END IF
! 463: IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
! 464: CALL DLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1),
! 465: $ RWORK(IU2CS+IMIN-1), R )
! 466: ELSE IF( NU .LT. MU ) THEN
! 467: CALL DLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU,
! 468: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
! 469: ELSE
! 470: CALL DLARTGS( B22D(IMIN), B22E(IMIN), MU,
! 471: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) )
! 472: END IF
! 473: RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1)
! 474: RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1)
! 475: *
! 476: TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) +
! 477: $ RWORK(IU1SN+IMIN-1)*B11D(IMIN+1)
! 478: B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) -
! 479: $ RWORK(IU1SN+IMIN-1)*B11E(IMIN)
! 480: B11E(IMIN) = TEMP
! 481: IF( IMAX .GT. IMIN+1 ) THEN
! 482: B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1)
! 483: B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1)
! 484: END IF
! 485: TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) +
! 486: $ RWORK(IU1SN+IMIN-1)*B12E(IMIN)
! 487: B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) -
! 488: $ RWORK(IU1SN+IMIN-1)*B12D(IMIN)
! 489: B12D(IMIN) = TEMP
! 490: B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1)
! 491: B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1)
! 492: TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) +
! 493: $ RWORK(IU2SN+IMIN-1)*B21D(IMIN+1)
! 494: B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) -
! 495: $ RWORK(IU2SN+IMIN-1)*B21E(IMIN)
! 496: B21E(IMIN) = TEMP
! 497: IF( IMAX .GT. IMIN+1 ) THEN
! 498: B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1)
! 499: B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1)
! 500: END IF
! 501: TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) +
! 502: $ RWORK(IU2SN+IMIN-1)*B22E(IMIN)
! 503: B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) -
! 504: $ RWORK(IU2SN+IMIN-1)*B22D(IMIN)
! 505: B22D(IMIN) = TEMP
! 506: B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1)
! 507: B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1)
! 508: *
! 509: * Inner loop: chase bulges from B11(IMIN,IMIN+2),
! 510: * B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
! 511: * bottom-right
! 512: *
! 513: DO I = IMIN+1, IMAX-1
! 514: *
! 515: * Compute PHI(I-1)
! 516: *
! 517: X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1)
! 518: X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE
! 519: Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1)
! 520: Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE
! 521: *
! 522: PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) )
! 523: *
! 524: * Determine if there are bulges to chase or if a new direct
! 525: * summand has been reached
! 526: *
! 527: RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
! 528: RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
! 529: RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
! 530: RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
! 531: *
! 532: * If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
! 533: * B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
! 534: * chasing by applying the original shift again.
! 535: *
! 536: IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN
! 537: CALL DLARTGP( X2, X1, RWORK(IV1TSN+I-1),
! 538: $ RWORK(IV1TCS+I-1), R )
! 539: ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN
! 540: CALL DLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1),
! 541: $ RWORK(IV1TCS+I-1), R )
! 542: ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN
! 543: CALL DLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1),
! 544: $ RWORK(IV1TCS+I-1), R )
! 545: ELSE IF( MU .LE. NU ) THEN
! 546: CALL DLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1),
! 547: $ RWORK(IV1TSN+I-1) )
! 548: ELSE
! 549: CALL DLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1),
! 550: $ RWORK(IV1TSN+I-1) )
! 551: END IF
! 552: RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1)
! 553: RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1)
! 554: IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
! 555: CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1),
! 556: $ RWORK(IV2TCS+I-1-1), R )
! 557: ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
! 558: CALL DLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1),
! 559: $ RWORK(IV2TCS+I-1-1), R )
! 560: ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
! 561: CALL DLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1),
! 562: $ RWORK(IV2TCS+I-1-1), R )
! 563: ELSE IF( NU .LT. MU ) THEN
! 564: CALL DLARTGS( B12E(I-1), B12D(I), NU,
! 565: $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
! 566: ELSE
! 567: CALL DLARTGS( B22E(I-1), B22D(I), MU,
! 568: $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) )
! 569: END IF
! 570: *
! 571: TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I)
! 572: B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) -
! 573: $ RWORK(IV1TSN+I-1)*B11D(I)
! 574: B11D(I) = TEMP
! 575: B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1)
! 576: B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1)
! 577: TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I)
! 578: B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) -
! 579: $ RWORK(IV1TSN+I-1)*B21D(I)
! 580: B21D(I) = TEMP
! 581: B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1)
! 582: B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1)
! 583: TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) +
! 584: $ RWORK(IV2TSN+I-1-1)*B12D(I)
! 585: B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) -
! 586: $ RWORK(IV2TSN+I-1-1)*B12E(I-1)
! 587: B12E(I-1) = TEMP
! 588: B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I)
! 589: B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I)
! 590: TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) +
! 591: $ RWORK(IV2TSN+I-1-1)*B22D(I)
! 592: B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) -
! 593: $ RWORK(IV2TSN+I-1-1)*B22E(I-1)
! 594: B22E(I-1) = TEMP
! 595: B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I)
! 596: B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I)
! 597: *
! 598: * Compute THETA(I)
! 599: *
! 600: X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1)
! 601: X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE
! 602: Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1)
! 603: Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE
! 604: *
! 605: THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) )
! 606: *
! 607: * Determine if there are bulges to chase or if a new direct
! 608: * summand has been reached
! 609: *
! 610: RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
! 611: RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
! 612: RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
! 613: RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
! 614: *
! 615: * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
! 616: * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
! 617: * chasing by applying the original shift again.
! 618: *
! 619: IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN
! 620: CALL DLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1),
! 621: $ R )
! 622: ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN
! 623: CALL DLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1),
! 624: $ RWORK(IU1CS+I-1), R )
! 625: ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN
! 626: CALL DLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1),
! 627: $ RWORK(IU1CS+I-1), R )
! 628: ELSE IF( MU .LE. NU ) THEN
! 629: CALL DLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1),
! 630: $ RWORK(IU1SN+I-1) )
! 631: ELSE
! 632: CALL DLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1),
! 633: $ RWORK(IU1SN+I-1) )
! 634: END IF
! 635: IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN
! 636: CALL DLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1),
! 637: $ R )
! 638: ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN
! 639: CALL DLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1),
! 640: $ RWORK(IU2CS+I-1), R )
! 641: ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN
! 642: CALL DLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1),
! 643: $ RWORK(IU2CS+I-1), R )
! 644: ELSE IF( NU .LT. MU ) THEN
! 645: CALL DLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1),
! 646: $ RWORK(IU2SN+I-1) )
! 647: ELSE
! 648: CALL DLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1),
! 649: $ RWORK(IU2SN+I-1) )
! 650: END IF
! 651: RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1)
! 652: RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1)
! 653: *
! 654: TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1)
! 655: B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) -
! 656: $ RWORK(IU1SN+I-1)*B11E(I)
! 657: B11E(I) = TEMP
! 658: IF( I .LT. IMAX - 1 ) THEN
! 659: B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1)
! 660: B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1)
! 661: END IF
! 662: TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1)
! 663: B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) -
! 664: $ RWORK(IU2SN+I-1)*B21E(I)
! 665: B21E(I) = TEMP
! 666: IF( I .LT. IMAX - 1 ) THEN
! 667: B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1)
! 668: B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1)
! 669: END IF
! 670: TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I)
! 671: B12E(I) = RWORK(IU1CS+I-1)*B12E(I) -
! 672: $ RWORK(IU1SN+I-1)*B12D(I)
! 673: B12D(I) = TEMP
! 674: B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1)
! 675: B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1)
! 676: TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I)
! 677: B22E(I) = RWORK(IU2CS+I-1)*B22E(I) -
! 678: $ RWORK(IU2SN+I-1)*B22D(I)
! 679: B22D(I) = TEMP
! 680: B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1)
! 681: B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1)
! 682: *
! 683: END DO
! 684: *
! 685: * Compute PHI(IMAX-1)
! 686: *
! 687: X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) +
! 688: $ COS(THETA(IMAX-1))*B21E(IMAX-1)
! 689: Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) +
! 690: $ COS(THETA(IMAX-1))*B22D(IMAX-1)
! 691: Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE
! 692: *
! 693: PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) )
! 694: *
! 695: * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
! 696: *
! 697: RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
! 698: RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
! 699: *
! 700: IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
! 701: CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1),
! 702: $ RWORK(IV2TCS+IMAX-1-1), R )
! 703: ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN
! 704: CALL DLARTGP( B12BULGE, B12D(IMAX-1),
! 705: $ RWORK(IV2TSN+IMAX-1-1),
! 706: $ RWORK(IV2TCS+IMAX-1-1), R )
! 707: ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN
! 708: CALL DLARTGP( B22BULGE, B22D(IMAX-1),
! 709: $ RWORK(IV2TSN+IMAX-1-1),
! 710: $ RWORK(IV2TCS+IMAX-1-1), R )
! 711: ELSE IF( NU .LT. MU ) THEN
! 712: CALL DLARTGS( B12E(IMAX-1), B12D(IMAX), NU,
! 713: $ RWORK(IV2TCS+IMAX-1-1),
! 714: $ RWORK(IV2TSN+IMAX-1-1) )
! 715: ELSE
! 716: CALL DLARTGS( B22E(IMAX-1), B22D(IMAX), MU,
! 717: $ RWORK(IV2TCS+IMAX-1-1),
! 718: $ RWORK(IV2TSN+IMAX-1-1) )
! 719: END IF
! 720: *
! 721: TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) +
! 722: $ RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX)
! 723: B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) -
! 724: $ RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1)
! 725: B12E(IMAX-1) = TEMP
! 726: TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) +
! 727: $ RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX)
! 728: B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) -
! 729: $ RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1)
! 730: B22E(IMAX-1) = TEMP
! 731: *
! 732: * Update singular vectors
! 733: *
! 734: IF( WANTU1 ) THEN
! 735: IF( COLMAJOR ) THEN
! 736: CALL ZLASR( 'R', 'V', 'F', P, IMAX-IMIN+1,
! 737: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
! 738: $ U1(1,IMIN), LDU1 )
! 739: ELSE
! 740: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, P,
! 741: $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1),
! 742: $ U1(IMIN,1), LDU1 )
! 743: END IF
! 744: END IF
! 745: IF( WANTU2 ) THEN
! 746: IF( COLMAJOR ) THEN
! 747: CALL ZLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1,
! 748: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
! 749: $ U2(1,IMIN), LDU2 )
! 750: ELSE
! 751: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P,
! 752: $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1),
! 753: $ U2(IMIN,1), LDU2 )
! 754: END IF
! 755: END IF
! 756: IF( WANTV1T ) THEN
! 757: IF( COLMAJOR ) THEN
! 758: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q,
! 759: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
! 760: $ V1T(IMIN,1), LDV1T )
! 761: ELSE
! 762: CALL ZLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1,
! 763: $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1),
! 764: $ V1T(1,IMIN), LDV1T )
! 765: END IF
! 766: END IF
! 767: IF( WANTV2T ) THEN
! 768: IF( COLMAJOR ) THEN
! 769: CALL ZLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q,
! 770: $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
! 771: $ V2T(IMIN,1), LDV2T )
! 772: ELSE
! 773: CALL ZLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1,
! 774: $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1),
! 775: $ V2T(1,IMIN), LDV2T )
! 776: END IF
! 777: END IF
! 778: *
! 779: * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
! 780: *
! 781: IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN
! 782: B11D(IMAX) = -B11D(IMAX)
! 783: B21D(IMAX) = -B21D(IMAX)
! 784: IF( WANTV1T ) THEN
! 785: IF( COLMAJOR ) THEN
! 786: CALL ZSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T )
! 787: ELSE
! 788: CALL ZSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 )
! 789: END IF
! 790: END IF
! 791: END IF
! 792: *
! 793: * Compute THETA(IMAX)
! 794: *
! 795: X1 = COS(PHI(IMAX-1))*B11D(IMAX) +
! 796: $ SIN(PHI(IMAX-1))*B12E(IMAX-1)
! 797: Y1 = COS(PHI(IMAX-1))*B21D(IMAX) +
! 798: $ SIN(PHI(IMAX-1))*B22E(IMAX-1)
! 799: *
! 800: THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) )
! 801: *
! 802: * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
! 803: * and B22(IMAX,IMAX-1)
! 804: *
! 805: IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN
! 806: B12D(IMAX) = -B12D(IMAX)
! 807: IF( WANTU1 ) THEN
! 808: IF( COLMAJOR ) THEN
! 809: CALL ZSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 )
! 810: ELSE
! 811: CALL ZSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 )
! 812: END IF
! 813: END IF
! 814: END IF
! 815: IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN
! 816: B22D(IMAX) = -B22D(IMAX)
! 817: IF( WANTU2 ) THEN
! 818: IF( COLMAJOR ) THEN
! 819: CALL ZSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 )
! 820: ELSE
! 821: CALL ZSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 )
! 822: END IF
! 823: END IF
! 824: END IF
! 825: *
! 826: * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
! 827: *
! 828: IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN
! 829: IF( WANTV2T ) THEN
! 830: IF( COLMAJOR ) THEN
! 831: CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T )
! 832: ELSE
! 833: CALL ZSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 )
! 834: END IF
! 835: END IF
! 836: END IF
! 837: *
! 838: * Test for negligible sines or cosines
! 839: *
! 840: DO I = IMIN, IMAX
! 841: IF( THETA(I) .LT. THRESH ) THEN
! 842: THETA(I) = ZERO
! 843: ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN
! 844: THETA(I) = PIOVER2
! 845: END IF
! 846: END DO
! 847: DO I = IMIN, IMAX-1
! 848: IF( PHI(I) .LT. THRESH ) THEN
! 849: PHI(I) = ZERO
! 850: ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN
! 851: PHI(I) = PIOVER2
! 852: END IF
! 853: END DO
! 854: *
! 855: * Deflate
! 856: *
! 857: IF (IMAX .GT. 1) THEN
! 858: DO WHILE( PHI(IMAX-1) .EQ. ZERO )
! 859: IMAX = IMAX - 1
! 860: IF (IMAX .LE. 1) EXIT
! 861: END DO
! 862: END IF
! 863: IF( IMIN .GT. IMAX - 1 )
! 864: $ IMIN = IMAX - 1
! 865: IF (IMIN .GT. 1) THEN
! 866: DO WHILE (PHI(IMIN-1) .NE. ZERO)
! 867: IMIN = IMIN - 1
! 868: IF (IMIN .LE. 1) EXIT
! 869: END DO
! 870: END IF
! 871: *
! 872: * Repeat main iteration loop
! 873: *
! 874: END DO
! 875: *
! 876: * Postprocessing: order THETA from least to greatest
! 877: *
! 878: DO I = 1, Q
! 879: *
! 880: MINI = I
! 881: THETAMIN = THETA(I)
! 882: DO J = I+1, Q
! 883: IF( THETA(J) .LT. THETAMIN ) THEN
! 884: MINI = J
! 885: THETAMIN = THETA(J)
! 886: END IF
! 887: END DO
! 888: *
! 889: IF( MINI .NE. I ) THEN
! 890: THETA(MINI) = THETA(I)
! 891: THETA(I) = THETAMIN
! 892: IF( COLMAJOR ) THEN
! 893: IF( WANTU1 )
! 894: $ CALL ZSWAP( P, U1(1,I), 1, U1(1,MINI), 1 )
! 895: IF( WANTU2 )
! 896: $ CALL ZSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 )
! 897: IF( WANTV1T )
! 898: $ CALL ZSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T )
! 899: IF( WANTV2T )
! 900: $ CALL ZSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1),
! 901: $ LDV2T )
! 902: ELSE
! 903: IF( WANTU1 )
! 904: $ CALL ZSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 )
! 905: IF( WANTU2 )
! 906: $ CALL ZSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 )
! 907: IF( WANTV1T )
! 908: $ CALL ZSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 )
! 909: IF( WANTV2T )
! 910: $ CALL ZSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 )
! 911: END IF
! 912: END IF
! 913: *
! 914: END DO
! 915: *
! 916: RETURN
! 917: *
! 918: * End of ZBBCSD
! 919: *
! 920: END
! 921:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zbbcsd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack