1*> \brief \b DGSVJ0 pre-processor for the routine dgesvj. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DGSVJ0 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgsvj0.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgsvj0.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgsvj0.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, 22* SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP 26* DOUBLE PRECISION EPS, SFMIN, TOL 27* CHARACTER*1 JOBV 28* .. 29* .. Array Arguments .. 30* DOUBLE PRECISION A( LDA, * ), SVA( N ), D( N ), V( LDV, * ), 31* $ WORK( LWORK ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> DGSVJ0 is called from DGESVJ as a pre-processor and that is its main 41*> purpose. It applies Jacobi rotations in the same way as DGESVJ does, but 42*> it does not check convergence (stopping criterion). Few tuning 43*> parameters (marked by [TP]) are available for the implementer. 44*> \endverbatim 45* 46* Arguments: 47* ========== 48* 49*> \param[in] JOBV 50*> \verbatim 51*> JOBV is CHARACTER*1 52*> Specifies whether the output from this procedure is used 53*> to compute the matrix V: 54*> = 'V': the product of the Jacobi rotations is accumulated 55*> by postmulyiplying the N-by-N array V. 56*> (See the description of V.) 57*> = 'A': the product of the Jacobi rotations is accumulated 58*> by postmulyiplying the MV-by-N array V. 59*> (See the descriptions of MV and V.) 60*> = 'N': the Jacobi rotations are not accumulated. 61*> \endverbatim 62*> 63*> \param[in] M 64*> \verbatim 65*> M is INTEGER 66*> The number of rows of the input matrix A. M >= 0. 67*> \endverbatim 68*> 69*> \param[in] N 70*> \verbatim 71*> N is INTEGER 72*> The number of columns of the input matrix A. 73*> M >= N >= 0. 74*> \endverbatim 75*> 76*> \param[in,out] A 77*> \verbatim 78*> A is DOUBLE PRECISION array, dimension (LDA,N) 79*> On entry, M-by-N matrix A, such that A*diag(D) represents 80*> the input matrix. 81*> On exit, 82*> A_onexit * D_onexit represents the input matrix A*diag(D) 83*> post-multiplied by a sequence of Jacobi rotations, where the 84*> rotation threshold and the total number of sweeps are given in 85*> TOL and NSWEEP, respectively. 86*> (See the descriptions of D, TOL and NSWEEP.) 87*> \endverbatim 88*> 89*> \param[in] LDA 90*> \verbatim 91*> LDA is INTEGER 92*> The leading dimension of the array A. LDA >= max(1,M). 93*> \endverbatim 94*> 95*> \param[in,out] D 96*> \verbatim 97*> D is DOUBLE PRECISION array, dimension (N) 98*> The array D accumulates the scaling factors from the fast scaled 99*> Jacobi rotations. 100*> On entry, A*diag(D) represents the input matrix. 101*> On exit, A_onexit*diag(D_onexit) represents the input matrix 102*> post-multiplied by a sequence of Jacobi rotations, where the 103*> rotation threshold and the total number of sweeps are given in 104*> TOL and NSWEEP, respectively. 105*> (See the descriptions of A, TOL and NSWEEP.) 106*> \endverbatim 107*> 108*> \param[in,out] SVA 109*> \verbatim 110*> SVA is DOUBLE PRECISION array, dimension (N) 111*> On entry, SVA contains the Euclidean norms of the columns of 112*> the matrix A*diag(D). 113*> On exit, SVA contains the Euclidean norms of the columns of 114*> the matrix onexit*diag(D_onexit). 115*> \endverbatim 116*> 117*> \param[in] MV 118*> \verbatim 119*> MV is INTEGER 120*> If JOBV = 'A', then MV rows of V are post-multipled by a 121*> sequence of Jacobi rotations. 122*> If JOBV = 'N', then MV is not referenced. 123*> \endverbatim 124*> 125*> \param[in,out] V 126*> \verbatim 127*> V is DOUBLE PRECISION array, dimension (LDV,N) 128*> If JOBV = 'V' then N rows of V are post-multipled by a 129*> sequence of Jacobi rotations. 130*> If JOBV = 'A' then MV rows of V are post-multipled by a 131*> sequence of Jacobi rotations. 132*> If JOBV = 'N', then V is not referenced. 133*> \endverbatim 134*> 135*> \param[in] LDV 136*> \verbatim 137*> LDV is INTEGER 138*> The leading dimension of the array V, LDV >= 1. 139*> If JOBV = 'V', LDV >= N. 140*> If JOBV = 'A', LDV >= MV. 141*> \endverbatim 142*> 143*> \param[in] EPS 144*> \verbatim 145*> EPS is DOUBLE PRECISION 146*> EPS = DLAMCH('Epsilon') 147*> \endverbatim 148*> 149*> \param[in] SFMIN 150*> \verbatim 151*> SFMIN is DOUBLE PRECISION 152*> SFMIN = DLAMCH('Safe Minimum') 153*> \endverbatim 154*> 155*> \param[in] TOL 156*> \verbatim 157*> TOL is DOUBLE PRECISION 158*> TOL is the threshold for Jacobi rotations. For a pair 159*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is 160*> applied only if DABS(COS(angle(A(:,p),A(:,q)))) > TOL. 161*> \endverbatim 162*> 163*> \param[in] NSWEEP 164*> \verbatim 165*> NSWEEP is INTEGER 166*> NSWEEP is the number of sweeps of Jacobi rotations to be 167*> performed. 168*> \endverbatim 169*> 170*> \param[out] WORK 171*> \verbatim 172*> WORK is DOUBLE PRECISION array, dimension (LWORK) 173*> \endverbatim 174*> 175*> \param[in] LWORK 176*> \verbatim 177*> LWORK is INTEGER 178*> LWORK is the dimension of WORK. LWORK >= M. 179*> \endverbatim 180*> 181*> \param[out] INFO 182*> \verbatim 183*> INFO is INTEGER 184*> = 0: successful exit. 185*> < 0: if INFO = -i, then the i-th argument had an illegal value 186*> \endverbatim 187* 188* Authors: 189* ======== 190* 191*> \author Univ. of Tennessee 192*> \author Univ. of California Berkeley 193*> \author Univ. of Colorado Denver 194*> \author NAG Ltd. 195* 196*> \ingroup doubleOTHERcomputational 197* 198*> \par Further Details: 199* ===================== 200*> 201*> DGSVJ0 is used just to enable DGESVJ to call a simplified version of 202*> itself to work on a submatrix of the original matrix. 203*> 204*> \par Contributors: 205* ================== 206*> 207*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) 208*> 209*> \par Bugs, Examples and Comments: 210* ================================= 211*> 212*> Please report all bugs and send interesting test examples and comments to 213*> drmac@math.hr. Thank you. 214* 215* ===================================================================== 216 SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, 217 $ SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) 218* 219* -- LAPACK computational routine -- 220* -- LAPACK is a software package provided by Univ. of Tennessee, -- 221* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 222* 223* .. Scalar Arguments .. 224 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP 225 DOUBLE PRECISION EPS, SFMIN, TOL 226 CHARACTER*1 JOBV 227* .. 228* .. Array Arguments .. 229 DOUBLE PRECISION A( LDA, * ), SVA( N ), D( N ), V( LDV, * ), 230 $ WORK( LWORK ) 231* .. 232* 233* ===================================================================== 234* 235* .. Local Parameters .. 236 DOUBLE PRECISION ZERO, HALF, ONE 237 PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) 238* .. 239* .. Local Scalars .. 240 DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, 241 $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS, 242 $ ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA, 243 $ THSIGN 244 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, 245 $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL, 246 $ NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND 247 LOGICAL APPLV, ROTOK, RSVEC 248* .. 249* .. Local Arrays .. 250 DOUBLE PRECISION FASTR( 5 ) 251* .. 252* .. Intrinsic Functions .. 253 INTRINSIC DABS, MAX, DBLE, MIN, DSIGN, DSQRT 254* .. 255* .. External Functions .. 256 DOUBLE PRECISION DDOT, DNRM2 257 INTEGER IDAMAX 258 LOGICAL LSAME 259 EXTERNAL IDAMAX, LSAME, DDOT, DNRM2 260* .. 261* .. External Subroutines .. 262 EXTERNAL DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP, 263 $ XERBLA 264* .. 265* .. Executable Statements .. 266* 267* Test the input parameters. 268* 269 APPLV = LSAME( JOBV, 'A' ) 270 RSVEC = LSAME( JOBV, 'V' ) 271 IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN 272 INFO = -1 273 ELSE IF( M.LT.0 ) THEN 274 INFO = -2 275 ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN 276 INFO = -3 277 ELSE IF( LDA.LT.M ) THEN 278 INFO = -5 279 ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN 280 INFO = -8 281 ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 282 $ ( APPLV.AND.( LDV.LT.MV ) ) ) THEN 283 INFO = -10 284 ELSE IF( TOL.LE.EPS ) THEN 285 INFO = -13 286 ELSE IF( NSWEEP.LT.0 ) THEN 287 INFO = -14 288 ELSE IF( LWORK.LT.M ) THEN 289 INFO = -16 290 ELSE 291 INFO = 0 292 END IF 293* 294* #:( 295 IF( INFO.NE.0 ) THEN 296 CALL XERBLA( 'DGSVJ0', -INFO ) 297 RETURN 298 END IF 299* 300 IF( RSVEC ) THEN 301 MVL = N 302 ELSE IF( APPLV ) THEN 303 MVL = MV 304 END IF 305 RSVEC = RSVEC .OR. APPLV 306 307 ROOTEPS = DSQRT( EPS ) 308 ROOTSFMIN = DSQRT( SFMIN ) 309 SMALL = SFMIN / EPS 310 BIG = ONE / SFMIN 311 ROOTBIG = ONE / ROOTSFMIN 312 BIGTHETA = ONE / ROOTEPS 313 ROOTTOL = DSQRT( TOL ) 314* 315* -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#- 316* 317 EMPTSW = ( N*( N-1 ) ) / 2 318 NOTROT = 0 319 FASTR( 1 ) = ZERO 320* 321* -#- Row-cyclic pivot strategy with de Rijk's pivoting -#- 322* 323 324 SWBAND = 0 325*[TP] SWBAND is a tuning parameter. It is meaningful and effective 326* if SGESVJ is used as a computational routine in the preconditioned 327* Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure 328* ...... 329 330 KBL = MIN( 8, N ) 331*[TP] KBL is a tuning parameter that defines the tile size in the 332* tiling of the p-q loops of pivot pairs. In general, an optimal 333* value of KBL depends on the matrix dimensions and on the 334* parameters of the computer's memory. 335* 336 NBL = N / KBL 337 IF( ( NBL*KBL ).NE.N )NBL = NBL + 1 338 339 BLSKIP = ( KBL**2 ) + 1 340*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL. 341 342 ROWSKIP = MIN( 5, KBL ) 343*[TP] ROWSKIP is a tuning parameter. 344 345 LKAHEAD = 1 346*[TP] LKAHEAD is a tuning parameter. 347 SWBAND = 0 348 PSKIPPED = 0 349* 350 DO 1993 i = 1, NSWEEP 351* .. go go go ... 352* 353 MXAAPQ = ZERO 354 MXSINJ = ZERO 355 ISWROT = 0 356* 357 NOTROT = 0 358 PSKIPPED = 0 359* 360 DO 2000 ibr = 1, NBL 361 362 igl = ( ibr-1 )*KBL + 1 363* 364 DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr ) 365* 366 igl = igl + ir1*KBL 367* 368 DO 2001 p = igl, MIN( igl+KBL-1, N-1 ) 369 370* .. de Rijk's pivoting 371 q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1 372 IF( p.NE.q ) THEN 373 CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 374 IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, 375 $ V( 1, q ), 1 ) 376 TEMP1 = SVA( p ) 377 SVA( p ) = SVA( q ) 378 SVA( q ) = TEMP1 379 TEMP1 = D( p ) 380 D( p ) = D( q ) 381 D( q ) = TEMP1 382 END IF 383* 384 IF( ir1.EQ.0 ) THEN 385* 386* Column norms are periodically updated by explicit 387* norm computation. 388* Caveat: 389* Some BLAS implementations compute DNRM2(M,A(1,p),1) 390* as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may result in 391* overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and 392* underflow for ||A(:,p)||_2 < DSQRT(underflow_threshold). 393* Hence, DNRM2 cannot be trusted, not even in the case when 394* the true norm is far from the under(over)flow boundaries. 395* If properly implemented DNRM2 is available, the IF-THEN-ELSE 396* below should read "AAPP = DNRM2( M, A(1,p), 1 ) * D(p)". 397* 398 IF( ( SVA( p ).LT.ROOTBIG ) .AND. 399 $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN 400 SVA( p ) = DNRM2( M, A( 1, p ), 1 )*D( p ) 401 ELSE 402 TEMP1 = ZERO 403 AAPP = ONE 404 CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) 405 SVA( p ) = TEMP1*DSQRT( AAPP )*D( p ) 406 END IF 407 AAPP = SVA( p ) 408 ELSE 409 AAPP = SVA( p ) 410 END IF 411 412* 413 IF( AAPP.GT.ZERO ) THEN 414* 415 PSKIPPED = 0 416* 417 DO 2002 q = p + 1, MIN( igl+KBL-1, N ) 418* 419 AAQQ = SVA( q ) 420 421 IF( AAQQ.GT.ZERO ) THEN 422* 423 AAPP0 = AAPP 424 IF( AAQQ.GE.ONE ) THEN 425 ROTOK = ( SMALL*AAPP ).LE.AAQQ 426 IF( AAPP.LT.( BIG / AAQQ ) ) THEN 427 AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, 428 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 429 $ / AAPP 430 ELSE 431 CALL DCOPY( M, A( 1, p ), 1, WORK, 1 ) 432 CALL DLASCL( 'G', 0, 0, AAPP, D( p ), 433 $ M, 1, WORK, LDA, IERR ) 434 AAPQ = DDOT( M, WORK, 1, A( 1, q ), 435 $ 1 )*D( q ) / AAQQ 436 END IF 437 ELSE 438 ROTOK = AAPP.LE.( AAQQ / SMALL ) 439 IF( AAPP.GT.( SMALL / AAQQ ) ) THEN 440 AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, 441 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 442 $ / AAPP 443 ELSE 444 CALL DCOPY( M, A( 1, q ), 1, WORK, 1 ) 445 CALL DLASCL( 'G', 0, 0, AAQQ, D( q ), 446 $ M, 1, WORK, LDA, IERR ) 447 AAPQ = DDOT( M, WORK, 1, A( 1, p ), 448 $ 1 )*D( p ) / AAPP 449 END IF 450 END IF 451* 452 MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) ) 453* 454* TO rotate or NOT to rotate, THAT is the question ... 455* 456 IF( DABS( AAPQ ).GT.TOL ) THEN 457* 458* .. rotate 459* ROTATED = ROTATED + ONE 460* 461 IF( ir1.EQ.0 ) THEN 462 NOTROT = 0 463 PSKIPPED = 0 464 ISWROT = ISWROT + 1 465 END IF 466* 467 IF( ROTOK ) THEN 468* 469 AQOAP = AAQQ / AAPP 470 APOAQ = AAPP / AAQQ 471 THETA = -HALF*DABS( AQOAP-APOAQ )/AAPQ 472* 473 IF( DABS( THETA ).GT.BIGTHETA ) THEN 474* 475 T = HALF / THETA 476 FASTR( 3 ) = T*D( p ) / D( q ) 477 FASTR( 4 ) = -T*D( q ) / D( p ) 478 CALL DROTM( M, A( 1, p ), 1, 479 $ A( 1, q ), 1, FASTR ) 480 IF( RSVEC )CALL DROTM( MVL, 481 $ V( 1, p ), 1, 482 $ V( 1, q ), 1, 483 $ FASTR ) 484 SVA( q ) = AAQQ*DSQRT( MAX( ZERO, 485 $ ONE+T*APOAQ*AAPQ ) ) 486 AAPP = AAPP*DSQRT( MAX( ZERO, 487 $ ONE-T*AQOAP*AAPQ ) ) 488 MXSINJ = MAX( MXSINJ, DABS( T ) ) 489* 490 ELSE 491* 492* .. choose correct signum for THETA and rotate 493* 494 THSIGN = -DSIGN( ONE, AAPQ ) 495 T = ONE / ( THETA+THSIGN* 496 $ DSQRT( ONE+THETA*THETA ) ) 497 CS = DSQRT( ONE / ( ONE+T*T ) ) 498 SN = T*CS 499* 500 MXSINJ = MAX( MXSINJ, DABS( SN ) ) 501 SVA( q ) = AAQQ*DSQRT( MAX( ZERO, 502 $ ONE+T*APOAQ*AAPQ ) ) 503 AAPP = AAPP*DSQRT( MAX( ZERO, 504 $ ONE-T*AQOAP*AAPQ ) ) 505* 506 APOAQ = D( p ) / D( q ) 507 AQOAP = D( q ) / D( p ) 508 IF( D( p ).GE.ONE ) THEN 509 IF( D( q ).GE.ONE ) THEN 510 FASTR( 3 ) = T*APOAQ 511 FASTR( 4 ) = -T*AQOAP 512 D( p ) = D( p )*CS 513 D( q ) = D( q )*CS 514 CALL DROTM( M, A( 1, p ), 1, 515 $ A( 1, q ), 1, 516 $ FASTR ) 517 IF( RSVEC )CALL DROTM( MVL, 518 $ V( 1, p ), 1, V( 1, q ), 519 $ 1, FASTR ) 520 ELSE 521 CALL DAXPY( M, -T*AQOAP, 522 $ A( 1, q ), 1, 523 $ A( 1, p ), 1 ) 524 CALL DAXPY( M, CS*SN*APOAQ, 525 $ A( 1, p ), 1, 526 $ A( 1, q ), 1 ) 527 D( p ) = D( p )*CS 528 D( q ) = D( q ) / CS 529 IF( RSVEC ) THEN 530 CALL DAXPY( MVL, -T*AQOAP, 531 $ V( 1, q ), 1, 532 $ V( 1, p ), 1 ) 533 CALL DAXPY( MVL, 534 $ CS*SN*APOAQ, 535 $ V( 1, p ), 1, 536 $ V( 1, q ), 1 ) 537 END IF 538 END IF 539 ELSE 540 IF( D( q ).GE.ONE ) THEN 541 CALL DAXPY( M, T*APOAQ, 542 $ A( 1, p ), 1, 543 $ A( 1, q ), 1 ) 544 CALL DAXPY( M, -CS*SN*AQOAP, 545 $ A( 1, q ), 1, 546 $ A( 1, p ), 1 ) 547 D( p ) = D( p ) / CS 548 D( q ) = D( q )*CS 549 IF( RSVEC ) THEN 550 CALL DAXPY( MVL, T*APOAQ, 551 $ V( 1, p ), 1, 552 $ V( 1, q ), 1 ) 553 CALL DAXPY( MVL, 554 $ -CS*SN*AQOAP, 555 $ V( 1, q ), 1, 556 $ V( 1, p ), 1 ) 557 END IF 558 ELSE 559 IF( D( p ).GE.D( q ) ) THEN 560 CALL DAXPY( M, -T*AQOAP, 561 $ A( 1, q ), 1, 562 $ A( 1, p ), 1 ) 563 CALL DAXPY( M, CS*SN*APOAQ, 564 $ A( 1, p ), 1, 565 $ A( 1, q ), 1 ) 566 D( p ) = D( p )*CS 567 D( q ) = D( q ) / CS 568 IF( RSVEC ) THEN 569 CALL DAXPY( MVL, 570 $ -T*AQOAP, 571 $ V( 1, q ), 1, 572 $ V( 1, p ), 1 ) 573 CALL DAXPY( MVL, 574 $ CS*SN*APOAQ, 575 $ V( 1, p ), 1, 576 $ V( 1, q ), 1 ) 577 END IF 578 ELSE 579 CALL DAXPY( M, T*APOAQ, 580 $ A( 1, p ), 1, 581 $ A( 1, q ), 1 ) 582 CALL DAXPY( M, 583 $ -CS*SN*AQOAP, 584 $ A( 1, q ), 1, 585 $ A( 1, p ), 1 ) 586 D( p ) = D( p ) / CS 587 D( q ) = D( q )*CS 588 IF( RSVEC ) THEN 589 CALL DAXPY( MVL, 590 $ T*APOAQ, V( 1, p ), 591 $ 1, V( 1, q ), 1 ) 592 CALL DAXPY( MVL, 593 $ -CS*SN*AQOAP, 594 $ V( 1, q ), 1, 595 $ V( 1, p ), 1 ) 596 END IF 597 END IF 598 END IF 599 END IF 600 END IF 601* 602 ELSE 603* .. have to use modified Gram-Schmidt like transformation 604 CALL DCOPY( M, A( 1, p ), 1, WORK, 1 ) 605 CALL DLASCL( 'G', 0, 0, AAPP, ONE, M, 606 $ 1, WORK, LDA, IERR ) 607 CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M, 608 $ 1, A( 1, q ), LDA, IERR ) 609 TEMP1 = -AAPQ*D( p ) / D( q ) 610 CALL DAXPY( M, TEMP1, WORK, 1, 611 $ A( 1, q ), 1 ) 612 CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M, 613 $ 1, A( 1, q ), LDA, IERR ) 614 SVA( q ) = AAQQ*DSQRT( MAX( ZERO, 615 $ ONE-AAPQ*AAPQ ) ) 616 MXSINJ = MAX( MXSINJ, SFMIN ) 617 END IF 618* END IF ROTOK THEN ... ELSE 619* 620* In the case of cancellation in updating SVA(q), SVA(p) 621* recompute SVA(q), SVA(p). 622 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 623 $ THEN 624 IF( ( AAQQ.LT.ROOTBIG ) .AND. 625 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN 626 SVA( q ) = DNRM2( M, A( 1, q ), 1 )* 627 $ D( q ) 628 ELSE 629 T = ZERO 630 AAQQ = ONE 631 CALL DLASSQ( M, A( 1, q ), 1, T, 632 $ AAQQ ) 633 SVA( q ) = T*DSQRT( AAQQ )*D( q ) 634 END IF 635 END IF 636 IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN 637 IF( ( AAPP.LT.ROOTBIG ) .AND. 638 $ ( AAPP.GT.ROOTSFMIN ) ) THEN 639 AAPP = DNRM2( M, A( 1, p ), 1 )* 640 $ D( p ) 641 ELSE 642 T = ZERO 643 AAPP = ONE 644 CALL DLASSQ( M, A( 1, p ), 1, T, 645 $ AAPP ) 646 AAPP = T*DSQRT( AAPP )*D( p ) 647 END IF 648 SVA( p ) = AAPP 649 END IF 650* 651 ELSE 652* A(:,p) and A(:,q) already numerically orthogonal 653 IF( ir1.EQ.0 )NOTROT = NOTROT + 1 654 PSKIPPED = PSKIPPED + 1 655 END IF 656 ELSE 657* A(:,q) is zero column 658 IF( ir1.EQ.0 )NOTROT = NOTROT + 1 659 PSKIPPED = PSKIPPED + 1 660 END IF 661* 662 IF( ( i.LE.SWBAND ) .AND. 663 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN 664 IF( ir1.EQ.0 )AAPP = -AAPP 665 NOTROT = 0 666 GO TO 2103 667 END IF 668* 669 2002 CONTINUE 670* END q-LOOP 671* 672 2103 CONTINUE 673* bailed out of q-loop 674 675 SVA( p ) = AAPP 676 677 ELSE 678 SVA( p ) = AAPP 679 IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) ) 680 $ NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p 681 END IF 682* 683 2001 CONTINUE 684* end of the p-loop 685* end of doing the block ( ibr, ibr ) 686 1002 CONTINUE 687* end of ir1-loop 688* 689*........................................................ 690* ... go to the off diagonal blocks 691* 692 igl = ( ibr-1 )*KBL + 1 693* 694 DO 2010 jbc = ibr + 1, NBL 695* 696 jgl = ( jbc-1 )*KBL + 1 697* 698* doing the block at ( ibr, jbc ) 699* 700 IJBLSK = 0 701 DO 2100 p = igl, MIN( igl+KBL-1, N ) 702* 703 AAPP = SVA( p ) 704* 705 IF( AAPP.GT.ZERO ) THEN 706* 707 PSKIPPED = 0 708* 709 DO 2200 q = jgl, MIN( jgl+KBL-1, N ) 710* 711 AAQQ = SVA( q ) 712* 713 IF( AAQQ.GT.ZERO ) THEN 714 AAPP0 = AAPP 715* 716* -#- M x 2 Jacobi SVD -#- 717* 718* -#- Safe Gram matrix computation -#- 719* 720 IF( AAQQ.GE.ONE ) THEN 721 IF( AAPP.GE.AAQQ ) THEN 722 ROTOK = ( SMALL*AAPP ).LE.AAQQ 723 ELSE 724 ROTOK = ( SMALL*AAQQ ).LE.AAPP 725 END IF 726 IF( AAPP.LT.( BIG / AAQQ ) ) THEN 727 AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, 728 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 729 $ / AAPP 730 ELSE 731 CALL DCOPY( M, A( 1, p ), 1, WORK, 1 ) 732 CALL DLASCL( 'G', 0, 0, AAPP, D( p ), 733 $ M, 1, WORK, LDA, IERR ) 734 AAPQ = DDOT( M, WORK, 1, A( 1, q ), 735 $ 1 )*D( q ) / AAQQ 736 END IF 737 ELSE 738 IF( AAPP.GE.AAQQ ) THEN 739 ROTOK = AAPP.LE.( AAQQ / SMALL ) 740 ELSE 741 ROTOK = AAQQ.LE.( AAPP / SMALL ) 742 END IF 743 IF( AAPP.GT.( SMALL / AAQQ ) ) THEN 744 AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, 745 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 746 $ / AAPP 747 ELSE 748 CALL DCOPY( M, A( 1, q ), 1, WORK, 1 ) 749 CALL DLASCL( 'G', 0, 0, AAQQ, D( q ), 750 $ M, 1, WORK, LDA, IERR ) 751 AAPQ = DDOT( M, WORK, 1, A( 1, p ), 752 $ 1 )*D( p ) / AAPP 753 END IF 754 END IF 755* 756 MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) ) 757* 758* TO rotate or NOT to rotate, THAT is the question ... 759* 760 IF( DABS( AAPQ ).GT.TOL ) THEN 761 NOTROT = 0 762* ROTATED = ROTATED + 1 763 PSKIPPED = 0 764 ISWROT = ISWROT + 1 765* 766 IF( ROTOK ) THEN 767* 768 AQOAP = AAQQ / AAPP 769 APOAQ = AAPP / AAQQ 770 THETA = -HALF*DABS( AQOAP-APOAQ )/AAPQ 771 IF( AAQQ.GT.AAPP0 )THETA = -THETA 772* 773 IF( DABS( THETA ).GT.BIGTHETA ) THEN 774 T = HALF / THETA 775 FASTR( 3 ) = T*D( p ) / D( q ) 776 FASTR( 4 ) = -T*D( q ) / D( p ) 777 CALL DROTM( M, A( 1, p ), 1, 778 $ A( 1, q ), 1, FASTR ) 779 IF( RSVEC )CALL DROTM( MVL, 780 $ V( 1, p ), 1, 781 $ V( 1, q ), 1, 782 $ FASTR ) 783 SVA( q ) = AAQQ*DSQRT( MAX( ZERO, 784 $ ONE+T*APOAQ*AAPQ ) ) 785 AAPP = AAPP*DSQRT( MAX( ZERO, 786 $ ONE-T*AQOAP*AAPQ ) ) 787 MXSINJ = MAX( MXSINJ, DABS( T ) ) 788 ELSE 789* 790* .. choose correct signum for THETA and rotate 791* 792 THSIGN = -DSIGN( ONE, AAPQ ) 793 IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN 794 T = ONE / ( THETA+THSIGN* 795 $ DSQRT( ONE+THETA*THETA ) ) 796 CS = DSQRT( ONE / ( ONE+T*T ) ) 797 SN = T*CS 798 MXSINJ = MAX( MXSINJ, DABS( SN ) ) 799 SVA( q ) = AAQQ*DSQRT( MAX( ZERO, 800 $ ONE+T*APOAQ*AAPQ ) ) 801 AAPP = AAPP*DSQRT( MAX( ZERO, 802 $ ONE-T*AQOAP*AAPQ ) ) 803* 804 APOAQ = D( p ) / D( q ) 805 AQOAP = D( q ) / D( p ) 806 IF( D( p ).GE.ONE ) THEN 807* 808 IF( D( q ).GE.ONE ) THEN 809 FASTR( 3 ) = T*APOAQ 810 FASTR( 4 ) = -T*AQOAP 811 D( p ) = D( p )*CS 812 D( q ) = D( q )*CS 813 CALL DROTM( M, A( 1, p ), 1, 814 $ A( 1, q ), 1, 815 $ FASTR ) 816 IF( RSVEC )CALL DROTM( MVL, 817 $ V( 1, p ), 1, V( 1, q ), 818 $ 1, FASTR ) 819 ELSE 820 CALL DAXPY( M, -T*AQOAP, 821 $ A( 1, q ), 1, 822 $ A( 1, p ), 1 ) 823 CALL DAXPY( M, CS*SN*APOAQ, 824 $ A( 1, p ), 1, 825 $ A( 1, q ), 1 ) 826 IF( RSVEC ) THEN 827 CALL DAXPY( MVL, -T*AQOAP, 828 $ V( 1, q ), 1, 829 $ V( 1, p ), 1 ) 830 CALL DAXPY( MVL, 831 $ CS*SN*APOAQ, 832 $ V( 1, p ), 1, 833 $ V( 1, q ), 1 ) 834 END IF 835 D( p ) = D( p )*CS 836 D( q ) = D( q ) / CS 837 END IF 838 ELSE 839 IF( D( q ).GE.ONE ) THEN 840 CALL DAXPY( M, T*APOAQ, 841 $ A( 1, p ), 1, 842 $ A( 1, q ), 1 ) 843 CALL DAXPY( M, -CS*SN*AQOAP, 844 $ A( 1, q ), 1, 845 $ A( 1, p ), 1 ) 846 IF( RSVEC ) THEN 847 CALL DAXPY( MVL, T*APOAQ, 848 $ V( 1, p ), 1, 849 $ V( 1, q ), 1 ) 850 CALL DAXPY( MVL, 851 $ -CS*SN*AQOAP, 852 $ V( 1, q ), 1, 853 $ V( 1, p ), 1 ) 854 END IF 855 D( p ) = D( p ) / CS 856 D( q ) = D( q )*CS 857 ELSE 858 IF( D( p ).GE.D( q ) ) THEN 859 CALL DAXPY( M, -T*AQOAP, 860 $ A( 1, q ), 1, 861 $ A( 1, p ), 1 ) 862 CALL DAXPY( M, CS*SN*APOAQ, 863 $ A( 1, p ), 1, 864 $ A( 1, q ), 1 ) 865 D( p ) = D( p )*CS 866 D( q ) = D( q ) / CS 867 IF( RSVEC ) THEN 868 CALL DAXPY( MVL, 869 $ -T*AQOAP, 870 $ V( 1, q ), 1, 871 $ V( 1, p ), 1 ) 872 CALL DAXPY( MVL, 873 $ CS*SN*APOAQ, 874 $ V( 1, p ), 1, 875 $ V( 1, q ), 1 ) 876 END IF 877 ELSE 878 CALL DAXPY( M, T*APOAQ, 879 $ A( 1, p ), 1, 880 $ A( 1, q ), 1 ) 881 CALL DAXPY( M, 882 $ -CS*SN*AQOAP, 883 $ A( 1, q ), 1, 884 $ A( 1, p ), 1 ) 885 D( p ) = D( p ) / CS 886 D( q ) = D( q )*CS 887 IF( RSVEC ) THEN 888 CALL DAXPY( MVL, 889 $ T*APOAQ, V( 1, p ), 890 $ 1, V( 1, q ), 1 ) 891 CALL DAXPY( MVL, 892 $ -CS*SN*AQOAP, 893 $ V( 1, q ), 1, 894 $ V( 1, p ), 1 ) 895 END IF 896 END IF 897 END IF 898 END IF 899 END IF 900* 901 ELSE 902 IF( AAPP.GT.AAQQ ) THEN 903 CALL DCOPY( M, A( 1, p ), 1, WORK, 904 $ 1 ) 905 CALL DLASCL( 'G', 0, 0, AAPP, ONE, 906 $ M, 1, WORK, LDA, IERR ) 907 CALL DLASCL( 'G', 0, 0, AAQQ, ONE, 908 $ M, 1, A( 1, q ), LDA, 909 $ IERR ) 910 TEMP1 = -AAPQ*D( p ) / D( q ) 911 CALL DAXPY( M, TEMP1, WORK, 1, 912 $ A( 1, q ), 1 ) 913 CALL DLASCL( 'G', 0, 0, ONE, AAQQ, 914 $ M, 1, A( 1, q ), LDA, 915 $ IERR ) 916 SVA( q ) = AAQQ*DSQRT( MAX( ZERO, 917 $ ONE-AAPQ*AAPQ ) ) 918 MXSINJ = MAX( MXSINJ, SFMIN ) 919 ELSE 920 CALL DCOPY( M, A( 1, q ), 1, WORK, 921 $ 1 ) 922 CALL DLASCL( 'G', 0, 0, AAQQ, ONE, 923 $ M, 1, WORK, LDA, IERR ) 924 CALL DLASCL( 'G', 0, 0, AAPP, ONE, 925 $ M, 1, A( 1, p ), LDA, 926 $ IERR ) 927 TEMP1 = -AAPQ*D( q ) / D( p ) 928 CALL DAXPY( M, TEMP1, WORK, 1, 929 $ A( 1, p ), 1 ) 930 CALL DLASCL( 'G', 0, 0, ONE, AAPP, 931 $ M, 1, A( 1, p ), LDA, 932 $ IERR ) 933 SVA( p ) = AAPP*DSQRT( MAX( ZERO, 934 $ ONE-AAPQ*AAPQ ) ) 935 MXSINJ = MAX( MXSINJ, SFMIN ) 936 END IF 937 END IF 938* END IF ROTOK THEN ... ELSE 939* 940* In the case of cancellation in updating SVA(q) 941* .. recompute SVA(q) 942 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 943 $ THEN 944 IF( ( AAQQ.LT.ROOTBIG ) .AND. 945 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN 946 SVA( q ) = DNRM2( M, A( 1, q ), 1 )* 947 $ D( q ) 948 ELSE 949 T = ZERO 950 AAQQ = ONE 951 CALL DLASSQ( M, A( 1, q ), 1, T, 952 $ AAQQ ) 953 SVA( q ) = T*DSQRT( AAQQ )*D( q ) 954 END IF 955 END IF 956 IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN 957 IF( ( AAPP.LT.ROOTBIG ) .AND. 958 $ ( AAPP.GT.ROOTSFMIN ) ) THEN 959 AAPP = DNRM2( M, A( 1, p ), 1 )* 960 $ D( p ) 961 ELSE 962 T = ZERO 963 AAPP = ONE 964 CALL DLASSQ( M, A( 1, p ), 1, T, 965 $ AAPP ) 966 AAPP = T*DSQRT( AAPP )*D( p ) 967 END IF 968 SVA( p ) = AAPP 969 END IF 970* end of OK rotation 971 ELSE 972 NOTROT = NOTROT + 1 973 PSKIPPED = PSKIPPED + 1 974 IJBLSK = IJBLSK + 1 975 END IF 976 ELSE 977 NOTROT = NOTROT + 1 978 PSKIPPED = PSKIPPED + 1 979 IJBLSK = IJBLSK + 1 980 END IF 981* 982 IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) 983 $ THEN 984 SVA( p ) = AAPP 985 NOTROT = 0 986 GO TO 2011 987 END IF 988 IF( ( i.LE.SWBAND ) .AND. 989 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN 990 AAPP = -AAPP 991 NOTROT = 0 992 GO TO 2203 993 END IF 994* 995 2200 CONTINUE 996* end of the q-loop 997 2203 CONTINUE 998* 999 SVA( p ) = AAPP 1000* 1001 ELSE 1002 IF( AAPP.EQ.ZERO )NOTROT = NOTROT + 1003 $ MIN( jgl+KBL-1, N ) - jgl + 1 1004 IF( AAPP.LT.ZERO )NOTROT = 0 1005 END IF 1006 1007 2100 CONTINUE 1008* end of the p-loop 1009 2010 CONTINUE 1010* end of the jbc-loop 1011 2011 CONTINUE 1012*2011 bailed out of the jbc-loop 1013 DO 2012 p = igl, MIN( igl+KBL-1, N ) 1014 SVA( p ) = DABS( SVA( p ) ) 1015 2012 CONTINUE 1016* 1017 2000 CONTINUE 1018*2000 :: end of the ibr-loop 1019* 1020* .. update SVA(N) 1021 IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) 1022 $ THEN 1023 SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N ) 1024 ELSE 1025 T = ZERO 1026 AAPP = ONE 1027 CALL DLASSQ( M, A( 1, N ), 1, T, AAPP ) 1028 SVA( N ) = T*DSQRT( AAPP )*D( N ) 1029 END IF 1030* 1031* Additional steering devices 1032* 1033 IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. 1034 $ ( ISWROT.LE.N ) ) )SWBAND = i 1035* 1036 IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND. 1037 $ ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN 1038 GO TO 1994 1039 END IF 1040* 1041 IF( NOTROT.GE.EMPTSW )GO TO 1994 1042 1043 1993 CONTINUE 1044* end i=1:NSWEEP loop 1045* #:) Reaching this point means that the procedure has completed the given 1046* number of iterations. 1047 INFO = NSWEEP - 1 1048 GO TO 1995 1049 1994 CONTINUE 1050* #:) Reaching this point means that during the i-th sweep all pivots were 1051* below the given tolerance, causing early exit. 1052* 1053 INFO = 0 1054* #:) INFO = 0 confirms successful iterations. 1055 1995 CONTINUE 1056* 1057* Sort the vector D. 1058 DO 5991 p = 1, N - 1 1059 q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1 1060 IF( p.NE.q ) THEN 1061 TEMP1 = SVA( p ) 1062 SVA( p ) = SVA( q ) 1063 SVA( q ) = TEMP1 1064 TEMP1 = D( p ) 1065 D( p ) = D( q ) 1066 D( q ) = TEMP1 1067 CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 1068 IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 ) 1069 END IF 1070 5991 CONTINUE 1071* 1072 RETURN 1073* .. 1074* .. END OF DGSVJ0 1075* .. 1076 END 1077