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