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