![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, 2: + EPS, 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: * .. 20: * .. Scalar Arguments .. 21: DOUBLE PRECISION EPS, SFMIN, TOL 22: INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP 23: CHARACTER*1 JOBV 24: * .. 25: * .. Array Arguments .. 26: DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ), 27: + WORK( LWORK ) 28: * .. 29: * 30: * Purpose 31: * ======= 32: * 33: * DGSVJ1 is called from SGESVJ as a pre-processor and that is its main 34: * purpose. It applies Jacobi rotations in the same way as SGESVJ does, but 35: * it targets only particular pivots and it does not check convergence 36: * (stopping criterion). Few tunning parameters (marked by [TP]) are 37: * available for the implementer. 38: * 39: * Further Details 40: * ~~~~~~~~~~~~~~~ 41: * DGSVJ1 applies few sweeps of Jacobi rotations in the column space of 42: * the input M-by-N matrix A. The pivot pairs are taken from the (1,2) 43: * off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The 44: * block-entries (tiles) of the (1,2) off-diagonal block are marked by the 45: * [x]'s in the following scheme: 46: * 47: * | * * * [x] [x] [x]| 48: * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. 49: * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. 50: * |[x] [x] [x] * * * | 51: * |[x] [x] [x] * * * | 52: * |[x] [x] [x] * * * | 53: * 54: * In terms of the columns of A, the first N1 columns are rotated 'against' 55: * the remaining N-N1 columns, trying to increase the angle between the 56: * corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is 57: * tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter. 58: * The number of sweeps is given in NSWEEP and the orthogonality threshold 59: * is given in TOL. 60: * 61: * Contributors 62: * ~~~~~~~~~~~~ 63: * Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) 64: * 65: * Arguments 66: * ========= 67: * 68: * JOBV (input) CHARACTER*1 69: * Specifies whether the output from this procedure is used 70: * to compute the matrix V: 71: * = 'V': the product of the Jacobi rotations is accumulated 72: * by postmulyiplying the N-by-N array V. 73: * (See the description of V.) 74: * = 'A': the product of the Jacobi rotations is accumulated 75: * by postmulyiplying the MV-by-N array V. 76: * (See the descriptions of MV and V.) 77: * = 'N': the Jacobi rotations are not accumulated. 78: * 79: * M (input) INTEGER 80: * The number of rows of the input matrix A. M >= 0. 81: * 82: * N (input) INTEGER 83: * The number of columns of the input matrix A. 84: * M >= N >= 0. 85: * 86: * N1 (input) INTEGER 87: * N1 specifies the 2 x 2 block partition, the first N1 columns are 88: * rotated 'against' the remaining N-N1 columns of A. 89: * 90: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 91: * On entry, M-by-N matrix A, such that A*diag(D) represents 92: * the input matrix. 93: * On exit, 94: * A_onexit * D_onexit represents the input matrix A*diag(D) 95: * post-multiplied by a sequence of Jacobi rotations, where the 96: * rotation threshold and the total number of sweeps are given in 97: * TOL and NSWEEP, respectively. 98: * (See the descriptions of N1, D, TOL and NSWEEP.) 99: * 100: * LDA (input) INTEGER 101: * The leading dimension of the array A. LDA >= max(1,M). 102: * 103: * D (input/workspace/output) DOUBLE PRECISION array, dimension (N) 104: * The array D accumulates the scaling factors from the fast scaled 105: * Jacobi rotations. 106: * On entry, A*diag(D) represents the input matrix. 107: * On exit, A_onexit*diag(D_onexit) represents the input matrix 108: * post-multiplied by a sequence of Jacobi rotations, where the 109: * rotation threshold and the total number of sweeps are given in 110: * TOL and NSWEEP, respectively. 111: * (See the descriptions of N1, A, TOL and NSWEEP.) 112: * 113: * SVA (input/workspace/output) DOUBLE PRECISION array, dimension (N) 114: * On entry, SVA contains the Euclidean norms of the columns of 115: * the matrix A*diag(D). 116: * On exit, SVA contains the Euclidean norms of the columns of 117: * the matrix onexit*diag(D_onexit). 118: * 119: * MV (input) 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: * 124: * V (input/output) DOUBLE PRECISION array, dimension (LDV,N) 125: * If JOBV .EQ. 'V' then N rows of V are post-multipled by a 126: * sequence of Jacobi rotations. 127: * If JOBV .EQ. 'A' then MV rows of V are post-multipled by a 128: * sequence of Jacobi rotations. 129: * If JOBV = 'N', then V is not referenced. 130: * 131: * LDV (input) INTEGER 132: * The leading dimension of the array V, LDV >= 1. 133: * If JOBV = 'V', LDV .GE. N. 134: * If JOBV = 'A', LDV .GE. MV. 135: * 136: * EPS (input) DOUBLE PRECISION 137: * EPS = DLAMCH('Epsilon') 138: * 139: * SFMIN (input) DOUBLE PRECISION 140: * SFMIN = DLAMCH('Safe Minimum') 141: * 142: * TOL (input) DOUBLE PRECISION 143: * TOL is the threshold for Jacobi rotations. For a pair 144: * A(:,p), A(:,q) of pivot columns, the Jacobi rotation is 145: * applied only if DABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL. 146: * 147: * NSWEEP (input) INTEGER 148: * NSWEEP is the number of sweeps of Jacobi rotations to be 149: * performed. 150: * 151: * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 152: * 153: * LWORK (input) INTEGER 154: * LWORK is the dimension of WORK. LWORK .GE. M. 155: * 156: * INFO (output) INTEGER 157: * = 0 : successful exit. 158: * < 0 : if INFO = -i, then the i-th argument had an illegal value 159: * 160: * ===================================================================== 161: * 162: * .. Local Parameters .. 163: DOUBLE PRECISION ZERO, HALF, ONE, TWO 164: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, 165: + TWO = 2.0D0 ) 166: * .. 167: * .. Local Scalars .. 168: DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, 169: + BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG, 170: + ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T, 171: + TEMP1, THETA, THSIGN 172: INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK, 173: + ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr, 174: + p, PSKIPPED, q, ROWSKIP, SWBAND 175: LOGICAL APPLV, ROTOK, RSVEC 176: * .. 177: * .. Local Arrays .. 178: DOUBLE PRECISION FASTR( 5 ) 179: * .. 180: * .. Intrinsic Functions .. 181: INTRINSIC DABS, DMAX1, DBLE, MIN0, DSIGN, DSQRT 182: * .. 183: * .. External Functions .. 184: DOUBLE PRECISION DDOT, DNRM2 185: INTEGER IDAMAX 186: LOGICAL LSAME 187: EXTERNAL IDAMAX, LSAME, DDOT, DNRM2 188: * .. 189: * .. External Subroutines .. 190: EXTERNAL DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP 191: * .. 192: * .. Executable Statements .. 193: * 194: * Test the input parameters. 195: * 196: APPLV = LSAME( JOBV, 'A' ) 197: RSVEC = LSAME( JOBV, 'V' ) 198: IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN 199: INFO = -1 200: ELSE IF( M.LT.0 ) THEN 201: INFO = -2 202: ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN 203: INFO = -3 204: ELSE IF( N1.LT.0 ) THEN 205: INFO = -4 206: ELSE IF( LDA.LT.M ) THEN 207: INFO = -6 208: ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN 209: INFO = -9 210: ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 211: & ( APPLV.AND.( LDV.LT.MV ) ) ) THEN 212: INFO = -11 213: ELSE IF( TOL.LE.EPS ) THEN 214: INFO = -14 215: ELSE IF( NSWEEP.LT.0 ) THEN 216: INFO = -15 217: ELSE IF( LWORK.LT.M ) THEN 218: INFO = -17 219: ELSE 220: INFO = 0 221: END IF 222: * 223: * #:( 224: IF( INFO.NE.0 ) THEN 225: CALL XERBLA( 'DGSVJ1', -INFO ) 226: RETURN 227: END IF 228: * 229: IF( RSVEC ) THEN 230: MVL = N 231: ELSE IF( APPLV ) THEN 232: MVL = MV 233: END IF 234: RSVEC = RSVEC .OR. APPLV 235: 236: ROOTEPS = DSQRT( EPS ) 237: ROOTSFMIN = DSQRT( SFMIN ) 238: SMALL = SFMIN / EPS 239: BIG = ONE / SFMIN 240: ROOTBIG = ONE / ROOTSFMIN 241: LARGE = BIG / DSQRT( DBLE( M*N ) ) 242: BIGTHETA = ONE / ROOTEPS 243: ROOTTOL = DSQRT( TOL ) 244: * 245: * .. Initialize the right singular vector matrix .. 246: * 247: * RSVEC = LSAME( JOBV, 'Y' ) 248: * 249: EMPTSW = N1*( N-N1 ) 250: NOTROT = 0 251: FASTR( 1 ) = ZERO 252: * 253: * .. Row-cyclic pivot strategy with de Rijk's pivoting .. 254: * 255: KBL = MIN0( 8, N ) 256: NBLR = N1 / KBL 257: IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1 258: 259: * .. the tiling is nblr-by-nblc [tiles] 260: 261: NBLC = ( N-N1 ) / KBL 262: IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1 263: BLSKIP = ( KBL**2 ) + 1 264: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL. 265: 266: ROWSKIP = MIN0( 5, KBL ) 267: *[TP] ROWSKIP is a tuning parameter. 268: SWBAND = 0 269: *[TP] SWBAND is a tuning parameter. It is meaningful and effective 270: * if SGESVJ is used as a computational routine in the preconditioned 271: * Jacobi SVD algorithm SGESVJ. 272: * 273: * 274: * | * * * [x] [x] [x]| 275: * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. 276: * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. 277: * |[x] [x] [x] * * * | 278: * |[x] [x] [x] * * * | 279: * |[x] [x] [x] * * * | 280: * 281: * 282: DO 1993 i = 1, NSWEEP 283: * .. go go go ... 284: * 285: MXAAPQ = ZERO 286: MXSINJ = ZERO 287: ISWROT = 0 288: * 289: NOTROT = 0 290: PSKIPPED = 0 291: * 292: DO 2000 ibr = 1, NBLR 293: 294: igl = ( ibr-1 )*KBL + 1 295: * 296: * 297: *........................................................ 298: * ... go to the off diagonal blocks 299: 300: igl = ( ibr-1 )*KBL + 1 301: 302: DO 2010 jbc = 1, NBLC 303: 304: jgl = N1 + ( jbc-1 )*KBL + 1 305: 306: * doing the block at ( ibr, jbc ) 307: 308: IJBLSK = 0 309: DO 2100 p = igl, MIN0( igl+KBL-1, N1 ) 310: 311: AAPP = SVA( p ) 312: 313: IF( AAPP.GT.ZERO ) THEN 314: 315: PSKIPPED = 0 316: 317: DO 2200 q = jgl, MIN0( jgl+KBL-1, N ) 318: * 319: AAQQ = SVA( q ) 320: 321: IF( AAQQ.GT.ZERO ) THEN 322: AAPP0 = AAPP 323: * 324: * .. M x 2 Jacobi SVD .. 325: * 326: * .. Safe Gram matrix computation .. 327: * 328: IF( AAQQ.GE.ONE ) THEN 329: IF( AAPP.GE.AAQQ ) THEN 330: ROTOK = ( SMALL*AAPP ).LE.AAQQ 331: ELSE 332: ROTOK = ( SMALL*AAQQ ).LE.AAPP 333: END IF 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: IF( AAPP.GE.AAQQ ) THEN 347: ROTOK = AAPP.LE.( AAQQ / SMALL ) 348: ELSE 349: ROTOK = AAQQ.LE.( AAPP / SMALL ) 350: END IF 351: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN 352: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, 353: + q ), 1 )*D( p )*D( q ) / AAQQ ) 354: + / AAPP 355: ELSE 356: CALL DCOPY( M, A( 1, q ), 1, WORK, 1 ) 357: CALL DLASCL( 'G', 0, 0, AAQQ, D( q ), 358: + M, 1, WORK, LDA, IERR ) 359: AAPQ = DDOT( M, WORK, 1, A( 1, p ), 360: + 1 )*D( p ) / AAPP 361: END IF 362: END IF 363: 364: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) ) 365: 366: * TO rotate or NOT to rotate, THAT is the question ... 367: * 368: IF( DABS( AAPQ ).GT.TOL ) THEN 369: NOTROT = 0 370: * ROTATED = ROTATED + 1 371: PSKIPPED = 0 372: ISWROT = ISWROT + 1 373: * 374: IF( ROTOK ) THEN 375: * 376: AQOAP = AAQQ / AAPP 377: APOAQ = AAPP / AAQQ 378: THETA = -HALF*DABS( AQOAP-APOAQ ) / 379: + AAPQ 380: IF( AAQQ.GT.AAPP0 )THETA = -THETA 381: 382: IF( DABS( THETA ).GT.BIGTHETA ) THEN 383: T = HALF / THETA 384: FASTR( 3 ) = T*D( p ) / D( q ) 385: FASTR( 4 ) = -T*D( q ) / D( p ) 386: CALL DROTM( M, A( 1, p ), 1, 387: + A( 1, q ), 1, FASTR ) 388: IF( RSVEC )CALL DROTM( MVL, 389: + V( 1, p ), 1, 390: + V( 1, q ), 1, 391: + FASTR ) 392: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, 393: + ONE+T*APOAQ*AAPQ ) ) 394: AAPP = AAPP*DSQRT( DMAX1( ZERO, 395: + ONE-T*AQOAP*AAPQ ) ) 396: MXSINJ = DMAX1( MXSINJ, DABS( T ) ) 397: ELSE 398: * 399: * .. choose correct signum for THETA and rotate 400: * 401: THSIGN = -DSIGN( ONE, AAPQ ) 402: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN 403: T = ONE / ( THETA+THSIGN* 404: + DSQRT( ONE+THETA*THETA ) ) 405: CS = DSQRT( ONE / ( ONE+T*T ) ) 406: SN = T*CS 407: MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) 408: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, 409: + ONE+T*APOAQ*AAPQ ) ) 410: AAPP = AAPP*DSQRT( DMAX1( ZERO, 411: + ONE-T*AQOAP*AAPQ ) ) 412: 413: APOAQ = D( p ) / D( q ) 414: AQOAP = D( q ) / D( p ) 415: IF( D( p ).GE.ONE ) THEN 416: * 417: IF( D( q ).GE.ONE ) THEN 418: FASTR( 3 ) = T*APOAQ 419: FASTR( 4 ) = -T*AQOAP 420: D( p ) = D( p )*CS 421: D( q ) = D( q )*CS 422: CALL DROTM( M, A( 1, p ), 1, 423: + A( 1, q ), 1, 424: + FASTR ) 425: IF( RSVEC )CALL DROTM( MVL, 426: + V( 1, p ), 1, V( 1, q ), 427: + 1, FASTR ) 428: ELSE 429: CALL DAXPY( M, -T*AQOAP, 430: + A( 1, q ), 1, 431: + A( 1, p ), 1 ) 432: CALL DAXPY( M, CS*SN*APOAQ, 433: + A( 1, p ), 1, 434: + A( 1, q ), 1 ) 435: IF( RSVEC ) THEN 436: CALL DAXPY( MVL, -T*AQOAP, 437: + V( 1, q ), 1, 438: + V( 1, p ), 1 ) 439: CALL DAXPY( MVL, 440: + CS*SN*APOAQ, 441: + V( 1, p ), 1, 442: + V( 1, q ), 1 ) 443: END IF 444: D( p ) = D( p )*CS 445: D( q ) = D( q ) / CS 446: END IF 447: ELSE 448: IF( D( q ).GE.ONE ) THEN 449: CALL DAXPY( M, T*APOAQ, 450: + A( 1, p ), 1, 451: + A( 1, q ), 1 ) 452: CALL DAXPY( M, -CS*SN*AQOAP, 453: + A( 1, q ), 1, 454: + A( 1, p ), 1 ) 455: IF( RSVEC ) THEN 456: CALL DAXPY( MVL, T*APOAQ, 457: + V( 1, p ), 1, 458: + V( 1, q ), 1 ) 459: CALL DAXPY( MVL, 460: + -CS*SN*AQOAP, 461: + V( 1, q ), 1, 462: + V( 1, p ), 1 ) 463: END IF 464: D( p ) = D( p ) / CS 465: D( q ) = D( q )*CS 466: ELSE 467: IF( D( p ).GE.D( q ) ) THEN 468: CALL DAXPY( M, -T*AQOAP, 469: + A( 1, q ), 1, 470: + A( 1, p ), 1 ) 471: CALL DAXPY( M, CS*SN*APOAQ, 472: + A( 1, p ), 1, 473: + A( 1, q ), 1 ) 474: D( p ) = D( p )*CS 475: D( q ) = D( q ) / CS 476: IF( RSVEC ) THEN 477: CALL DAXPY( MVL, 478: + -T*AQOAP, 479: + V( 1, q ), 1, 480: + V( 1, p ), 1 ) 481: CALL DAXPY( MVL, 482: + CS*SN*APOAQ, 483: + V( 1, p ), 1, 484: + V( 1, q ), 1 ) 485: END IF 486: ELSE 487: CALL DAXPY( M, T*APOAQ, 488: + A( 1, p ), 1, 489: + A( 1, q ), 1 ) 490: CALL DAXPY( M, 491: + -CS*SN*AQOAP, 492: + A( 1, q ), 1, 493: + A( 1, p ), 1 ) 494: D( p ) = D( p ) / CS 495: D( q ) = D( q )*CS 496: IF( RSVEC ) THEN 497: CALL DAXPY( MVL, 498: + T*APOAQ, V( 1, p ), 499: + 1, V( 1, q ), 1 ) 500: CALL DAXPY( MVL, 501: + -CS*SN*AQOAP, 502: + V( 1, q ), 1, 503: + V( 1, p ), 1 ) 504: END IF 505: END IF 506: END IF 507: END IF 508: END IF 509: 510: ELSE 511: IF( AAPP.GT.AAQQ ) THEN 512: CALL DCOPY( M, A( 1, p ), 1, WORK, 513: + 1 ) 514: CALL DLASCL( 'G', 0, 0, AAPP, ONE, 515: + M, 1, WORK, LDA, IERR ) 516: CALL DLASCL( 'G', 0, 0, AAQQ, ONE, 517: + M, 1, A( 1, q ), LDA, 518: + IERR ) 519: TEMP1 = -AAPQ*D( p ) / D( q ) 520: CALL DAXPY( M, TEMP1, WORK, 1, 521: + A( 1, q ), 1 ) 522: CALL DLASCL( 'G', 0, 0, ONE, AAQQ, 523: + M, 1, A( 1, q ), LDA, 524: + IERR ) 525: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, 526: + ONE-AAPQ*AAPQ ) ) 527: MXSINJ = DMAX1( MXSINJ, SFMIN ) 528: ELSE 529: CALL DCOPY( M, A( 1, q ), 1, WORK, 530: + 1 ) 531: CALL DLASCL( 'G', 0, 0, AAQQ, ONE, 532: + M, 1, WORK, LDA, IERR ) 533: CALL DLASCL( 'G', 0, 0, AAPP, ONE, 534: + M, 1, A( 1, p ), LDA, 535: + IERR ) 536: TEMP1 = -AAPQ*D( q ) / D( p ) 537: CALL DAXPY( M, TEMP1, WORK, 1, 538: + A( 1, p ), 1 ) 539: CALL DLASCL( 'G', 0, 0, ONE, AAPP, 540: + M, 1, A( 1, p ), LDA, 541: + IERR ) 542: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO, 543: + ONE-AAPQ*AAPQ ) ) 544: MXSINJ = DMAX1( MXSINJ, SFMIN ) 545: END IF 546: END IF 547: * END IF ROTOK THEN ... ELSE 548: * 549: * In the case of cancellation in updating SVA(q) 550: * .. recompute SVA(q) 551: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 552: + THEN 553: IF( ( AAQQ.LT.ROOTBIG ) .AND. 554: + ( AAQQ.GT.ROOTSFMIN ) ) THEN 555: SVA( q ) = DNRM2( M, A( 1, q ), 1 )* 556: + D( q ) 557: ELSE 558: T = ZERO 559: AAQQ = ONE 560: CALL DLASSQ( M, A( 1, q ), 1, T, 561: + AAQQ ) 562: SVA( q ) = T*DSQRT( AAQQ )*D( q ) 563: END IF 564: END IF 565: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN 566: IF( ( AAPP.LT.ROOTBIG ) .AND. 567: + ( AAPP.GT.ROOTSFMIN ) ) THEN 568: AAPP = DNRM2( M, A( 1, p ), 1 )* 569: + D( p ) 570: ELSE 571: T = ZERO 572: AAPP = ONE 573: CALL DLASSQ( M, A( 1, p ), 1, T, 574: + AAPP ) 575: AAPP = T*DSQRT( AAPP )*D( p ) 576: END IF 577: SVA( p ) = AAPP 578: END IF 579: * end of OK rotation 580: ELSE 581: NOTROT = NOTROT + 1 582: * SKIPPED = SKIPPED + 1 583: PSKIPPED = PSKIPPED + 1 584: IJBLSK = IJBLSK + 1 585: END IF 586: ELSE 587: NOTROT = NOTROT + 1 588: PSKIPPED = PSKIPPED + 1 589: IJBLSK = IJBLSK + 1 590: END IF 591: 592: * IF ( NOTROT .GE. EMPTSW ) GO TO 2011 593: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) 594: + THEN 595: SVA( p ) = AAPP 596: NOTROT = 0 597: GO TO 2011 598: END IF 599: IF( ( i.LE.SWBAND ) .AND. 600: + ( PSKIPPED.GT.ROWSKIP ) ) THEN 601: AAPP = -AAPP 602: NOTROT = 0 603: GO TO 2203 604: END IF 605: 606: * 607: 2200 CONTINUE 608: * end of the q-loop 609: 2203 CONTINUE 610: 611: SVA( p ) = AAPP 612: * 613: ELSE 614: IF( AAPP.EQ.ZERO )NOTROT = NOTROT + 615: + MIN0( jgl+KBL-1, N ) - jgl + 1 616: IF( AAPP.LT.ZERO )NOTROT = 0 617: *** IF ( NOTROT .GE. EMPTSW ) GO TO 2011 618: END IF 619: 620: 2100 CONTINUE 621: * end of the p-loop 622: 2010 CONTINUE 623: * end of the jbc-loop 624: 2011 CONTINUE 625: *2011 bailed out of the jbc-loop 626: DO 2012 p = igl, MIN0( igl+KBL-1, N ) 627: SVA( p ) = DABS( SVA( p ) ) 628: 2012 CONTINUE 629: *** IF ( NOTROT .GE. EMPTSW ) GO TO 1994 630: 2000 CONTINUE 631: *2000 :: end of the ibr-loop 632: * 633: * .. update SVA(N) 634: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) 635: + THEN 636: SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N ) 637: ELSE 638: T = ZERO 639: AAPP = ONE 640: CALL DLASSQ( M, A( 1, N ), 1, T, AAPP ) 641: SVA( N ) = T*DSQRT( AAPP )*D( N ) 642: END IF 643: * 644: * Additional steering devices 645: * 646: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. 647: + ( ISWROT.LE.N ) ) )SWBAND = i 648: 649: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND. 650: + ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN 651: GO TO 1994 652: END IF 653: 654: * 655: IF( NOTROT.GE.EMPTSW )GO TO 1994 656: 657: 1993 CONTINUE 658: * end i=1:NSWEEP loop 659: * #:) Reaching this point means that the procedure has completed the given 660: * number of sweeps. 661: INFO = NSWEEP - 1 662: GO TO 1995 663: 1994 CONTINUE 664: * #:) Reaching this point means that during the i-th sweep all pivots were 665: * below the given threshold, causing early exit. 666: 667: INFO = 0 668: * #:) INFO = 0 confirms successful iterations. 669: 1995 CONTINUE 670: * 671: * Sort the vector D 672: * 673: DO 5991 p = 1, N - 1 674: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1 675: IF( p.NE.q ) THEN 676: TEMP1 = SVA( p ) 677: SVA( p ) = SVA( q ) 678: SVA( q ) = TEMP1 679: TEMP1 = D( p ) 680: D( p ) = D( q ) 681: D( q ) = TEMP1 682: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 683: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 ) 684: END IF 685: 5991 CONTINUE 686: * 687: RETURN 688: * .. 689: * .. END OF DGSVJ1 690: * .. 691: END