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