1*> \brief \b CBBCSD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CBBCSD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cbbcsd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cbbcsd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cbbcsd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CBBCSD( 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* REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ), 32* $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 33* $ PHI( * ), THETA( * ), RWORK( * ) 34* COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 35* $ V2T( LDV2T, * ) 36* .. 37* 38* 39*> \par Purpose: 40* ============= 41*> 42*> \verbatim 43*> 44*> CBBCSD 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 CUNCSD 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 REAL 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 REAL 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 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 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 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 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 REAL array, dimension (Q) 209*> When CBBCSD converges, B11D contains the cosines of THETA(1), 210*> ..., THETA(Q). If CBBCSD 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 REAL array, dimension (Q-1) 218*> When CBBCSD converges, B11E contains zeros. If CBBCSD 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 REAL array, dimension (Q) 226*> When CBBCSD converges, B12D contains the negative sines of 227*> THETA(1), ..., THETA(Q). If CBBCSD 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 REAL array, dimension (Q-1) 235*> When CBBCSD converges, B12E contains zeros. If CBBCSD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 CBBCSD 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 REAL, 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 complexOTHERcomputational 328* 329* ===================================================================== 330 SUBROUTINE CBBCSD( 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 REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ), 346 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ), 347 $ PHI( * ), THETA( * ), RWORK( * ) 348 COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 349 $ V2T( LDV2T, * ) 350* .. 351* 352* =================================================================== 353* 354* .. Parameters .. 355 INTEGER MAXITR 356 PARAMETER ( MAXITR = 6 ) 357 REAL HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO 358 PARAMETER ( HUNDRED = 100.0E0, MEIGHTH = -0.125E0, 359 $ ONE = 1.0E0, PIOVER2 = 1.57079632679489662E0, 360 $ TEN = 10.0E0, ZERO = 0.0E0 ) 361 COMPLEX NEGONECOMPLEX 362 PARAMETER ( NEGONECOMPLEX = (-1.0E0,0.0E0) ) 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 REAL 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 Subroutines .. 377 EXTERNAL CLASR, CSCAL, CSWAP, SLARTGP, SLARTGS, SLAS2, 378 $ XERBLA 379* .. 380* .. External Functions .. 381 REAL SLAMCH 382 LOGICAL LSAME 383 EXTERNAL LSAME, SLAMCH 384* .. 385* .. Intrinsic Functions .. 386 INTRINSIC ABS, ATAN2, COS, MAX, MIN, SIN, SQRT 387* .. 388* .. Executable Statements .. 389* 390* Test input arguments 391* 392 INFO = 0 393 LQUERY = LRWORK .EQ. -1 394 WANTU1 = LSAME( JOBU1, 'Y' ) 395 WANTU2 = LSAME( JOBU2, 'Y' ) 396 WANTV1T = LSAME( JOBV1T, 'Y' ) 397 WANTV2T = LSAME( JOBV2T, 'Y' ) 398 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 399* 400 IF( M .LT. 0 ) THEN 401 INFO = -6 402 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 403 INFO = -7 404 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 405 INFO = -8 406 ELSE IF( Q .GT. P .OR. Q .GT. M-P .OR. Q .GT. M-Q ) THEN 407 INFO = -8 408 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 409 INFO = -12 410 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 411 INFO = -14 412 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 413 INFO = -16 414 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 415 INFO = -18 416 END IF 417* 418* Quick return if Q = 0 419* 420 IF( INFO .EQ. 0 .AND. Q .EQ. 0 ) THEN 421 LRWORKMIN = 1 422 RWORK(1) = LRWORKMIN 423 RETURN 424 END IF 425* 426* Compute workspace 427* 428 IF( INFO .EQ. 0 ) THEN 429 IU1CS = 1 430 IU1SN = IU1CS + Q 431 IU2CS = IU1SN + Q 432 IU2SN = IU2CS + Q 433 IV1TCS = IU2SN + Q 434 IV1TSN = IV1TCS + Q 435 IV2TCS = IV1TSN + Q 436 IV2TSN = IV2TCS + Q 437 LRWORKOPT = IV2TSN + Q - 1 438 LRWORKMIN = LRWORKOPT 439 RWORK(1) = LRWORKOPT 440 IF( LRWORK .LT. LRWORKMIN .AND. .NOT. LQUERY ) THEN 441 INFO = -28 442 END IF 443 END IF 444* 445 IF( INFO .NE. 0 ) THEN 446 CALL XERBLA( 'CBBCSD', -INFO ) 447 RETURN 448 ELSE IF( LQUERY ) THEN 449 RETURN 450 END IF 451* 452* Get machine constants 453* 454 EPS = SLAMCH( 'Epsilon' ) 455 UNFL = SLAMCH( 'Safe minimum' ) 456 TOLMUL = MAX( TEN, MIN( HUNDRED, EPS**MEIGHTH ) ) 457 TOL = TOLMUL*EPS 458 THRESH = MAX( TOL, MAXITR*Q*Q*UNFL ) 459* 460* Test for negligible sines or cosines 461* 462 DO I = 1, Q 463 IF( THETA(I) .LT. THRESH ) THEN 464 THETA(I) = ZERO 465 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 466 THETA(I) = PIOVER2 467 END IF 468 END DO 469 DO I = 1, Q-1 470 IF( PHI(I) .LT. THRESH ) THEN 471 PHI(I) = ZERO 472 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 473 PHI(I) = PIOVER2 474 END IF 475 END DO 476* 477* Initial deflation 478* 479 IMAX = Q 480 DO WHILE( IMAX .GT. 1 ) 481 IF( PHI(IMAX-1) .NE. ZERO ) THEN 482 EXIT 483 END IF 484 IMAX = IMAX - 1 485 END DO 486 IMIN = IMAX - 1 487 IF ( IMIN .GT. 1 ) THEN 488 DO WHILE( PHI(IMIN-1) .NE. ZERO ) 489 IMIN = IMIN - 1 490 IF ( IMIN .LE. 1 ) EXIT 491 END DO 492 END IF 493* 494* Initialize iteration counter 495* 496 MAXIT = MAXITR*Q*Q 497 ITER = 0 498* 499* Begin main iteration loop 500* 501 DO WHILE( IMAX .GT. 1 ) 502* 503* Compute the matrix entries 504* 505 B11D(IMIN) = COS( THETA(IMIN) ) 506 B21D(IMIN) = -SIN( THETA(IMIN) ) 507 DO I = IMIN, IMAX - 1 508 B11E(I) = -SIN( THETA(I) ) * SIN( PHI(I) ) 509 B11D(I+1) = COS( THETA(I+1) ) * COS( PHI(I) ) 510 B12D(I) = SIN( THETA(I) ) * COS( PHI(I) ) 511 B12E(I) = COS( THETA(I+1) ) * SIN( PHI(I) ) 512 B21E(I) = -COS( THETA(I) ) * SIN( PHI(I) ) 513 B21D(I+1) = -SIN( THETA(I+1) ) * COS( PHI(I) ) 514 B22D(I) = COS( THETA(I) ) * COS( PHI(I) ) 515 B22E(I) = -SIN( THETA(I+1) ) * SIN( PHI(I) ) 516 END DO 517 B12D(IMAX) = SIN( THETA(IMAX) ) 518 B22D(IMAX) = COS( THETA(IMAX) ) 519* 520* Abort if not converging; otherwise, increment ITER 521* 522 IF( ITER .GT. MAXIT ) THEN 523 INFO = 0 524 DO I = 1, Q 525 IF( PHI(I) .NE. ZERO ) 526 $ INFO = INFO + 1 527 END DO 528 RETURN 529 END IF 530* 531 ITER = ITER + IMAX - IMIN 532* 533* Compute shifts 534* 535 THETAMAX = THETA(IMIN) 536 THETAMIN = THETA(IMIN) 537 DO I = IMIN+1, IMAX 538 IF( THETA(I) > THETAMAX ) 539 $ THETAMAX = THETA(I) 540 IF( THETA(I) < THETAMIN ) 541 $ THETAMIN = THETA(I) 542 END DO 543* 544 IF( THETAMAX .GT. PIOVER2 - THRESH ) THEN 545* 546* Zero on diagonals of B11 and B22; induce deflation with a 547* zero shift 548* 549 MU = ZERO 550 NU = ONE 551* 552 ELSE IF( THETAMIN .LT. THRESH ) THEN 553* 554* Zero on diagonals of B12 and B22; induce deflation with a 555* zero shift 556* 557 MU = ONE 558 NU = ZERO 559* 560 ELSE 561* 562* Compute shifts for B11 and B21 and use the lesser 563* 564 CALL SLAS2( B11D(IMAX-1), B11E(IMAX-1), B11D(IMAX), SIGMA11, 565 $ DUMMY ) 566 CALL SLAS2( B21D(IMAX-1), B21E(IMAX-1), B21D(IMAX), SIGMA21, 567 $ DUMMY ) 568* 569 IF( SIGMA11 .LE. SIGMA21 ) THEN 570 MU = SIGMA11 571 NU = SQRT( ONE - MU**2 ) 572 IF( MU .LT. THRESH ) THEN 573 MU = ZERO 574 NU = ONE 575 END IF 576 ELSE 577 NU = SIGMA21 578 MU = SQRT( 1.0 - NU**2 ) 579 IF( NU .LT. THRESH ) THEN 580 MU = ONE 581 NU = ZERO 582 END IF 583 END IF 584 END IF 585* 586* Rotate to produce bulges in B11 and B21 587* 588 IF( MU .LE. NU ) THEN 589 CALL SLARTGS( B11D(IMIN), B11E(IMIN), MU, 590 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 591 ELSE 592 CALL SLARTGS( B21D(IMIN), B21E(IMIN), NU, 593 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1) ) 594 END IF 595* 596 TEMP = RWORK(IV1TCS+IMIN-1)*B11D(IMIN) + 597 $ RWORK(IV1TSN+IMIN-1)*B11E(IMIN) 598 B11E(IMIN) = RWORK(IV1TCS+IMIN-1)*B11E(IMIN) - 599 $ RWORK(IV1TSN+IMIN-1)*B11D(IMIN) 600 B11D(IMIN) = TEMP 601 B11BULGE = RWORK(IV1TSN+IMIN-1)*B11D(IMIN+1) 602 B11D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B11D(IMIN+1) 603 TEMP = RWORK(IV1TCS+IMIN-1)*B21D(IMIN) + 604 $ RWORK(IV1TSN+IMIN-1)*B21E(IMIN) 605 B21E(IMIN) = RWORK(IV1TCS+IMIN-1)*B21E(IMIN) - 606 $ RWORK(IV1TSN+IMIN-1)*B21D(IMIN) 607 B21D(IMIN) = TEMP 608 B21BULGE = RWORK(IV1TSN+IMIN-1)*B21D(IMIN+1) 609 B21D(IMIN+1) = RWORK(IV1TCS+IMIN-1)*B21D(IMIN+1) 610* 611* Compute THETA(IMIN) 612* 613 THETA( IMIN ) = ATAN2( SQRT( B21D(IMIN)**2+B21BULGE**2 ), 614 $ SQRT( B11D(IMIN)**2+B11BULGE**2 ) ) 615* 616* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN) 617* 618 IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN 619 CALL SLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1), 620 $ RWORK(IU1CS+IMIN-1), R ) 621 ELSE IF( MU .LE. NU ) THEN 622 CALL SLARTGS( B11E( IMIN ), B11D( IMIN + 1 ), MU, 623 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 624 ELSE 625 CALL SLARTGS( B12D( IMIN ), B12E( IMIN ), NU, 626 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) ) 627 END IF 628 IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN 629 CALL SLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1), 630 $ RWORK(IU2CS+IMIN-1), R ) 631 ELSE IF( NU .LT. MU ) THEN 632 CALL SLARTGS( B21E( IMIN ), B21D( IMIN + 1 ), NU, 633 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 634 ELSE 635 CALL SLARTGS( B22D(IMIN), B22E(IMIN), MU, 636 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1) ) 637 END IF 638 RWORK(IU2CS+IMIN-1) = -RWORK(IU2CS+IMIN-1) 639 RWORK(IU2SN+IMIN-1) = -RWORK(IU2SN+IMIN-1) 640* 641 TEMP = RWORK(IU1CS+IMIN-1)*B11E(IMIN) + 642 $ RWORK(IU1SN+IMIN-1)*B11D(IMIN+1) 643 B11D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11D(IMIN+1) - 644 $ RWORK(IU1SN+IMIN-1)*B11E(IMIN) 645 B11E(IMIN) = TEMP 646 IF( IMAX .GT. IMIN+1 ) THEN 647 B11BULGE = RWORK(IU1SN+IMIN-1)*B11E(IMIN+1) 648 B11E(IMIN+1) = RWORK(IU1CS+IMIN-1)*B11E(IMIN+1) 649 END IF 650 TEMP = RWORK(IU1CS+IMIN-1)*B12D(IMIN) + 651 $ RWORK(IU1SN+IMIN-1)*B12E(IMIN) 652 B12E(IMIN) = RWORK(IU1CS+IMIN-1)*B12E(IMIN) - 653 $ RWORK(IU1SN+IMIN-1)*B12D(IMIN) 654 B12D(IMIN) = TEMP 655 B12BULGE = RWORK(IU1SN+IMIN-1)*B12D(IMIN+1) 656 B12D(IMIN+1) = RWORK(IU1CS+IMIN-1)*B12D(IMIN+1) 657 TEMP = RWORK(IU2CS+IMIN-1)*B21E(IMIN) + 658 $ RWORK(IU2SN+IMIN-1)*B21D(IMIN+1) 659 B21D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21D(IMIN+1) - 660 $ RWORK(IU2SN+IMIN-1)*B21E(IMIN) 661 B21E(IMIN) = TEMP 662 IF( IMAX .GT. IMIN+1 ) THEN 663 B21BULGE = RWORK(IU2SN+IMIN-1)*B21E(IMIN+1) 664 B21E(IMIN+1) = RWORK(IU2CS+IMIN-1)*B21E(IMIN+1) 665 END IF 666 TEMP = RWORK(IU2CS+IMIN-1)*B22D(IMIN) + 667 $ RWORK(IU2SN+IMIN-1)*B22E(IMIN) 668 B22E(IMIN) = RWORK(IU2CS+IMIN-1)*B22E(IMIN) - 669 $ RWORK(IU2SN+IMIN-1)*B22D(IMIN) 670 B22D(IMIN) = TEMP 671 B22BULGE = RWORK(IU2SN+IMIN-1)*B22D(IMIN+1) 672 B22D(IMIN+1) = RWORK(IU2CS+IMIN-1)*B22D(IMIN+1) 673* 674* Inner loop: chase bulges from B11(IMIN,IMIN+2), 675* B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to 676* bottom-right 677* 678 DO I = IMIN+1, IMAX-1 679* 680* Compute PHI(I-1) 681* 682 X1 = SIN(THETA(I-1))*B11E(I-1) + COS(THETA(I-1))*B21E(I-1) 683 X2 = SIN(THETA(I-1))*B11BULGE + COS(THETA(I-1))*B21BULGE 684 Y1 = SIN(THETA(I-1))*B12D(I-1) + COS(THETA(I-1))*B22D(I-1) 685 Y2 = SIN(THETA(I-1))*B12BULGE + COS(THETA(I-1))*B22BULGE 686* 687 PHI(I-1) = ATAN2( SQRT(X1**2+X2**2), SQRT(Y1**2+Y2**2) ) 688* 689* Determine if there are bulges to chase or if a new direct 690* summand has been reached 691* 692 RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2 693 RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2 694 RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2 695 RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2 696* 697* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I), 698* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge- 699* chasing by applying the original shift again. 700* 701 IF( .NOT. RESTART11 .AND. .NOT. RESTART21 ) THEN 702 CALL SLARTGP( X2, X1, RWORK(IV1TSN+I-1), 703 $ RWORK(IV1TCS+I-1), R ) 704 ELSE IF( .NOT. RESTART11 .AND. RESTART21 ) THEN 705 CALL SLARTGP( B11BULGE, B11E(I-1), RWORK(IV1TSN+I-1), 706 $ RWORK(IV1TCS+I-1), R ) 707 ELSE IF( RESTART11 .AND. .NOT. RESTART21 ) THEN 708 CALL SLARTGP( B21BULGE, B21E(I-1), RWORK(IV1TSN+I-1), 709 $ RWORK(IV1TCS+I-1), R ) 710 ELSE IF( MU .LE. NU ) THEN 711 CALL SLARTGS( B11D(I), B11E(I), MU, RWORK(IV1TCS+I-1), 712 $ RWORK(IV1TSN+I-1) ) 713 ELSE 714 CALL SLARTGS( B21D(I), B21E(I), NU, RWORK(IV1TCS+I-1), 715 $ RWORK(IV1TSN+I-1) ) 716 END IF 717 RWORK(IV1TCS+I-1) = -RWORK(IV1TCS+I-1) 718 RWORK(IV1TSN+I-1) = -RWORK(IV1TSN+I-1) 719 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 720 CALL SLARTGP( Y2, Y1, RWORK(IV2TSN+I-1-1), 721 $ RWORK(IV2TCS+I-1-1), R ) 722 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 723 CALL SLARTGP( B12BULGE, B12D(I-1), RWORK(IV2TSN+I-1-1), 724 $ RWORK(IV2TCS+I-1-1), R ) 725 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 726 CALL SLARTGP( B22BULGE, B22D(I-1), RWORK(IV2TSN+I-1-1), 727 $ RWORK(IV2TCS+I-1-1), R ) 728 ELSE IF( NU .LT. MU ) THEN 729 CALL SLARTGS( B12E(I-1), B12D(I), NU, 730 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 731 ELSE 732 CALL SLARTGS( B22E(I-1), B22D(I), MU, 733 $ RWORK(IV2TCS+I-1-1), RWORK(IV2TSN+I-1-1) ) 734 END IF 735* 736 TEMP = RWORK(IV1TCS+I-1)*B11D(I) + RWORK(IV1TSN+I-1)*B11E(I) 737 B11E(I) = RWORK(IV1TCS+I-1)*B11E(I) - 738 $ RWORK(IV1TSN+I-1)*B11D(I) 739 B11D(I) = TEMP 740 B11BULGE = RWORK(IV1TSN+I-1)*B11D(I+1) 741 B11D(I+1) = RWORK(IV1TCS+I-1)*B11D(I+1) 742 TEMP = RWORK(IV1TCS+I-1)*B21D(I) + RWORK(IV1TSN+I-1)*B21E(I) 743 B21E(I) = RWORK(IV1TCS+I-1)*B21E(I) - 744 $ RWORK(IV1TSN+I-1)*B21D(I) 745 B21D(I) = TEMP 746 B21BULGE = RWORK(IV1TSN+I-1)*B21D(I+1) 747 B21D(I+1) = RWORK(IV1TCS+I-1)*B21D(I+1) 748 TEMP = RWORK(IV2TCS+I-1-1)*B12E(I-1) + 749 $ RWORK(IV2TSN+I-1-1)*B12D(I) 750 B12D(I) = RWORK(IV2TCS+I-1-1)*B12D(I) - 751 $ RWORK(IV2TSN+I-1-1)*B12E(I-1) 752 B12E(I-1) = TEMP 753 B12BULGE = RWORK(IV2TSN+I-1-1)*B12E(I) 754 B12E(I) = RWORK(IV2TCS+I-1-1)*B12E(I) 755 TEMP = RWORK(IV2TCS+I-1-1)*B22E(I-1) + 756 $ RWORK(IV2TSN+I-1-1)*B22D(I) 757 B22D(I) = RWORK(IV2TCS+I-1-1)*B22D(I) - 758 $ RWORK(IV2TSN+I-1-1)*B22E(I-1) 759 B22E(I-1) = TEMP 760 B22BULGE = RWORK(IV2TSN+I-1-1)*B22E(I) 761 B22E(I) = RWORK(IV2TCS+I-1-1)*B22E(I) 762* 763* Compute THETA(I) 764* 765 X1 = COS(PHI(I-1))*B11D(I) + SIN(PHI(I-1))*B12E(I-1) 766 X2 = COS(PHI(I-1))*B11BULGE + SIN(PHI(I-1))*B12BULGE 767 Y1 = COS(PHI(I-1))*B21D(I) + SIN(PHI(I-1))*B22E(I-1) 768 Y2 = COS(PHI(I-1))*B21BULGE + SIN(PHI(I-1))*B22BULGE 769* 770 THETA(I) = ATAN2( SQRT(Y1**2+Y2**2), SQRT(X1**2+X2**2) ) 771* 772* Determine if there are bulges to chase or if a new direct 773* summand has been reached 774* 775 RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2 776 RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2 777 RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2 778 RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2 779* 780* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1), 781* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge- 782* chasing by applying the original shift again. 783* 784 IF( .NOT. RESTART11 .AND. .NOT. RESTART12 ) THEN 785 CALL SLARTGP( X2, X1, RWORK(IU1SN+I-1), RWORK(IU1CS+I-1), 786 $ R ) 787 ELSE IF( .NOT. RESTART11 .AND. RESTART12 ) THEN 788 CALL SLARTGP( B11BULGE, B11D(I), RWORK(IU1SN+I-1), 789 $ RWORK(IU1CS+I-1), R ) 790 ELSE IF( RESTART11 .AND. .NOT. RESTART12 ) THEN 791 CALL SLARTGP( B12BULGE, B12E(I-1), RWORK(IU1SN+I-1), 792 $ RWORK(IU1CS+I-1), R ) 793 ELSE IF( MU .LE. NU ) THEN 794 CALL SLARTGS( B11E(I), B11D(I+1), MU, RWORK(IU1CS+I-1), 795 $ RWORK(IU1SN+I-1) ) 796 ELSE 797 CALL SLARTGS( B12D(I), B12E(I), NU, RWORK(IU1CS+I-1), 798 $ RWORK(IU1SN+I-1) ) 799 END IF 800 IF( .NOT. RESTART21 .AND. .NOT. RESTART22 ) THEN 801 CALL SLARTGP( Y2, Y1, RWORK(IU2SN+I-1), RWORK(IU2CS+I-1), 802 $ R ) 803 ELSE IF( .NOT. RESTART21 .AND. RESTART22 ) THEN 804 CALL SLARTGP( B21BULGE, B21D(I), RWORK(IU2SN+I-1), 805 $ RWORK(IU2CS+I-1), R ) 806 ELSE IF( RESTART21 .AND. .NOT. RESTART22 ) THEN 807 CALL SLARTGP( B22BULGE, B22E(I-1), RWORK(IU2SN+I-1), 808 $ RWORK(IU2CS+I-1), R ) 809 ELSE IF( NU .LT. MU ) THEN 810 CALL SLARTGS( B21E(I), B21E(I+1), NU, RWORK(IU2CS+I-1), 811 $ RWORK(IU2SN+I-1) ) 812 ELSE 813 CALL SLARTGS( B22D(I), B22E(I), MU, RWORK(IU2CS+I-1), 814 $ RWORK(IU2SN+I-1) ) 815 END IF 816 RWORK(IU2CS+I-1) = -RWORK(IU2CS+I-1) 817 RWORK(IU2SN+I-1) = -RWORK(IU2SN+I-1) 818* 819 TEMP = RWORK(IU1CS+I-1)*B11E(I) + RWORK(IU1SN+I-1)*B11D(I+1) 820 B11D(I+1) = RWORK(IU1CS+I-1)*B11D(I+1) - 821 $ RWORK(IU1SN+I-1)*B11E(I) 822 B11E(I) = TEMP 823 IF( I .LT. IMAX - 1 ) THEN 824 B11BULGE = RWORK(IU1SN+I-1)*B11E(I+1) 825 B11E(I+1) = RWORK(IU1CS+I-1)*B11E(I+1) 826 END IF 827 TEMP = RWORK(IU2CS+I-1)*B21E(I) + RWORK(IU2SN+I-1)*B21D(I+1) 828 B21D(I+1) = RWORK(IU2CS+I-1)*B21D(I+1) - 829 $ RWORK(IU2SN+I-1)*B21E(I) 830 B21E(I) = TEMP 831 IF( I .LT. IMAX - 1 ) THEN 832 B21BULGE = RWORK(IU2SN+I-1)*B21E(I+1) 833 B21E(I+1) = RWORK(IU2CS+I-1)*B21E(I+1) 834 END IF 835 TEMP = RWORK(IU1CS+I-1)*B12D(I) + RWORK(IU1SN+I-1)*B12E(I) 836 B12E(I) = RWORK(IU1CS+I-1)*B12E(I) - 837 $ RWORK(IU1SN+I-1)*B12D(I) 838 B12D(I) = TEMP 839 B12BULGE = RWORK(IU1SN+I-1)*B12D(I+1) 840 B12D(I+1) = RWORK(IU1CS+I-1)*B12D(I+1) 841 TEMP = RWORK(IU2CS+I-1)*B22D(I) + RWORK(IU2SN+I-1)*B22E(I) 842 B22E(I) = RWORK(IU2CS+I-1)*B22E(I) - 843 $ RWORK(IU2SN+I-1)*B22D(I) 844 B22D(I) = TEMP 845 B22BULGE = RWORK(IU2SN+I-1)*B22D(I+1) 846 B22D(I+1) = RWORK(IU2CS+I-1)*B22D(I+1) 847* 848 END DO 849* 850* Compute PHI(IMAX-1) 851* 852 X1 = SIN(THETA(IMAX-1))*B11E(IMAX-1) + 853 $ COS(THETA(IMAX-1))*B21E(IMAX-1) 854 Y1 = SIN(THETA(IMAX-1))*B12D(IMAX-1) + 855 $ COS(THETA(IMAX-1))*B22D(IMAX-1) 856 Y2 = SIN(THETA(IMAX-1))*B12BULGE + COS(THETA(IMAX-1))*B22BULGE 857* 858 PHI(IMAX-1) = ATAN2( ABS(X1), SQRT(Y1**2+Y2**2) ) 859* 860* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX) 861* 862 RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2 863 RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2 864* 865 IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN 866 CALL SLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1), 867 $ RWORK(IV2TCS+IMAX-1-1), R ) 868 ELSE IF( .NOT. RESTART12 .AND. RESTART22 ) THEN 869 CALL SLARTGP( B12BULGE, B12D(IMAX-1), 870 $ RWORK(IV2TSN+IMAX-1-1), 871 $ RWORK(IV2TCS+IMAX-1-1), R ) 872 ELSE IF( RESTART12 .AND. .NOT. RESTART22 ) THEN 873 CALL SLARTGP( B22BULGE, B22D(IMAX-1), 874 $ RWORK(IV2TSN+IMAX-1-1), 875 $ RWORK(IV2TCS+IMAX-1-1), R ) 876 ELSE IF( NU .LT. MU ) THEN 877 CALL SLARTGS( B12E(IMAX-1), B12D(IMAX), NU, 878 $ RWORK(IV2TCS+IMAX-1-1), 879 $ RWORK(IV2TSN+IMAX-1-1) ) 880 ELSE 881 CALL SLARTGS( B22E(IMAX-1), B22D(IMAX), MU, 882 $ RWORK(IV2TCS+IMAX-1-1), 883 $ RWORK(IV2TSN+IMAX-1-1) ) 884 END IF 885* 886 TEMP = RWORK(IV2TCS+IMAX-1-1)*B12E(IMAX-1) + 887 $ RWORK(IV2TSN+IMAX-1-1)*B12D(IMAX) 888 B12D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B12D(IMAX) - 889 $ RWORK(IV2TSN+IMAX-1-1)*B12E(IMAX-1) 890 B12E(IMAX-1) = TEMP 891 TEMP = RWORK(IV2TCS+IMAX-1-1)*B22E(IMAX-1) + 892 $ RWORK(IV2TSN+IMAX-1-1)*B22D(IMAX) 893 B22D(IMAX) = RWORK(IV2TCS+IMAX-1-1)*B22D(IMAX) - 894 $ RWORK(IV2TSN+IMAX-1-1)*B22E(IMAX-1) 895 B22E(IMAX-1) = TEMP 896* 897* Update singular vectors 898* 899 IF( WANTU1 ) THEN 900 IF( COLMAJOR ) THEN 901 CALL CLASR( 'R', 'V', 'F', P, IMAX-IMIN+1, 902 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 903 $ U1(1,IMIN), LDU1 ) 904 ELSE 905 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, P, 906 $ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1), 907 $ U1(IMIN,1), LDU1 ) 908 END IF 909 END IF 910 IF( WANTU2 ) THEN 911 IF( COLMAJOR ) THEN 912 CALL CLASR( 'R', 'V', 'F', M-P, IMAX-IMIN+1, 913 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 914 $ U2(1,IMIN), LDU2 ) 915 ELSE 916 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-P, 917 $ RWORK(IU2CS+IMIN-1), RWORK(IU2SN+IMIN-1), 918 $ U2(IMIN,1), LDU2 ) 919 END IF 920 END IF 921 IF( WANTV1T ) THEN 922 IF( COLMAJOR ) THEN 923 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, Q, 924 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 925 $ V1T(IMIN,1), LDV1T ) 926 ELSE 927 CALL CLASR( 'R', 'V', 'F', Q, IMAX-IMIN+1, 928 $ RWORK(IV1TCS+IMIN-1), RWORK(IV1TSN+IMIN-1), 929 $ V1T(1,IMIN), LDV1T ) 930 END IF 931 END IF 932 IF( WANTV2T ) THEN 933 IF( COLMAJOR ) THEN 934 CALL CLASR( 'L', 'V', 'F', IMAX-IMIN+1, M-Q, 935 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 936 $ V2T(IMIN,1), LDV2T ) 937 ELSE 938 CALL CLASR( 'R', 'V', 'F', M-Q, IMAX-IMIN+1, 939 $ RWORK(IV2TCS+IMIN-1), RWORK(IV2TSN+IMIN-1), 940 $ V2T(1,IMIN), LDV2T ) 941 END IF 942 END IF 943* 944* Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX) 945* 946 IF( B11E(IMAX-1)+B21E(IMAX-1) .GT. 0 ) THEN 947 B11D(IMAX) = -B11D(IMAX) 948 B21D(IMAX) = -B21D(IMAX) 949 IF( WANTV1T ) THEN 950 IF( COLMAJOR ) THEN 951 CALL CSCAL( Q, NEGONECOMPLEX, V1T(IMAX,1), LDV1T ) 952 ELSE 953 CALL CSCAL( Q, NEGONECOMPLEX, V1T(1,IMAX), 1 ) 954 END IF 955 END IF 956 END IF 957* 958* Compute THETA(IMAX) 959* 960 X1 = COS(PHI(IMAX-1))*B11D(IMAX) + 961 $ SIN(PHI(IMAX-1))*B12E(IMAX-1) 962 Y1 = COS(PHI(IMAX-1))*B21D(IMAX) + 963 $ SIN(PHI(IMAX-1))*B22E(IMAX-1) 964* 965 THETA(IMAX) = ATAN2( ABS(Y1), ABS(X1) ) 966* 967* Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX), 968* and B22(IMAX,IMAX-1) 969* 970 IF( B11D(IMAX)+B12E(IMAX-1) .LT. 0 ) THEN 971 B12D(IMAX) = -B12D(IMAX) 972 IF( WANTU1 ) THEN 973 IF( COLMAJOR ) THEN 974 CALL CSCAL( P, NEGONECOMPLEX, U1(1,IMAX), 1 ) 975 ELSE 976 CALL CSCAL( P, NEGONECOMPLEX, U1(IMAX,1), LDU1 ) 977 END IF 978 END IF 979 END IF 980 IF( B21D(IMAX)+B22E(IMAX-1) .GT. 0 ) THEN 981 B22D(IMAX) = -B22D(IMAX) 982 IF( WANTU2 ) THEN 983 IF( COLMAJOR ) THEN 984 CALL CSCAL( M-P, NEGONECOMPLEX, U2(1,IMAX), 1 ) 985 ELSE 986 CALL CSCAL( M-P, NEGONECOMPLEX, U2(IMAX,1), LDU2 ) 987 END IF 988 END IF 989 END IF 990* 991* Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX) 992* 993 IF( B12D(IMAX)+B22D(IMAX) .LT. 0 ) THEN 994 IF( WANTV2T ) THEN 995 IF( COLMAJOR ) THEN 996 CALL CSCAL( M-Q, NEGONECOMPLEX, V2T(IMAX,1), LDV2T ) 997 ELSE 998 CALL CSCAL( M-Q, NEGONECOMPLEX, V2T(1,IMAX), 1 ) 999 END IF 1000 END IF 1001 END IF 1002* 1003* Test for negligible sines or cosines 1004* 1005 DO I = IMIN, IMAX 1006 IF( THETA(I) .LT. THRESH ) THEN 1007 THETA(I) = ZERO 1008 ELSE IF( THETA(I) .GT. PIOVER2-THRESH ) THEN 1009 THETA(I) = PIOVER2 1010 END IF 1011 END DO 1012 DO I = IMIN, IMAX-1 1013 IF( PHI(I) .LT. THRESH ) THEN 1014 PHI(I) = ZERO 1015 ELSE IF( PHI(I) .GT. PIOVER2-THRESH ) THEN 1016 PHI(I) = PIOVER2 1017 END IF 1018 END DO 1019* 1020* Deflate 1021* 1022 IF (IMAX .GT. 1) THEN 1023 DO WHILE( PHI(IMAX-1) .EQ. ZERO ) 1024 IMAX = IMAX - 1 1025 IF (IMAX .LE. 1) EXIT 1026 END DO 1027 END IF 1028 IF( IMIN .GT. IMAX - 1 ) 1029 $ IMIN = IMAX - 1 1030 IF (IMIN .GT. 1) THEN 1031 DO WHILE (PHI(IMIN-1) .NE. ZERO) 1032 IMIN = IMIN - 1 1033 IF (IMIN .LE. 1) EXIT 1034 END DO 1035 END IF 1036* 1037* Repeat main iteration loop 1038* 1039 END DO 1040* 1041* Postprocessing: order THETA from least to greatest 1042* 1043 DO I = 1, Q 1044* 1045 MINI = I 1046 THETAMIN = THETA(I) 1047 DO J = I+1, Q 1048 IF( THETA(J) .LT. THETAMIN ) THEN 1049 MINI = J 1050 THETAMIN = THETA(J) 1051 END IF 1052 END DO 1053* 1054 IF( MINI .NE. I ) THEN 1055 THETA(MINI) = THETA(I) 1056 THETA(I) = THETAMIN 1057 IF( COLMAJOR ) THEN 1058 IF( WANTU1 ) 1059 $ CALL CSWAP( P, U1(1,I), 1, U1(1,MINI), 1 ) 1060 IF( WANTU2 ) 1061 $ CALL CSWAP( M-P, U2(1,I), 1, U2(1,MINI), 1 ) 1062 IF( WANTV1T ) 1063 $ CALL CSWAP( Q, V1T(I,1), LDV1T, V1T(MINI,1), LDV1T ) 1064 IF( WANTV2T ) 1065 $ CALL CSWAP( M-Q, V2T(I,1), LDV2T, V2T(MINI,1), 1066 $ LDV2T ) 1067 ELSE 1068 IF( WANTU1 ) 1069 $ CALL CSWAP( P, U1(I,1), LDU1, U1(MINI,1), LDU1 ) 1070 IF( WANTU2 ) 1071 $ CALL CSWAP( M-P, U2(I,1), LDU2, U2(MINI,1), LDU2 ) 1072 IF( WANTV1T ) 1073 $ CALL CSWAP( Q, V1T(1,I), 1, V1T(1,MINI), 1 ) 1074 IF( WANTV2T ) 1075 $ CALL CSWAP( M-Q, V2T(1,I), 1, V2T(1,MINI), 1 ) 1076 END IF 1077 END IF 1078* 1079 END DO 1080* 1081 RETURN 1082* 1083* End of CBBCSD 1084* 1085 END 1086 1087