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