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, a P-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, LDU1 >= MAX(1,P). 161*> \endverbatim 162*> 163*> \param[in,out] U2 164*> \verbatim 165*> U2 is COMPLEX*16 array, dimension (LDU2,M-P) 166*> On entry, an (M-P)-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, LDU2 >= MAX(1,M-P). 175*> \endverbatim 176*> 177*> \param[in,out] V1T 178*> \verbatim 179*> V1T is COMPLEX*16 array, dimension (LDV1T,Q) 180*> On entry, a Q-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, LDV1T >= MAX(1,Q). 189*> \endverbatim 190*> 191*> \param[in,out] V2T 192*> \verbatim 193*> V2T is COMPLEX*16 array, dimension (LDV2T,M-Q) 194*> On entry, an (M-Q)-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, LDV2T >= MAX(1,M-Q). 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 ZBBCSD converges, B21D contains the negative sines of 244*> THETA(1), ..., THETA(Q). If ZBBCSD 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 ZBBCSD converges, B21E contains zeros. If ZBBCSD 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 ZBBCSD converges, B22D contains the negative sines of 261*> THETA(1), ..., THETA(Q). If ZBBCSD 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 ZBBCSD converges, B22E contains zeros. If ZBBCSD 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,LRWORK)) 277*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 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 June 2016 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.7.1) -- 336* -- LAPACK is a software package provided by Univ. of Tennessee, -- 337* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 338* June 2016 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 1086