Annotation of rpl/lapack/lapack/dgsvj1.f, revision 1.1
1.1 ! bertrand 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.2.2) --
! 5: *
! 6: * -- Contributed by Zlatko Drmac of the University of Zagreb and --
! 7: * -- Kresimir Veselic of the Fernuniversitaet Hagen --
! 8: * -- June 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( MV.LT.0 ) THEN
! 209: INFO = -9
! 210: ELSE IF( LDV.LT.M ) THEN
! 211: INFO = -11
! 212: ELSE IF( TOL.LE.EPS ) THEN
! 213: INFO = -14
! 214: ELSE IF( NSWEEP.LT.0 ) THEN
! 215: INFO = -15
! 216: ELSE IF( LWORK.LT.M ) THEN
! 217: INFO = -17
! 218: ELSE
! 219: INFO = 0
! 220: END IF
! 221: *
! 222: * #:(
! 223: IF( INFO.NE.0 ) THEN
! 224: CALL XERBLA( 'DGSVJ1', -INFO )
! 225: RETURN
! 226: END IF
! 227: *
! 228: IF( RSVEC ) THEN
! 229: MVL = N
! 230: ELSE IF( APPLV ) THEN
! 231: MVL = MV
! 232: END IF
! 233: RSVEC = RSVEC .OR. APPLV
! 234:
! 235: ROOTEPS = DSQRT( EPS )
! 236: ROOTSFMIN = DSQRT( SFMIN )
! 237: SMALL = SFMIN / EPS
! 238: BIG = ONE / SFMIN
! 239: ROOTBIG = ONE / ROOTSFMIN
! 240: LARGE = BIG / DSQRT( DBLE( M*N ) )
! 241: BIGTHETA = ONE / ROOTEPS
! 242: ROOTTOL = DSQRT( TOL )
! 243: *
! 244: * .. Initialize the right singular vector matrix ..
! 245: *
! 246: * RSVEC = LSAME( JOBV, 'Y' )
! 247: *
! 248: EMPTSW = N1*( N-N1 )
! 249: NOTROT = 0
! 250: FASTR( 1 ) = ZERO
! 251: *
! 252: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
! 253: *
! 254: KBL = MIN0( 8, N )
! 255: NBLR = N1 / KBL
! 256: IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
! 257:
! 258: * .. the tiling is nblr-by-nblc [tiles]
! 259:
! 260: NBLC = ( N-N1 ) / KBL
! 261: IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
! 262: BLSKIP = ( KBL**2 ) + 1
! 263: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
! 264:
! 265: ROWSKIP = MIN0( 5, KBL )
! 266: *[TP] ROWSKIP is a tuning parameter.
! 267: SWBAND = 0
! 268: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
! 269: * if SGESVJ is used as a computational routine in the preconditioned
! 270: * Jacobi SVD algorithm SGESVJ.
! 271: *
! 272: *
! 273: * | * * * [x] [x] [x]|
! 274: * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
! 275: * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
! 276: * |[x] [x] [x] * * * |
! 277: * |[x] [x] [x] * * * |
! 278: * |[x] [x] [x] * * * |
! 279: *
! 280: *
! 281: DO 1993 i = 1, NSWEEP
! 282: * .. go go go ...
! 283: *
! 284: MXAAPQ = ZERO
! 285: MXSINJ = ZERO
! 286: ISWROT = 0
! 287: *
! 288: NOTROT = 0
! 289: PSKIPPED = 0
! 290: *
! 291: DO 2000 ibr = 1, NBLR
! 292:
! 293: igl = ( ibr-1 )*KBL + 1
! 294: *
! 295: *
! 296: *........................................................
! 297: * ... go to the off diagonal blocks
! 298:
! 299: igl = ( ibr-1 )*KBL + 1
! 300:
! 301: DO 2010 jbc = 1, NBLC
! 302:
! 303: jgl = N1 + ( jbc-1 )*KBL + 1
! 304:
! 305: * doing the block at ( ibr, jbc )
! 306:
! 307: IJBLSK = 0
! 308: DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
! 309:
! 310: AAPP = SVA( p )
! 311:
! 312: IF( AAPP.GT.ZERO ) THEN
! 313:
! 314: PSKIPPED = 0
! 315:
! 316: DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
! 317: *
! 318: AAQQ = SVA( q )
! 319:
! 320: IF( AAQQ.GT.ZERO ) THEN
! 321: AAPP0 = AAPP
! 322: *
! 323: * .. M x 2 Jacobi SVD ..
! 324: *
! 325: * .. Safe Gram matrix computation ..
! 326: *
! 327: IF( AAQQ.GE.ONE ) THEN
! 328: IF( AAPP.GE.AAQQ ) THEN
! 329: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 330: ELSE
! 331: ROTOK = ( SMALL*AAQQ ).LE.AAPP
! 332: END IF
! 333: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
! 334: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 335: + q ), 1 )*D( p )*D( q ) / AAQQ )
! 336: + / AAPP
! 337: ELSE
! 338: CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
! 339: CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
! 340: + M, 1, WORK, LDA, IERR )
! 341: AAPQ = DDOT( M, WORK, 1, A( 1, q ),
! 342: + 1 )*D( q ) / AAQQ
! 343: END IF
! 344: ELSE
! 345: IF( AAPP.GE.AAQQ ) THEN
! 346: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 347: ELSE
! 348: ROTOK = AAQQ.LE.( AAPP / SMALL )
! 349: END IF
! 350: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 351: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 352: + q ), 1 )*D( p )*D( q ) / AAQQ )
! 353: + / AAPP
! 354: ELSE
! 355: CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
! 356: CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
! 357: + M, 1, WORK, LDA, IERR )
! 358: AAPQ = DDOT( M, WORK, 1, A( 1, p ),
! 359: + 1 )*D( p ) / AAPP
! 360: END IF
! 361: END IF
! 362:
! 363: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
! 364:
! 365: * TO rotate or NOT to rotate, THAT is the question ...
! 366: *
! 367: IF( DABS( AAPQ ).GT.TOL ) THEN
! 368: NOTROT = 0
! 369: * ROTATED = ROTATED + 1
! 370: PSKIPPED = 0
! 371: ISWROT = ISWROT + 1
! 372: *
! 373: IF( ROTOK ) THEN
! 374: *
! 375: AQOAP = AAQQ / AAPP
! 376: APOAQ = AAPP / AAQQ
! 377: THETA = -HALF*DABS( AQOAP-APOAQ ) /
! 378: + AAPQ
! 379: IF( AAQQ.GT.AAPP0 )THETA = -THETA
! 380:
! 381: IF( DABS( THETA ).GT.BIGTHETA ) THEN
! 382: T = HALF / THETA
! 383: FASTR( 3 ) = T*D( p ) / D( q )
! 384: FASTR( 4 ) = -T*D( q ) / D( p )
! 385: CALL DROTM( M, A( 1, p ), 1,
! 386: + A( 1, q ), 1, FASTR )
! 387: IF( RSVEC )CALL DROTM( MVL,
! 388: + V( 1, p ), 1,
! 389: + V( 1, q ), 1,
! 390: + FASTR )
! 391: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 392: + ONE+T*APOAQ*AAPQ ) )
! 393: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 394: + ONE-T*AQOAP*AAPQ ) )
! 395: MXSINJ = DMAX1( MXSINJ, DABS( T ) )
! 396: ELSE
! 397: *
! 398: * .. choose correct signum for THETA and rotate
! 399: *
! 400: THSIGN = -DSIGN( ONE, AAPQ )
! 401: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
! 402: T = ONE / ( THETA+THSIGN*
! 403: + DSQRT( ONE+THETA*THETA ) )
! 404: CS = DSQRT( ONE / ( ONE+T*T ) )
! 405: SN = T*CS
! 406: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
! 407: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 408: + ONE+T*APOAQ*AAPQ ) )
! 409: AAPP = AAPP*DSQRT( ONE-T*AQOAP*
! 410: + AAPQ )
! 411:
! 412: APOAQ = D( p ) / D( q )
! 413: AQOAP = D( q ) / D( p )
! 414: IF( D( p ).GE.ONE ) THEN
! 415: *
! 416: IF( D( q ).GE.ONE ) THEN
! 417: FASTR( 3 ) = T*APOAQ
! 418: FASTR( 4 ) = -T*AQOAP
! 419: D( p ) = D( p )*CS
! 420: D( q ) = D( q )*CS
! 421: CALL DROTM( M, A( 1, p ), 1,
! 422: + A( 1, q ), 1,
! 423: + FASTR )
! 424: IF( RSVEC )CALL DROTM( MVL,
! 425: + V( 1, p ), 1, V( 1, q ),
! 426: + 1, FASTR )
! 427: ELSE
! 428: CALL DAXPY( M, -T*AQOAP,
! 429: + A( 1, q ), 1,
! 430: + A( 1, p ), 1 )
! 431: CALL DAXPY( M, CS*SN*APOAQ,
! 432: + A( 1, p ), 1,
! 433: + A( 1, q ), 1 )
! 434: IF( RSVEC ) THEN
! 435: CALL DAXPY( MVL, -T*AQOAP,
! 436: + V( 1, q ), 1,
! 437: + V( 1, p ), 1 )
! 438: CALL DAXPY( MVL,
! 439: + CS*SN*APOAQ,
! 440: + V( 1, p ), 1,
! 441: + V( 1, q ), 1 )
! 442: END IF
! 443: D( p ) = D( p )*CS
! 444: D( q ) = D( q ) / CS
! 445: END IF
! 446: ELSE
! 447: IF( D( q ).GE.ONE ) THEN
! 448: CALL DAXPY( M, T*APOAQ,
! 449: + A( 1, p ), 1,
! 450: + A( 1, q ), 1 )
! 451: CALL DAXPY( M, -CS*SN*AQOAP,
! 452: + A( 1, q ), 1,
! 453: + A( 1, p ), 1 )
! 454: IF( RSVEC ) THEN
! 455: CALL DAXPY( MVL, T*APOAQ,
! 456: + V( 1, p ), 1,
! 457: + V( 1, q ), 1 )
! 458: CALL DAXPY( MVL,
! 459: + -CS*SN*AQOAP,
! 460: + V( 1, q ), 1,
! 461: + V( 1, p ), 1 )
! 462: END IF
! 463: D( p ) = D( p ) / CS
! 464: D( q ) = D( q )*CS
! 465: ELSE
! 466: IF( D( p ).GE.D( q ) ) THEN
! 467: CALL DAXPY( M, -T*AQOAP,
! 468: + A( 1, q ), 1,
! 469: + A( 1, p ), 1 )
! 470: CALL DAXPY( M, CS*SN*APOAQ,
! 471: + A( 1, p ), 1,
! 472: + A( 1, q ), 1 )
! 473: D( p ) = D( p )*CS
! 474: D( q ) = D( q ) / CS
! 475: IF( RSVEC ) THEN
! 476: CALL DAXPY( MVL,
! 477: + -T*AQOAP,
! 478: + V( 1, q ), 1,
! 479: + V( 1, p ), 1 )
! 480: CALL DAXPY( MVL,
! 481: + CS*SN*APOAQ,
! 482: + V( 1, p ), 1,
! 483: + V( 1, q ), 1 )
! 484: END IF
! 485: ELSE
! 486: CALL DAXPY( M, T*APOAQ,
! 487: + A( 1, p ), 1,
! 488: + A( 1, q ), 1 )
! 489: CALL DAXPY( M,
! 490: + -CS*SN*AQOAP,
! 491: + A( 1, q ), 1,
! 492: + A( 1, p ), 1 )
! 493: D( p ) = D( p ) / CS
! 494: D( q ) = D( q )*CS
! 495: IF( RSVEC ) THEN
! 496: CALL DAXPY( MVL,
! 497: + T*APOAQ, V( 1, p ),
! 498: + 1, V( 1, q ), 1 )
! 499: CALL DAXPY( MVL,
! 500: + -CS*SN*AQOAP,
! 501: + V( 1, q ), 1,
! 502: + V( 1, p ), 1 )
! 503: END IF
! 504: END IF
! 505: END IF
! 506: END IF
! 507: END IF
! 508:
! 509: ELSE
! 510: IF( AAPP.GT.AAQQ ) THEN
! 511: CALL DCOPY( M, A( 1, p ), 1, WORK,
! 512: + 1 )
! 513: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
! 514: + M, 1, WORK, LDA, IERR )
! 515: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
! 516: + M, 1, A( 1, q ), LDA,
! 517: + 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,
! 522: + M, 1, A( 1, q ), LDA,
! 523: + IERR )
! 524: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 525: + ONE-AAPQ*AAPQ ) )
! 526: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 527: ELSE
! 528: CALL DCOPY( M, A( 1, q ), 1, WORK,
! 529: + 1 )
! 530: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
! 531: + M, 1, WORK, LDA, IERR )
! 532: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
! 533: + M, 1, A( 1, p ), LDA,
! 534: + IERR )
! 535: TEMP1 = -AAPQ*D( q ) / D( p )
! 536: CALL DAXPY( M, TEMP1, WORK, 1,
! 537: + A( 1, p ), 1 )
! 538: CALL DLASCL( 'G', 0, 0, ONE, AAPP,
! 539: + M, 1, A( 1, p ), LDA,
! 540: + IERR )
! 541: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
! 542: + ONE-AAPQ*AAPQ ) )
! 543: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 544: END IF
! 545: END IF
! 546: * END IF ROTOK THEN ... ELSE
! 547: *
! 548: * In the case of cancellation in updating SVA(q)
! 549: * .. recompute SVA(q)
! 550: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 551: + THEN
! 552: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 553: + ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 554: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
! 555: + D( q )
! 556: ELSE
! 557: T = ZERO
! 558: AAQQ = ZERO
! 559: CALL DLASSQ( M, A( 1, q ), 1, T,
! 560: + AAQQ )
! 561: SVA( q ) = T*DSQRT( AAQQ )*D( q )
! 562: END IF
! 563: END IF
! 564: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
! 565: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 566: + ( AAPP.GT.ROOTSFMIN ) ) THEN
! 567: AAPP = DNRM2( M, A( 1, p ), 1 )*
! 568: + D( p )
! 569: ELSE
! 570: T = ZERO
! 571: AAPP = ZERO
! 572: CALL DLASSQ( M, A( 1, p ), 1, T,
! 573: + AAPP )
! 574: AAPP = T*DSQRT( AAPP )*D( p )
! 575: END IF
! 576: SVA( p ) = AAPP
! 577: END IF
! 578: * end of OK rotation
! 579: ELSE
! 580: NOTROT = NOTROT + 1
! 581: * SKIPPED = SKIPPED + 1
! 582: PSKIPPED = PSKIPPED + 1
! 583: IJBLSK = IJBLSK + 1
! 584: END IF
! 585: ELSE
! 586: NOTROT = NOTROT + 1
! 587: PSKIPPED = PSKIPPED + 1
! 588: IJBLSK = IJBLSK + 1
! 589: END IF
! 590:
! 591: * IF ( NOTROT .GE. EMPTSW ) GO TO 2011
! 592: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
! 593: + THEN
! 594: SVA( p ) = AAPP
! 595: NOTROT = 0
! 596: GO TO 2011
! 597: END IF
! 598: IF( ( i.LE.SWBAND ) .AND.
! 599: + ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 600: AAPP = -AAPP
! 601: NOTROT = 0
! 602: GO TO 2203
! 603: END IF
! 604:
! 605: *
! 606: 2200 CONTINUE
! 607: * end of the q-loop
! 608: 2203 CONTINUE
! 609:
! 610: SVA( p ) = AAPP
! 611: *
! 612: ELSE
! 613: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
! 614: + MIN0( jgl+KBL-1, N ) - jgl + 1
! 615: IF( AAPP.LT.ZERO )NOTROT = 0
! 616: *** IF ( NOTROT .GE. EMPTSW ) GO TO 2011
! 617: END IF
! 618:
! 619: 2100 CONTINUE
! 620: * end of the p-loop
! 621: 2010 CONTINUE
! 622: * end of the jbc-loop
! 623: 2011 CONTINUE
! 624: *2011 bailed out of the jbc-loop
! 625: DO 2012 p = igl, MIN0( igl+KBL-1, N )
! 626: SVA( p ) = DABS( SVA( p ) )
! 627: 2012 CONTINUE
! 628: *** IF ( NOTROT .GE. EMPTSW ) GO TO 1994
! 629: 2000 CONTINUE
! 630: *2000 :: end of the ibr-loop
! 631: *
! 632: * .. update SVA(N)
! 633: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
! 634: + THEN
! 635: SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
! 636: ELSE
! 637: T = ZERO
! 638: AAPP = ZERO
! 639: CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
! 640: SVA( N ) = T*DSQRT( AAPP )*D( N )
! 641: END IF
! 642: *
! 643: * Additional steering devices
! 644: *
! 645: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
! 646: + ( ISWROT.LE.N ) ) )SWBAND = i
! 647:
! 648: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
! 649: + ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
! 650: GO TO 1994
! 651: END IF
! 652:
! 653: *
! 654: IF( NOTROT.GE.EMPTSW )GO TO 1994
! 655:
! 656: 1993 CONTINUE
! 657: * end i=1:NSWEEP loop
! 658: * #:) Reaching this point means that the procedure has completed the given
! 659: * number of sweeps.
! 660: INFO = NSWEEP - 1
! 661: GO TO 1995
! 662: 1994 CONTINUE
! 663: * #:) Reaching this point means that during the i-th sweep all pivots were
! 664: * below the given threshold, causing early exit.
! 665:
! 666: INFO = 0
! 667: * #:) INFO = 0 confirms successful iterations.
! 668: 1995 CONTINUE
! 669: *
! 670: * Sort the vector D
! 671: *
! 672: DO 5991 p = 1, N - 1
! 673: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 674: IF( p.NE.q ) THEN
! 675: TEMP1 = SVA( p )
! 676: SVA( p ) = SVA( q )
! 677: SVA( q ) = TEMP1
! 678: TEMP1 = D( p )
! 679: D( p ) = D( q )
! 680: D( q ) = TEMP1
! 681: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 682: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
! 683: END IF
! 684: 5991 CONTINUE
! 685: *
! 686: RETURN
! 687: * ..
! 688: * .. END OF DGSVJ1
! 689: * ..
! 690: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgsvj1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack