Annotation of rpl/lapack/lapack/dgsvj0.f, revision 1.1
1.1 ! bertrand 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.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: * .. 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( MV.LT.0 ) THEN
! 187: INFO = -8
! 188: ELSE IF( LDV.LT.M ) THEN
! 189: INFO = -10
! 190: ELSE IF( TOL.LE.EPS ) THEN
! 191: INFO = -13
! 192: ELSE IF( NSWEEP.LT.0 ) THEN
! 193: INFO = -14
! 194: ELSE IF( LWORK.LT.M ) THEN
! 195: INFO = -16
! 196: ELSE
! 197: INFO = 0
! 198: END IF
! 199: *
! 200: * #:(
! 201: IF( INFO.NE.0 ) THEN
! 202: CALL XERBLA( 'DGSVJ0', -INFO )
! 203: RETURN
! 204: END IF
! 205: *
! 206: IF( RSVEC ) THEN
! 207: MVL = N
! 208: ELSE IF( APPLV ) THEN
! 209: MVL = MV
! 210: END IF
! 211: RSVEC = RSVEC .OR. APPLV
! 212:
! 213: ROOTEPS = DSQRT( EPS )
! 214: ROOTSFMIN = DSQRT( SFMIN )
! 215: SMALL = SFMIN / EPS
! 216: BIG = ONE / SFMIN
! 217: ROOTBIG = ONE / ROOTSFMIN
! 218: BIGTHETA = ONE / ROOTEPS
! 219: ROOTTOL = DSQRT( TOL )
! 220: *
! 221: *
! 222: * -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#-
! 223: *
! 224: EMPTSW = ( N*( N-1 ) ) / 2
! 225: NOTROT = 0
! 226: FASTR( 1 ) = ZERO
! 227: *
! 228: * -#- Row-cyclic pivot strategy with de Rijk's pivoting -#-
! 229: *
! 230:
! 231: SWBAND = 0
! 232: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
! 233: * if SGESVJ is used as a computational routine in the preconditioned
! 234: * Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
! 235: * ......
! 236:
! 237: KBL = MIN0( 8, N )
! 238: *[TP] KBL is a tuning parameter that defines the tile size in the
! 239: * tiling of the p-q loops of pivot pairs. In general, an optimal
! 240: * value of KBL depends on the matrix dimensions and on the
! 241: * parameters of the computer's memory.
! 242: *
! 243: NBL = N / KBL
! 244: IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
! 245:
! 246: BLSKIP = ( KBL**2 ) + 1
! 247: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
! 248:
! 249: ROWSKIP = MIN0( 5, KBL )
! 250: *[TP] ROWSKIP is a tuning parameter.
! 251:
! 252: LKAHEAD = 1
! 253: *[TP] LKAHEAD is a tuning parameter.
! 254: SWBAND = 0
! 255: PSKIPPED = 0
! 256: *
! 257: DO 1993 i = 1, NSWEEP
! 258: * .. go go go ...
! 259: *
! 260: MXAAPQ = ZERO
! 261: MXSINJ = ZERO
! 262: ISWROT = 0
! 263: *
! 264: NOTROT = 0
! 265: PSKIPPED = 0
! 266: *
! 267: DO 2000 ibr = 1, NBL
! 268:
! 269: igl = ( ibr-1 )*KBL + 1
! 270: *
! 271: DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
! 272: *
! 273: igl = igl + ir1*KBL
! 274: *
! 275: DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
! 276:
! 277: * .. de Rijk's pivoting
! 278: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 279: IF( p.NE.q ) THEN
! 280: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 281: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,
! 282: + V( 1, q ), 1 )
! 283: TEMP1 = SVA( p )
! 284: SVA( p ) = SVA( q )
! 285: SVA( q ) = TEMP1
! 286: TEMP1 = D( p )
! 287: D( p ) = D( q )
! 288: D( q ) = TEMP1
! 289: END IF
! 290: *
! 291: IF( ir1.EQ.0 ) THEN
! 292: *
! 293: * Column norms are periodically updated by explicit
! 294: * norm computation.
! 295: * Caveat:
! 296: * Some BLAS implementations compute DNRM2(M,A(1,p),1)
! 297: * as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may result in
! 298: * overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and
! 299: * undeflow for ||A(:,p)||_2 < DSQRT(underflow_threshold).
! 300: * Hence, DNRM2 cannot be trusted, not even in the case when
! 301: * the true norm is far from the under(over)flow boundaries.
! 302: * If properly implemented DNRM2 is available, the IF-THEN-ELSE
! 303: * below should read "AAPP = DNRM2( M, A(1,p), 1 ) * D(p)".
! 304: *
! 305: IF( ( SVA( p ).LT.ROOTBIG ) .AND.
! 306: + ( SVA( p ).GT.ROOTSFMIN ) ) THEN
! 307: SVA( p ) = DNRM2( M, A( 1, p ), 1 )*D( p )
! 308: ELSE
! 309: TEMP1 = ZERO
! 310: AAPP = ZERO
! 311: CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
! 312: SVA( p ) = TEMP1*DSQRT( AAPP )*D( p )
! 313: END IF
! 314: AAPP = SVA( p )
! 315: ELSE
! 316: AAPP = SVA( p )
! 317: END IF
! 318:
! 319: *
! 320: IF( AAPP.GT.ZERO ) THEN
! 321: *
! 322: PSKIPPED = 0
! 323: *
! 324: DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
! 325: *
! 326: AAQQ = SVA( q )
! 327:
! 328: IF( AAQQ.GT.ZERO ) THEN
! 329: *
! 330: AAPP0 = AAPP
! 331: IF( AAQQ.GE.ONE ) THEN
! 332: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 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: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 346: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 347: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 348: + q ), 1 )*D( p )*D( q ) / AAQQ )
! 349: + / AAPP
! 350: ELSE
! 351: CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
! 352: CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
! 353: + M, 1, WORK, LDA, IERR )
! 354: AAPQ = DDOT( M, WORK, 1, A( 1, p ),
! 355: + 1 )*D( p ) / AAPP
! 356: END IF
! 357: END IF
! 358: *
! 359: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
! 360: *
! 361: * TO rotate or NOT to rotate, THAT is the question ...
! 362: *
! 363: IF( DABS( AAPQ ).GT.TOL ) THEN
! 364: *
! 365: * .. rotate
! 366: * ROTATED = ROTATED + ONE
! 367: *
! 368: IF( ir1.EQ.0 ) THEN
! 369: NOTROT = 0
! 370: PSKIPPED = 0
! 371: ISWROT = ISWROT + 1
! 372: END IF
! 373: *
! 374: IF( ROTOK ) THEN
! 375: *
! 376: AQOAP = AAQQ / AAPP
! 377: APOAQ = AAPP / AAQQ
! 378: THETA = -HALF*DABS( AQOAP-APOAQ ) /
! 379: + AAPQ
! 380: *
! 381: IF( DABS( THETA ).GT.BIGTHETA ) THEN
! 382: *
! 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( ONE-T*AQOAP*
! 395: + AAPQ )
! 396: MXSINJ = DMAX1( MXSINJ, DABS( T ) )
! 397: *
! 398: ELSE
! 399: *
! 400: * .. choose correct signum for THETA and rotate
! 401: *
! 402: THSIGN = -DSIGN( ONE, AAPQ )
! 403: T = ONE / ( THETA+THSIGN*
! 404: + DSQRT( ONE+THETA*THETA ) )
! 405: CS = DSQRT( ONE / ( ONE+T*T ) )
! 406: SN = T*CS
! 407: *
! 408: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
! 409: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 410: + ONE+T*APOAQ*AAPQ ) )
! 411: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 412: + ONE-T*AQOAP*AAPQ ) )
! 413: *
! 414: APOAQ = D( p ) / D( q )
! 415: AQOAP = D( q ) / D( p )
! 416: IF( D( p ).GE.ONE ) THEN
! 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: D( p ) = D( p )*CS
! 436: D( q ) = D( q ) / CS
! 437: IF( RSVEC ) THEN
! 438: CALL DAXPY( MVL, -T*AQOAP,
! 439: + V( 1, q ), 1,
! 440: + V( 1, p ), 1 )
! 441: CALL DAXPY( MVL,
! 442: + CS*SN*APOAQ,
! 443: + V( 1, p ), 1,
! 444: + V( 1, q ), 1 )
! 445: END IF
! 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: D( p ) = D( p ) / CS
! 456: D( q ) = D( q )*CS
! 457: IF( RSVEC ) THEN
! 458: CALL DAXPY( MVL, T*APOAQ,
! 459: + V( 1, p ), 1,
! 460: + V( 1, q ), 1 )
! 461: CALL DAXPY( MVL,
! 462: + -CS*SN*AQOAP,
! 463: + V( 1, q ), 1,
! 464: + V( 1, p ), 1 )
! 465: END IF
! 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: * .. have to use modified Gram-Schmidt like transformation
! 512: CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
! 513: CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
! 514: + 1, WORK, LDA, IERR )
! 515: CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
! 516: + 1, A( 1, q ), LDA, IERR )
! 517: TEMP1 = -AAPQ*D( p ) / D( q )
! 518: CALL DAXPY( M, TEMP1, WORK, 1,
! 519: + A( 1, q ), 1 )
! 520: CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
! 521: + 1, A( 1, q ), LDA, IERR )
! 522: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 523: + ONE-AAPQ*AAPQ ) )
! 524: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 525: END IF
! 526: * END IF ROTOK THEN ... ELSE
! 527: *
! 528: * In the case of cancellation in updating SVA(q), SVA(p)
! 529: * recompute SVA(q), SVA(p).
! 530: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 531: + THEN
! 532: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 533: + ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 534: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
! 535: + D( q )
! 536: ELSE
! 537: T = ZERO
! 538: AAQQ = ZERO
! 539: CALL DLASSQ( M, A( 1, q ), 1, T,
! 540: + AAQQ )
! 541: SVA( q ) = T*DSQRT( AAQQ )*D( q )
! 542: END IF
! 543: END IF
! 544: IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
! 545: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 546: + ( AAPP.GT.ROOTSFMIN ) ) THEN
! 547: AAPP = DNRM2( M, A( 1, p ), 1 )*
! 548: + D( p )
! 549: ELSE
! 550: T = ZERO
! 551: AAPP = ZERO
! 552: CALL DLASSQ( M, A( 1, p ), 1, T,
! 553: + AAPP )
! 554: AAPP = T*DSQRT( AAPP )*D( p )
! 555: END IF
! 556: SVA( p ) = AAPP
! 557: END IF
! 558: *
! 559: ELSE
! 560: * A(:,p) and A(:,q) already numerically orthogonal
! 561: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
! 562: PSKIPPED = PSKIPPED + 1
! 563: END IF
! 564: ELSE
! 565: * A(:,q) is zero column
! 566: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
! 567: PSKIPPED = PSKIPPED + 1
! 568: END IF
! 569: *
! 570: IF( ( i.LE.SWBAND ) .AND.
! 571: + ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 572: IF( ir1.EQ.0 )AAPP = -AAPP
! 573: NOTROT = 0
! 574: GO TO 2103
! 575: END IF
! 576: *
! 577: 2002 CONTINUE
! 578: * END q-LOOP
! 579: *
! 580: 2103 CONTINUE
! 581: * bailed out of q-loop
! 582:
! 583: SVA( p ) = AAPP
! 584:
! 585: ELSE
! 586: SVA( p ) = AAPP
! 587: IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
! 588: + NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
! 589: END IF
! 590: *
! 591: 2001 CONTINUE
! 592: * end of the p-loop
! 593: * end of doing the block ( ibr, ibr )
! 594: 1002 CONTINUE
! 595: * end of ir1-loop
! 596: *
! 597: *........................................................
! 598: * ... go to the off diagonal blocks
! 599: *
! 600: igl = ( ibr-1 )*KBL + 1
! 601: *
! 602: DO 2010 jbc = ibr + 1, NBL
! 603: *
! 604: jgl = ( jbc-1 )*KBL + 1
! 605: *
! 606: * doing the block at ( ibr, jbc )
! 607: *
! 608: IJBLSK = 0
! 609: DO 2100 p = igl, MIN0( igl+KBL-1, N )
! 610: *
! 611: AAPP = SVA( p )
! 612: *
! 613: IF( AAPP.GT.ZERO ) THEN
! 614: *
! 615: PSKIPPED = 0
! 616: *
! 617: DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
! 618: *
! 619: AAQQ = SVA( q )
! 620: *
! 621: IF( AAQQ.GT.ZERO ) THEN
! 622: AAPP0 = AAPP
! 623: *
! 624: * -#- M x 2 Jacobi SVD -#-
! 625: *
! 626: * -#- Safe Gram matrix computation -#-
! 627: *
! 628: IF( AAQQ.GE.ONE ) THEN
! 629: IF( AAPP.GE.AAQQ ) THEN
! 630: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 631: ELSE
! 632: ROTOK = ( SMALL*AAQQ ).LE.AAPP
! 633: END IF
! 634: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
! 635: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 636: + q ), 1 )*D( p )*D( q ) / AAQQ )
! 637: + / AAPP
! 638: ELSE
! 639: CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
! 640: CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
! 641: + M, 1, WORK, LDA, IERR )
! 642: AAPQ = DDOT( M, WORK, 1, A( 1, q ),
! 643: + 1 )*D( q ) / AAQQ
! 644: END IF
! 645: ELSE
! 646: IF( AAPP.GE.AAQQ ) THEN
! 647: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 648: ELSE
! 649: ROTOK = AAQQ.LE.( AAPP / SMALL )
! 650: END IF
! 651: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 652: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 653: + q ), 1 )*D( p )*D( q ) / AAQQ )
! 654: + / AAPP
! 655: ELSE
! 656: CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
! 657: CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
! 658: + M, 1, WORK, LDA, IERR )
! 659: AAPQ = DDOT( M, WORK, 1, A( 1, p ),
! 660: + 1 )*D( p ) / AAPP
! 661: END IF
! 662: END IF
! 663: *
! 664: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
! 665: *
! 666: * TO rotate or NOT to rotate, THAT is the question ...
! 667: *
! 668: IF( DABS( AAPQ ).GT.TOL ) THEN
! 669: NOTROT = 0
! 670: * ROTATED = ROTATED + 1
! 671: PSKIPPED = 0
! 672: ISWROT = ISWROT + 1
! 673: *
! 674: IF( ROTOK ) THEN
! 675: *
! 676: AQOAP = AAQQ / AAPP
! 677: APOAQ = AAPP / AAQQ
! 678: THETA = -HALF*DABS( AQOAP-APOAQ ) /
! 679: + AAPQ
! 680: IF( AAQQ.GT.AAPP0 )THETA = -THETA
! 681: *
! 682: IF( DABS( THETA ).GT.BIGTHETA ) THEN
! 683: T = HALF / THETA
! 684: FASTR( 3 ) = T*D( p ) / D( q )
! 685: FASTR( 4 ) = -T*D( q ) / D( p )
! 686: CALL DROTM( M, A( 1, p ), 1,
! 687: + A( 1, q ), 1, FASTR )
! 688: IF( RSVEC )CALL DROTM( MVL,
! 689: + V( 1, p ), 1,
! 690: + V( 1, q ), 1,
! 691: + FASTR )
! 692: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 693: + ONE+T*APOAQ*AAPQ ) )
! 694: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 695: + ONE-T*AQOAP*AAPQ ) )
! 696: MXSINJ = DMAX1( MXSINJ, DABS( T ) )
! 697: ELSE
! 698: *
! 699: * .. choose correct signum for THETA and rotate
! 700: *
! 701: THSIGN = -DSIGN( ONE, AAPQ )
! 702: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
! 703: T = ONE / ( THETA+THSIGN*
! 704: + DSQRT( ONE+THETA*THETA ) )
! 705: CS = DSQRT( ONE / ( ONE+T*T ) )
! 706: SN = T*CS
! 707: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
! 708: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 709: + ONE+T*APOAQ*AAPQ ) )
! 710: AAPP = AAPP*DSQRT( ONE-T*AQOAP*
! 711: + AAPQ )
! 712: *
! 713: APOAQ = D( p ) / D( q )
! 714: AQOAP = D( q ) / D( p )
! 715: IF( D( p ).GE.ONE ) THEN
! 716: *
! 717: IF( D( q ).GE.ONE ) THEN
! 718: FASTR( 3 ) = T*APOAQ
! 719: FASTR( 4 ) = -T*AQOAP
! 720: D( p ) = D( p )*CS
! 721: D( q ) = D( q )*CS
! 722: CALL DROTM( M, A( 1, p ), 1,
! 723: + A( 1, q ), 1,
! 724: + FASTR )
! 725: IF( RSVEC )CALL DROTM( MVL,
! 726: + V( 1, p ), 1, V( 1, q ),
! 727: + 1, FASTR )
! 728: ELSE
! 729: CALL DAXPY( M, -T*AQOAP,
! 730: + A( 1, q ), 1,
! 731: + A( 1, p ), 1 )
! 732: CALL DAXPY( M, CS*SN*APOAQ,
! 733: + A( 1, p ), 1,
! 734: + A( 1, q ), 1 )
! 735: IF( RSVEC ) THEN
! 736: CALL DAXPY( MVL, -T*AQOAP,
! 737: + V( 1, q ), 1,
! 738: + V( 1, p ), 1 )
! 739: CALL DAXPY( MVL,
! 740: + CS*SN*APOAQ,
! 741: + V( 1, p ), 1,
! 742: + V( 1, q ), 1 )
! 743: END IF
! 744: D( p ) = D( p )*CS
! 745: D( q ) = D( q ) / CS
! 746: END IF
! 747: ELSE
! 748: IF( D( q ).GE.ONE ) THEN
! 749: CALL DAXPY( M, T*APOAQ,
! 750: + A( 1, p ), 1,
! 751: + A( 1, q ), 1 )
! 752: CALL DAXPY( M, -CS*SN*AQOAP,
! 753: + A( 1, q ), 1,
! 754: + A( 1, p ), 1 )
! 755: IF( RSVEC ) THEN
! 756: CALL DAXPY( MVL, T*APOAQ,
! 757: + V( 1, p ), 1,
! 758: + V( 1, q ), 1 )
! 759: CALL DAXPY( MVL,
! 760: + -CS*SN*AQOAP,
! 761: + V( 1, q ), 1,
! 762: + V( 1, p ), 1 )
! 763: END IF
! 764: D( p ) = D( p ) / CS
! 765: D( q ) = D( q )*CS
! 766: ELSE
! 767: IF( D( p ).GE.D( q ) ) THEN
! 768: CALL DAXPY( M, -T*AQOAP,
! 769: + A( 1, q ), 1,
! 770: + A( 1, p ), 1 )
! 771: CALL DAXPY( M, CS*SN*APOAQ,
! 772: + A( 1, p ), 1,
! 773: + A( 1, q ), 1 )
! 774: D( p ) = D( p )*CS
! 775: D( q ) = D( q ) / CS
! 776: IF( RSVEC ) THEN
! 777: CALL DAXPY( MVL,
! 778: + -T*AQOAP,
! 779: + V( 1, q ), 1,
! 780: + V( 1, p ), 1 )
! 781: CALL DAXPY( MVL,
! 782: + CS*SN*APOAQ,
! 783: + V( 1, p ), 1,
! 784: + V( 1, q ), 1 )
! 785: END IF
! 786: ELSE
! 787: CALL DAXPY( M, T*APOAQ,
! 788: + A( 1, p ), 1,
! 789: + A( 1, q ), 1 )
! 790: CALL DAXPY( M,
! 791: + -CS*SN*AQOAP,
! 792: + A( 1, q ), 1,
! 793: + A( 1, p ), 1 )
! 794: D( p ) = D( p ) / CS
! 795: D( q ) = D( q )*CS
! 796: IF( RSVEC ) THEN
! 797: CALL DAXPY( MVL,
! 798: + T*APOAQ, V( 1, p ),
! 799: + 1, V( 1, q ), 1 )
! 800: CALL DAXPY( MVL,
! 801: + -CS*SN*AQOAP,
! 802: + V( 1, q ), 1,
! 803: + V( 1, p ), 1 )
! 804: END IF
! 805: END IF
! 806: END IF
! 807: END IF
! 808: END IF
! 809: *
! 810: ELSE
! 811: IF( AAPP.GT.AAQQ ) THEN
! 812: CALL DCOPY( M, A( 1, p ), 1, WORK,
! 813: + 1 )
! 814: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
! 815: + M, 1, WORK, LDA, IERR )
! 816: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
! 817: + M, 1, A( 1, q ), LDA,
! 818: + IERR )
! 819: TEMP1 = -AAPQ*D( p ) / D( q )
! 820: CALL DAXPY( M, TEMP1, WORK, 1,
! 821: + A( 1, q ), 1 )
! 822: CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
! 823: + M, 1, A( 1, q ), LDA,
! 824: + IERR )
! 825: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 826: + ONE-AAPQ*AAPQ ) )
! 827: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 828: ELSE
! 829: CALL DCOPY( M, A( 1, q ), 1, WORK,
! 830: + 1 )
! 831: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
! 832: + M, 1, WORK, LDA, IERR )
! 833: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
! 834: + M, 1, A( 1, p ), LDA,
! 835: + IERR )
! 836: TEMP1 = -AAPQ*D( q ) / D( p )
! 837: CALL DAXPY( M, TEMP1, WORK, 1,
! 838: + A( 1, p ), 1 )
! 839: CALL DLASCL( 'G', 0, 0, ONE, AAPP,
! 840: + M, 1, A( 1, p ), LDA,
! 841: + IERR )
! 842: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
! 843: + ONE-AAPQ*AAPQ ) )
! 844: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 845: END IF
! 846: END IF
! 847: * END IF ROTOK THEN ... ELSE
! 848: *
! 849: * In the case of cancellation in updating SVA(q)
! 850: * .. recompute SVA(q)
! 851: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 852: + THEN
! 853: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 854: + ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 855: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
! 856: + D( q )
! 857: ELSE
! 858: T = ZERO
! 859: AAQQ = ZERO
! 860: CALL DLASSQ( M, A( 1, q ), 1, T,
! 861: + AAQQ )
! 862: SVA( q ) = T*DSQRT( AAQQ )*D( q )
! 863: END IF
! 864: END IF
! 865: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
! 866: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 867: + ( AAPP.GT.ROOTSFMIN ) ) THEN
! 868: AAPP = DNRM2( M, A( 1, p ), 1 )*
! 869: + D( p )
! 870: ELSE
! 871: T = ZERO
! 872: AAPP = ZERO
! 873: CALL DLASSQ( M, A( 1, p ), 1, T,
! 874: + AAPP )
! 875: AAPP = T*DSQRT( AAPP )*D( p )
! 876: END IF
! 877: SVA( p ) = AAPP
! 878: END IF
! 879: * end of OK rotation
! 880: ELSE
! 881: NOTROT = NOTROT + 1
! 882: PSKIPPED = PSKIPPED + 1
! 883: IJBLSK = IJBLSK + 1
! 884: END IF
! 885: ELSE
! 886: NOTROT = NOTROT + 1
! 887: PSKIPPED = PSKIPPED + 1
! 888: IJBLSK = IJBLSK + 1
! 889: END IF
! 890: *
! 891: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
! 892: + THEN
! 893: SVA( p ) = AAPP
! 894: NOTROT = 0
! 895: GO TO 2011
! 896: END IF
! 897: IF( ( i.LE.SWBAND ) .AND.
! 898: + ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 899: AAPP = -AAPP
! 900: NOTROT = 0
! 901: GO TO 2203
! 902: END IF
! 903: *
! 904: 2200 CONTINUE
! 905: * end of the q-loop
! 906: 2203 CONTINUE
! 907: *
! 908: SVA( p ) = AAPP
! 909: *
! 910: ELSE
! 911: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
! 912: + MIN0( jgl+KBL-1, N ) - jgl + 1
! 913: IF( AAPP.LT.ZERO )NOTROT = 0
! 914: END IF
! 915:
! 916: 2100 CONTINUE
! 917: * end of the p-loop
! 918: 2010 CONTINUE
! 919: * end of the jbc-loop
! 920: 2011 CONTINUE
! 921: *2011 bailed out of the jbc-loop
! 922: DO 2012 p = igl, MIN0( igl+KBL-1, N )
! 923: SVA( p ) = DABS( SVA( p ) )
! 924: 2012 CONTINUE
! 925: *
! 926: 2000 CONTINUE
! 927: *2000 :: end of the ibr-loop
! 928: *
! 929: * .. update SVA(N)
! 930: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
! 931: + THEN
! 932: SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
! 933: ELSE
! 934: T = ZERO
! 935: AAPP = ZERO
! 936: CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
! 937: SVA( N ) = T*DSQRT( AAPP )*D( N )
! 938: END IF
! 939: *
! 940: * Additional steering devices
! 941: *
! 942: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
! 943: + ( ISWROT.LE.N ) ) )SWBAND = i
! 944: *
! 945: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
! 946: + ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
! 947: GO TO 1994
! 948: END IF
! 949: *
! 950: IF( NOTROT.GE.EMPTSW )GO TO 1994
! 951:
! 952: 1993 CONTINUE
! 953: * end i=1:NSWEEP loop
! 954: * #:) Reaching this point means that the procedure has comleted the given
! 955: * number of iterations.
! 956: INFO = NSWEEP - 1
! 957: GO TO 1995
! 958: 1994 CONTINUE
! 959: * #:) Reaching this point means that during the i-th sweep all pivots were
! 960: * below the given tolerance, causing early exit.
! 961: *
! 962: INFO = 0
! 963: * #:) INFO = 0 confirms successful iterations.
! 964: 1995 CONTINUE
! 965: *
! 966: * Sort the vector D.
! 967: DO 5991 p = 1, N - 1
! 968: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 969: IF( p.NE.q ) THEN
! 970: TEMP1 = SVA( p )
! 971: SVA( p ) = SVA( q )
! 972: SVA( q ) = TEMP1
! 973: TEMP1 = D( p )
! 974: D( p ) = D( q )
! 975: D( q ) = TEMP1
! 976: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 977: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
! 978: END IF
! 979: 5991 CONTINUE
! 980: *
! 981: RETURN
! 982: * ..
! 983: * .. END OF DGSVJ0
! 984: * ..
! 985: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgsvj0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack