Annotation of rpl/lapack/lapack/dgesvj.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
! 2: + LDV, 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: INTEGER INFO, LDA, LDV, LWORK, M, MV, N
! 22: CHARACTER*1 JOBA, JOBU, JOBV
! 23: * ..
! 24: * .. Array Arguments ..
! 25: DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ),
! 26: + WORK( LWORK )
! 27: * ..
! 28: *
! 29: * Purpose
! 30: * =======
! 31: *
! 32: * DGESVJ computes the singular value decomposition (SVD) of a real
! 33: * M-by-N matrix A, where M >= N. The SVD of A is written as
! 34: * [++] [xx] [x0] [xx]
! 35: * A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx]
! 36: * [++] [xx]
! 37: * where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
! 38: * matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
! 39: * of SIGMA are the singular values of A. The columns of U and V are the
! 40: * left and the right singular vectors of A, respectively.
! 41: *
! 42: * Further Details
! 43: * ~~~~~~~~~~~~~~~
! 44: * The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
! 45: * rotations. The rotations are implemented as fast scaled rotations of
! 46: * Anda and Park [1]. In the case of underflow of the Jacobi angle, a
! 47: * modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
! 48: * column interchanges of de Rijk [2]. The relative accuracy of the computed
! 49: * singular values and the accuracy of the computed singular vectors (in
! 50: * angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
! 51: * The condition number that determines the accuracy in the full rank case
! 52: * is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
! 53: * spectral condition number. The best performance of this Jacobi SVD
! 54: * procedure is achieved if used in an accelerated version of Drmac and
! 55: * Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
! 56: * Some tunning parameters (marked with [TP]) are available for the
! 57: * implementer.
! 58: * The computational range for the nonzero singular values is the machine
! 59: * number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
! 60: * denormalized singular values can be computed with the corresponding
! 61: * gradual loss of accurate digits.
! 62: *
! 63: * Contributors
! 64: * ~~~~~~~~~~~~
! 65: * Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
! 66: *
! 67: * References
! 68: * ~~~~~~~~~~
! 69: * [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
! 70: * SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
! 71: * [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
! 72: * singular value decomposition on a vector computer.
! 73: * SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
! 74: * [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
! 75: * [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
! 76: * value computation in floating point arithmetic.
! 77: * SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
! 78: * [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
! 79: * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
! 80: * LAPACK Working note 169.
! 81: * [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
! 82: * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
! 83: * LAPACK Working note 170.
! 84: * [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
! 85: * QSVD, (H,K)-SVD computations.
! 86: * Department of Mathematics, University of Zagreb, 2008.
! 87: *
! 88: * Bugs, Examples and Comments
! 89: * ~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 90: * Please report all bugs and send interesting test examples and comments to
! 91: * drmac@math.hr. Thank you.
! 92: *
! 93: * Arguments
! 94: * =========
! 95: *
! 96: * JOBA (input) CHARACTER* 1
! 97: * Specifies the structure of A.
! 98: * = 'L': The input matrix A is lower triangular;
! 99: * = 'U': The input matrix A is upper triangular;
! 100: * = 'G': The input matrix A is general M-by-N matrix, M >= N.
! 101: *
! 102: * JOBU (input) CHARACTER*1
! 103: * Specifies whether to compute the left singular vectors
! 104: * (columns of U):
! 105: * = 'U': The left singular vectors corresponding to the nonzero
! 106: * singular values are computed and returned in the leading
! 107: * columns of A. See more details in the description of A.
! 108: * The default numerical orthogonality threshold is set to
! 109: * approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E').
! 110: * = 'C': Analogous to JOBU='U', except that user can control the
! 111: * level of numerical orthogonality of the computed left
! 112: * singular vectors. TOL can be set to TOL = CTOL*EPS, where
! 113: * CTOL is given on input in the array WORK.
! 114: * No CTOL smaller than ONE is allowed. CTOL greater
! 115: * than 1 / EPS is meaningless. The option 'C'
! 116: * can be used if M*EPS is satisfactory orthogonality
! 117: * of the computed left singular vectors, so CTOL=M could
! 118: * save few sweeps of Jacobi rotations.
! 119: * See the descriptions of A and WORK(1).
! 120: * = 'N': The matrix U is not computed. However, see the
! 121: * description of A.
! 122: *
! 123: * JOBV (input) CHARACTER*1
! 124: * Specifies whether to compute the right singular vectors, that
! 125: * is, the matrix V:
! 126: * = 'V' : the matrix V is computed and returned in the array V
! 127: * = 'A' : the Jacobi rotations are applied to the MV-by-N
! 128: * array V. In other words, the right singular vector
! 129: * matrix V is not computed explicitly, instead it is
! 130: * applied to an MV-by-N matrix initially stored in the
! 131: * first MV rows of V.
! 132: * = 'N' : the matrix V is not computed and the array V is not
! 133: * referenced
! 134: *
! 135: * M (input) INTEGER
! 136: * The number of rows of the input matrix A. M >= 0.
! 137: *
! 138: * N (input) INTEGER
! 139: * The number of columns of the input matrix A.
! 140: * M >= N >= 0.
! 141: *
! 142: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! 143: * On entry, the M-by-N matrix A.
! 144: * On exit :
! 145: * If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :
! 146: * If INFO .EQ. 0 :
! 147: * RANKA orthonormal columns of U are returned in the
! 148: * leading RANKA columns of the array A. Here RANKA <= N
! 149: * is the number of computed singular values of A that are
! 150: * above the underflow threshold DLAMCH('S'). The singular
! 151: * vectors corresponding to underflowed or zero singular
! 152: * values are not computed. The value of RANKA is returned
! 153: * in the array WORK as RANKA=NINT(WORK(2)). Also see the
! 154: * descriptions of SVA and WORK. The computed columns of U
! 155: * are mutually numerically orthogonal up to approximately
! 156: * TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
! 157: * see the description of JOBU.
! 158: * If INFO .GT. 0 :
! 159: * the procedure DGESVJ did not converge in the given number
! 160: * of iterations (sweeps). In that case, the computed
! 161: * columns of U may not be orthogonal up to TOL. The output
! 162: * U (stored in A), SIGMA (given by the computed singular
! 163: * values in SVA(1:N)) and V is still a decomposition of the
! 164: * input matrix A in the sense that the residual
! 165: * ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
! 166: *
! 167: * If JOBU .EQ. 'N' :
! 168: * If INFO .EQ. 0 :
! 169: * Note that the left singular vectors are 'for free' in the
! 170: * one-sided Jacobi SVD algorithm. However, if only the
! 171: * singular values are needed, the level of numerical
! 172: * orthogonality of U is not an issue and iterations are
! 173: * stopped when the columns of the iterated matrix are
! 174: * numerically orthogonal up to approximately M*EPS. Thus,
! 175: * on exit, A contains the columns of U scaled with the
! 176: * corresponding singular values.
! 177: * If INFO .GT. 0 :
! 178: * the procedure DGESVJ did not converge in the given number
! 179: * of iterations (sweeps).
! 180: *
! 181: * LDA (input) INTEGER
! 182: * The leading dimension of the array A. LDA >= max(1,M).
! 183: *
! 184: * SVA (workspace/output) DOUBLE PRECISION array, dimension (N)
! 185: * On exit :
! 186: * If INFO .EQ. 0 :
! 187: * depending on the value SCALE = WORK(1), we have:
! 188: * If SCALE .EQ. ONE :
! 189: * SVA(1:N) contains the computed singular values of A.
! 190: * During the computation SVA contains the Euclidean column
! 191: * norms of the iterated matrices in the array A.
! 192: * If SCALE .NE. ONE :
! 193: * The singular values of A are SCALE*SVA(1:N), and this
! 194: * factored representation is due to the fact that some of the
! 195: * singular values of A might underflow or overflow.
! 196: * If INFO .GT. 0 :
! 197: * the procedure DGESVJ did not converge in the given number of
! 198: * iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
! 199: *
! 200: * MV (input) INTEGER
! 201: * If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ
! 202: * is applied to the first MV rows of V. See the description of JOBV.
! 203: *
! 204: * V (input/output) DOUBLE PRECISION array, dimension (LDV,N)
! 205: * If JOBV = 'V', then V contains on exit the N-by-N matrix of
! 206: * the right singular vectors;
! 207: * If JOBV = 'A', then V contains the product of the computed right
! 208: * singular vector matrix and the initial matrix in
! 209: * the array V.
! 210: * If JOBV = 'N', then V is not referenced.
! 211: *
! 212: * LDV (input) INTEGER
! 213: * The leading dimension of the array V, LDV .GE. 1.
! 214: * If JOBV .EQ. 'V', then LDV .GE. max(1,N).
! 215: * If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
! 216: *
! 217: * WORK (input/workspace/output) DOUBLE PRECISION array, dimension max(4,M+N).
! 218: * On entry :
! 219: * If JOBU .EQ. 'C' :
! 220: * WORK(1) = CTOL, where CTOL defines the threshold for convergence.
! 221: * The process stops if all columns of A are mutually
! 222: * orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
! 223: * It is required that CTOL >= ONE, i.e. it is not
! 224: * allowed to force the routine to obtain orthogonality
! 225: * below EPSILON.
! 226: * On exit :
! 227: * WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
! 228: * are the computed singular values of A.
! 229: * (See description of SVA().)
! 230: * WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
! 231: * singular values.
! 232: * WORK(3) = NINT(WORK(3)) is the number of the computed singular
! 233: * values that are larger than the underflow threshold.
! 234: * WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
! 235: * rotations needed for numerical convergence.
! 236: * WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
! 237: * This is useful information in cases when DGESVJ did
! 238: * not converge, as it can be used to estimate whether
! 239: * the output is stil useful and for post festum analysis.
! 240: * WORK(6) = the largest absolute value over all sines of the
! 241: * Jacobi rotation angles in the last sweep. It can be
! 242: * useful for a post festum analysis.
! 243: *
! 244: * LWORK (input) INTEGER
! 245: * length of WORK, WORK >= MAX(6,M+N)
! 246: *
! 247: * INFO (output) INTEGER
! 248: * = 0 : successful exit.
! 249: * < 0 : if INFO = -i, then the i-th argument had an illegal value
! 250: * > 0 : DGESVJ did not converge in the maximal allowed number (30)
! 251: * of sweeps. The output may still be useful. See the
! 252: * description of WORK.
! 253: *
! 254: * =====================================================================
! 255: *
! 256: * .. Local Parameters ..
! 257: DOUBLE PRECISION ZERO, HALF, ONE, TWO
! 258: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
! 259: + TWO = 2.0D0 )
! 260: INTEGER NSWEEP
! 261: PARAMETER ( NSWEEP = 30 )
! 262: * ..
! 263: * .. Local Scalars ..
! 264: DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
! 265: + BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,
! 266: + MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
! 267: + SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,
! 268: + THSIGN, TOL
! 269: INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
! 270: + ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
! 271: + N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
! 272: + SWBAND
! 273: LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
! 274: + RSVEC, UCTOL, UPPER
! 275: * ..
! 276: * .. Local Arrays ..
! 277: DOUBLE PRECISION FASTR( 5 )
! 278: * ..
! 279: * .. Intrinsic Functions ..
! 280: INTRINSIC DABS, DMAX1, DMIN1, DBLE, MIN0, DSIGN, DSQRT
! 281: * ..
! 282: * .. External Functions ..
! 283: * ..
! 284: * from BLAS
! 285: DOUBLE PRECISION DDOT, DNRM2
! 286: EXTERNAL DDOT, DNRM2
! 287: INTEGER IDAMAX
! 288: EXTERNAL IDAMAX
! 289: * from LAPACK
! 290: DOUBLE PRECISION DLAMCH
! 291: EXTERNAL DLAMCH
! 292: LOGICAL LSAME
! 293: EXTERNAL LSAME
! 294: * ..
! 295: * .. External Subroutines ..
! 296: * ..
! 297: * from BLAS
! 298: EXTERNAL DAXPY, DCOPY, DROTM, DSCAL, DSWAP
! 299: * from LAPACK
! 300: EXTERNAL DLASCL, DLASET, DLASSQ, XERBLA
! 301: *
! 302: EXTERNAL DGSVJ0, DGSVJ1
! 303: * ..
! 304: * .. Executable Statements ..
! 305: *
! 306: * Test the input arguments
! 307: *
! 308: LSVEC = LSAME( JOBU, 'U' )
! 309: UCTOL = LSAME( JOBU, 'C' )
! 310: RSVEC = LSAME( JOBV, 'V' )
! 311: APPLV = LSAME( JOBV, 'A' )
! 312: UPPER = LSAME( JOBA, 'U' )
! 313: LOWER = LSAME( JOBA, 'L' )
! 314: *
! 315: IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
! 316: INFO = -1
! 317: ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
! 318: INFO = -2
! 319: ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
! 320: INFO = -3
! 321: ELSE IF( M.LT.0 ) THEN
! 322: INFO = -4
! 323: ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
! 324: INFO = -5
! 325: ELSE IF( LDA.LT.M ) THEN
! 326: INFO = -7
! 327: ELSE IF( MV.LT.0 ) THEN
! 328: INFO = -9
! 329: ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
! 330: + ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
! 331: INFO = -11
! 332: ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
! 333: INFO = -12
! 334: ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN
! 335: INFO = -13
! 336: ELSE
! 337: INFO = 0
! 338: END IF
! 339: *
! 340: * #:(
! 341: IF( INFO.NE.0 ) THEN
! 342: CALL XERBLA( 'DGESVJ', -INFO )
! 343: RETURN
! 344: END IF
! 345: *
! 346: * #:) Quick return for void matrix
! 347: *
! 348: IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
! 349: *
! 350: * Set numerical parameters
! 351: * The stopping criterion for Jacobi rotations is
! 352: *
! 353: * max_{i<>j}|A(:,i)^T * A(:,j)|/(||A(:,i)||*||A(:,j)||) < CTOL*EPS
! 354: *
! 355: * where EPS is the round-off and CTOL is defined as follows:
! 356: *
! 357: IF( UCTOL ) THEN
! 358: * ... user controlled
! 359: CTOL = WORK( 1 )
! 360: ELSE
! 361: * ... default
! 362: IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
! 363: CTOL = DSQRT( DBLE( M ) )
! 364: ELSE
! 365: CTOL = DBLE( M )
! 366: END IF
! 367: END IF
! 368: * ... and the machine dependent parameters are
! 369: *[!] (Make sure that DLAMCH() works properly on the target machine.)
! 370: *
! 371: EPSILON = DLAMCH( 'Epsilon' )
! 372: ROOTEPS = DSQRT( EPSILON )
! 373: SFMIN = DLAMCH( 'SafeMinimum' )
! 374: ROOTSFMIN = DSQRT( SFMIN )
! 375: SMALL = SFMIN / EPSILON
! 376: BIG = DLAMCH( 'Overflow' )
! 377: * BIG = ONE / SFMIN
! 378: ROOTBIG = ONE / ROOTSFMIN
! 379: LARGE = BIG / DSQRT( DBLE( M*N ) )
! 380: BIGTHETA = ONE / ROOTEPS
! 381: *
! 382: TOL = CTOL*EPSILON
! 383: ROOTTOL = DSQRT( TOL )
! 384: *
! 385: IF( DBLE( M )*EPSILON.GE.ONE ) THEN
! 386: INFO = -5
! 387: CALL XERBLA( 'DGESVJ', -INFO )
! 388: RETURN
! 389: END IF
! 390: *
! 391: * Initialize the right singular vector matrix.
! 392: *
! 393: IF( RSVEC ) THEN
! 394: MVL = N
! 395: CALL DLASET( 'A', MVL, N, ZERO, ONE, V, LDV )
! 396: ELSE IF( APPLV ) THEN
! 397: MVL = MV
! 398: END IF
! 399: RSVEC = RSVEC .OR. APPLV
! 400: *
! 401: * Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
! 402: *(!) If necessary, scale A to protect the largest singular value
! 403: * from overflow. It is possible that saving the largest singular
! 404: * value destroys the information about the small ones.
! 405: * This initial scaling is almost minimal in the sense that the
! 406: * goal is to make sure that no column norm overflows, and that
! 407: * DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
! 408: * in A are detected, the procedure returns with INFO=-6.
! 409: *
! 410: SCALE = ONE / DSQRT( DBLE( M )*DBLE( N ) )
! 411: NOSCALE = .TRUE.
! 412: GOSCALE = .TRUE.
! 413: *
! 414: IF( LOWER ) THEN
! 415: * the input matrix is M-by-N lower triangular (trapezoidal)
! 416: DO 1874 p = 1, N
! 417: AAPP = ZERO
! 418: AAQQ = ZERO
! 419: CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
! 420: IF( AAPP.GT.BIG ) THEN
! 421: INFO = -6
! 422: CALL XERBLA( 'DGESVJ', -INFO )
! 423: RETURN
! 424: END IF
! 425: AAQQ = DSQRT( AAQQ )
! 426: IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
! 427: SVA( p ) = AAPP*AAQQ
! 428: ELSE
! 429: NOSCALE = .FALSE.
! 430: SVA( p ) = AAPP*( AAQQ*SCALE )
! 431: IF( GOSCALE ) THEN
! 432: GOSCALE = .FALSE.
! 433: DO 1873 q = 1, p - 1
! 434: SVA( q ) = SVA( q )*SCALE
! 435: 1873 CONTINUE
! 436: END IF
! 437: END IF
! 438: 1874 CONTINUE
! 439: ELSE IF( UPPER ) THEN
! 440: * the input matrix is M-by-N upper triangular (trapezoidal)
! 441: DO 2874 p = 1, N
! 442: AAPP = ZERO
! 443: AAQQ = ZERO
! 444: CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
! 445: IF( AAPP.GT.BIG ) THEN
! 446: INFO = -6
! 447: CALL XERBLA( 'DGESVJ', -INFO )
! 448: RETURN
! 449: END IF
! 450: AAQQ = DSQRT( AAQQ )
! 451: IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
! 452: SVA( p ) = AAPP*AAQQ
! 453: ELSE
! 454: NOSCALE = .FALSE.
! 455: SVA( p ) = AAPP*( AAQQ*SCALE )
! 456: IF( GOSCALE ) THEN
! 457: GOSCALE = .FALSE.
! 458: DO 2873 q = 1, p - 1
! 459: SVA( q ) = SVA( q )*SCALE
! 460: 2873 CONTINUE
! 461: END IF
! 462: END IF
! 463: 2874 CONTINUE
! 464: ELSE
! 465: * the input matrix is M-by-N general dense
! 466: DO 3874 p = 1, N
! 467: AAPP = ZERO
! 468: AAQQ = ZERO
! 469: CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
! 470: IF( AAPP.GT.BIG ) THEN
! 471: INFO = -6
! 472: CALL XERBLA( 'DGESVJ', -INFO )
! 473: RETURN
! 474: END IF
! 475: AAQQ = DSQRT( AAQQ )
! 476: IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
! 477: SVA( p ) = AAPP*AAQQ
! 478: ELSE
! 479: NOSCALE = .FALSE.
! 480: SVA( p ) = AAPP*( AAQQ*SCALE )
! 481: IF( GOSCALE ) THEN
! 482: GOSCALE = .FALSE.
! 483: DO 3873 q = 1, p - 1
! 484: SVA( q ) = SVA( q )*SCALE
! 485: 3873 CONTINUE
! 486: END IF
! 487: END IF
! 488: 3874 CONTINUE
! 489: END IF
! 490: *
! 491: IF( NOSCALE )SCALE = ONE
! 492: *
! 493: * Move the smaller part of the spectrum from the underflow threshold
! 494: *(!) Start by determining the position of the nonzero entries of the
! 495: * array SVA() relative to ( SFMIN, BIG ).
! 496: *
! 497: AAPP = ZERO
! 498: AAQQ = BIG
! 499: DO 4781 p = 1, N
! 500: IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) )
! 501: AAPP = DMAX1( AAPP, SVA( p ) )
! 502: 4781 CONTINUE
! 503: *
! 504: * #:) Quick return for zero matrix
! 505: *
! 506: IF( AAPP.EQ.ZERO ) THEN
! 507: IF( LSVEC )CALL DLASET( 'G', M, N, ZERO, ONE, A, LDA )
! 508: WORK( 1 ) = ONE
! 509: WORK( 2 ) = ZERO
! 510: WORK( 3 ) = ZERO
! 511: WORK( 4 ) = ZERO
! 512: WORK( 5 ) = ZERO
! 513: WORK( 6 ) = ZERO
! 514: RETURN
! 515: END IF
! 516: *
! 517: * #:) Quick return for one-column matrix
! 518: *
! 519: IF( N.EQ.1 ) THEN
! 520: IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
! 521: + A( 1, 1 ), LDA, IERR )
! 522: WORK( 1 ) = ONE / SCALE
! 523: IF( SVA( 1 ).GE.SFMIN ) THEN
! 524: WORK( 2 ) = ONE
! 525: ELSE
! 526: WORK( 2 ) = ZERO
! 527: END IF
! 528: WORK( 3 ) = ZERO
! 529: WORK( 4 ) = ZERO
! 530: WORK( 5 ) = ZERO
! 531: WORK( 6 ) = ZERO
! 532: RETURN
! 533: END IF
! 534: *
! 535: * Protect small singular values from underflow, and try to
! 536: * avoid underflows/overflows in computing Jacobi rotations.
! 537: *
! 538: SN = DSQRT( SFMIN / EPSILON )
! 539: TEMP1 = DSQRT( BIG / DBLE( N ) )
! 540: IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
! 541: + ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
! 542: TEMP1 = DMIN1( BIG, TEMP1 / AAPP )
! 543: * AAQQ = AAQQ*TEMP1
! 544: * AAPP = AAPP*TEMP1
! 545: ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
! 546: TEMP1 = DMIN1( SN / AAQQ, BIG / ( AAPP*DSQRT( DBLE( N ) ) ) )
! 547: * AAQQ = AAQQ*TEMP1
! 548: * AAPP = AAPP*TEMP1
! 549: ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
! 550: TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP )
! 551: * AAQQ = AAQQ*TEMP1
! 552: * AAPP = AAPP*TEMP1
! 553: ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
! 554: TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) )
! 555: * AAQQ = AAQQ*TEMP1
! 556: * AAPP = AAPP*TEMP1
! 557: ELSE
! 558: TEMP1 = ONE
! 559: END IF
! 560: *
! 561: * Scale, if necessary
! 562: *
! 563: IF( TEMP1.NE.ONE ) THEN
! 564: CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
! 565: END IF
! 566: SCALE = TEMP1*SCALE
! 567: IF( SCALE.NE.ONE ) THEN
! 568: CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
! 569: SCALE = ONE / SCALE
! 570: END IF
! 571: *
! 572: * Row-cyclic Jacobi SVD algorithm with column pivoting
! 573: *
! 574: EMPTSW = ( N*( N-1 ) ) / 2
! 575: NOTROT = 0
! 576: FASTR( 1 ) = ZERO
! 577: *
! 578: * A is represented in factored form A = A * diag(WORK), where diag(WORK)
! 579: * is initialized to identity. WORK is updated during fast scaled
! 580: * rotations.
! 581: *
! 582: DO 1868 q = 1, N
! 583: WORK( q ) = ONE
! 584: 1868 CONTINUE
! 585: *
! 586: *
! 587: SWBAND = 3
! 588: *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
! 589: * if DGESVJ is used as a computational routine in the preconditioned
! 590: * Jacobi SVD algorithm DGESVJ. For sweeps i=1:SWBAND the procedure
! 591: * works on pivots inside a band-like region around the diagonal.
! 592: * The boundaries are determined dynamically, based on the number of
! 593: * pivots above a threshold.
! 594: *
! 595: KBL = MIN0( 8, N )
! 596: *[TP] KBL is a tuning parameter that defines the tile size in the
! 597: * tiling of the p-q loops of pivot pairs. In general, an optimal
! 598: * value of KBL depends on the matrix dimensions and on the
! 599: * parameters of the computer's memory.
! 600: *
! 601: NBL = N / KBL
! 602: IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
! 603: *
! 604: BLSKIP = KBL**2
! 605: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
! 606: *
! 607: ROWSKIP = MIN0( 5, KBL )
! 608: *[TP] ROWSKIP is a tuning parameter.
! 609: *
! 610: LKAHEAD = 1
! 611: *[TP] LKAHEAD is a tuning parameter.
! 612: *
! 613: * Quasi block transformations, using the lower (upper) triangular
! 614: * structure of the input matrix. The quasi-block-cycling usually
! 615: * invokes cubic convergence. Big part of this cycle is done inside
! 616: * canonical subspaces of dimensions less than M.
! 617: *
! 618: IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN
! 619: *[TP] The number of partition levels and the actual partition are
! 620: * tuning parameters.
! 621: N4 = N / 4
! 622: N2 = N / 2
! 623: N34 = 3*N4
! 624: IF( APPLV ) THEN
! 625: q = 0
! 626: ELSE
! 627: q = 1
! 628: END IF
! 629: *
! 630: IF( LOWER ) THEN
! 631: *
! 632: * This works very well on lower triangular matrices, in particular
! 633: * in the framework of the preconditioned Jacobi SVD (xGEJSV).
! 634: * The idea is simple:
! 635: * [+ 0 0 0] Note that Jacobi transformations of [0 0]
! 636: * [+ + 0 0] [0 0]
! 637: * [+ + x 0] actually work on [x 0] [x 0]
! 638: * [+ + x x] [x x]. [x x]
! 639: *
! 640: CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
! 641: + WORK( N34+1 ), SVA( N34+1 ), MVL,
! 642: + V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,
! 643: + 2, WORK( N+1 ), LWORK-N, IERR )
! 644: *
! 645: CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
! 646: + WORK( N2+1 ), SVA( N2+1 ), MVL,
! 647: + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,
! 648: + WORK( N+1 ), LWORK-N, IERR )
! 649: *
! 650: CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
! 651: + WORK( N2+1 ), SVA( N2+1 ), MVL,
! 652: + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
! 653: + WORK( N+1 ), LWORK-N, IERR )
! 654: *
! 655: CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
! 656: + WORK( N4+1 ), SVA( N4+1 ), MVL,
! 657: + V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,
! 658: + WORK( N+1 ), LWORK-N, IERR )
! 659: *
! 660: CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
! 661: + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
! 662: + IERR )
! 663: *
! 664: CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
! 665: + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
! 666: + LWORK-N, IERR )
! 667: *
! 668: *
! 669: ELSE IF( UPPER ) THEN
! 670: *
! 671: *
! 672: CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
! 673: + EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
! 674: + IERR )
! 675: *
! 676: CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
! 677: + SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
! 678: + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
! 679: + IERR )
! 680: *
! 681: CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
! 682: + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
! 683: + LWORK-N, IERR )
! 684: *
! 685: CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
! 686: + WORK( N2+1 ), SVA( N2+1 ), MVL,
! 687: + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
! 688: + WORK( N+1 ), LWORK-N, IERR )
! 689:
! 690: END IF
! 691: *
! 692: END IF
! 693: *
! 694: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
! 695: *
! 696: DO 1993 i = 1, NSWEEP
! 697: *
! 698: * .. go go go ...
! 699: *
! 700: MXAAPQ = ZERO
! 701: MXSINJ = ZERO
! 702: ISWROT = 0
! 703: *
! 704: NOTROT = 0
! 705: PSKIPPED = 0
! 706: *
! 707: * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
! 708: * 1 <= p < q <= N. This is the first step toward a blocked implementation
! 709: * of the rotations. New implementation, based on block transformations,
! 710: * is under development.
! 711: *
! 712: DO 2000 ibr = 1, NBL
! 713: *
! 714: igl = ( ibr-1 )*KBL + 1
! 715: *
! 716: DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
! 717: *
! 718: igl = igl + ir1*KBL
! 719: *
! 720: DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
! 721: *
! 722: * .. de Rijk's pivoting
! 723: *
! 724: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 725: IF( p.NE.q ) THEN
! 726: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 727: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,
! 728: + V( 1, q ), 1 )
! 729: TEMP1 = SVA( p )
! 730: SVA( p ) = SVA( q )
! 731: SVA( q ) = TEMP1
! 732: TEMP1 = WORK( p )
! 733: WORK( p ) = WORK( q )
! 734: WORK( q ) = TEMP1
! 735: END IF
! 736: *
! 737: IF( ir1.EQ.0 ) THEN
! 738: *
! 739: * Column norms are periodically updated by explicit
! 740: * norm computation.
! 741: * Caveat:
! 742: * Unfortunately, some BLAS implementations compute DNRM2(M,A(1,p),1)
! 743: * as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to
! 744: * overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and to
! 745: * underflow for ||A(:,p)||_2 < DSQRT(underflow_threshold).
! 746: * Hence, DNRM2 cannot be trusted, not even in the case when
! 747: * the true norm is far from the under(over)flow boundaries.
! 748: * If properly implemented DNRM2 is available, the IF-THEN-ELSE
! 749: * below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)".
! 750: *
! 751: IF( ( SVA( p ).LT.ROOTBIG ) .AND.
! 752: + ( SVA( p ).GT.ROOTSFMIN ) ) THEN
! 753: SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )
! 754: ELSE
! 755: TEMP1 = ZERO
! 756: AAPP = ZERO
! 757: CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
! 758: SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )
! 759: END IF
! 760: AAPP = SVA( p )
! 761: ELSE
! 762: AAPP = SVA( p )
! 763: END IF
! 764: *
! 765: IF( AAPP.GT.ZERO ) THEN
! 766: *
! 767: PSKIPPED = 0
! 768: *
! 769: DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
! 770: *
! 771: AAQQ = SVA( q )
! 772: *
! 773: IF( AAQQ.GT.ZERO ) THEN
! 774: *
! 775: AAPP0 = AAPP
! 776: IF( AAQQ.GE.ONE ) THEN
! 777: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 778: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
! 779: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 780: + q ), 1 )*WORK( p )*WORK( q ) /
! 781: + AAQQ ) / AAPP
! 782: ELSE
! 783: CALL DCOPY( M, A( 1, p ), 1,
! 784: + WORK( N+1 ), 1 )
! 785: CALL DLASCL( 'G', 0, 0, AAPP,
! 786: + WORK( p ), M, 1,
! 787: + WORK( N+1 ), LDA, IERR )
! 788: AAPQ = DDOT( M, WORK( N+1 ), 1,
! 789: + A( 1, q ), 1 )*WORK( q ) / AAQQ
! 790: END IF
! 791: ELSE
! 792: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 793: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 794: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 795: + q ), 1 )*WORK( p )*WORK( q ) /
! 796: + AAQQ ) / AAPP
! 797: ELSE
! 798: CALL DCOPY( M, A( 1, q ), 1,
! 799: + WORK( N+1 ), 1 )
! 800: CALL DLASCL( 'G', 0, 0, AAQQ,
! 801: + WORK( q ), M, 1,
! 802: + WORK( N+1 ), LDA, IERR )
! 803: AAPQ = DDOT( M, WORK( N+1 ), 1,
! 804: + A( 1, p ), 1 )*WORK( p ) / AAPP
! 805: END IF
! 806: END IF
! 807: *
! 808: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
! 809: *
! 810: * TO rotate or NOT to rotate, THAT is the question ...
! 811: *
! 812: IF( DABS( AAPQ ).GT.TOL ) THEN
! 813: *
! 814: * .. rotate
! 815: *[RTD] ROTATED = ROTATED + ONE
! 816: *
! 817: IF( ir1.EQ.0 ) THEN
! 818: NOTROT = 0
! 819: PSKIPPED = 0
! 820: ISWROT = ISWROT + 1
! 821: END IF
! 822: *
! 823: IF( ROTOK ) THEN
! 824: *
! 825: AQOAP = AAQQ / AAPP
! 826: APOAQ = AAPP / AAQQ
! 827: THETA = -HALF*DABS( AQOAP-APOAQ ) /
! 828: + AAPQ
! 829: *
! 830: IF( DABS( THETA ).GT.BIGTHETA ) THEN
! 831: *
! 832: T = HALF / THETA
! 833: FASTR( 3 ) = T*WORK( p ) / WORK( q )
! 834: FASTR( 4 ) = -T*WORK( q ) /
! 835: + WORK( p )
! 836: CALL DROTM( M, A( 1, p ), 1,
! 837: + A( 1, q ), 1, FASTR )
! 838: IF( RSVEC )CALL DROTM( MVL,
! 839: + V( 1, p ), 1,
! 840: + V( 1, q ), 1,
! 841: + FASTR )
! 842: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 843: + ONE+T*APOAQ*AAPQ ) )
! 844: AAPP = AAPP*DSQRT( ONE-T*AQOAP*
! 845: + AAPQ )
! 846: MXSINJ = DMAX1( MXSINJ, DABS( T ) )
! 847: *
! 848: ELSE
! 849: *
! 850: * .. choose correct signum for THETA and rotate
! 851: *
! 852: THSIGN = -DSIGN( ONE, AAPQ )
! 853: T = ONE / ( THETA+THSIGN*
! 854: + DSQRT( ONE+THETA*THETA ) )
! 855: CS = DSQRT( ONE / ( ONE+T*T ) )
! 856: SN = T*CS
! 857: *
! 858: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
! 859: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 860: + ONE+T*APOAQ*AAPQ ) )
! 861: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 862: + ONE-T*AQOAP*AAPQ ) )
! 863: *
! 864: APOAQ = WORK( p ) / WORK( q )
! 865: AQOAP = WORK( q ) / WORK( p )
! 866: IF( WORK( p ).GE.ONE ) THEN
! 867: IF( WORK( q ).GE.ONE ) THEN
! 868: FASTR( 3 ) = T*APOAQ
! 869: FASTR( 4 ) = -T*AQOAP
! 870: WORK( p ) = WORK( p )*CS
! 871: WORK( q ) = WORK( q )*CS
! 872: CALL DROTM( M, A( 1, p ), 1,
! 873: + A( 1, q ), 1,
! 874: + FASTR )
! 875: IF( RSVEC )CALL DROTM( MVL,
! 876: + V( 1, p ), 1, V( 1, q ),
! 877: + 1, FASTR )
! 878: ELSE
! 879: CALL DAXPY( M, -T*AQOAP,
! 880: + A( 1, q ), 1,
! 881: + A( 1, p ), 1 )
! 882: CALL DAXPY( M, CS*SN*APOAQ,
! 883: + A( 1, p ), 1,
! 884: + A( 1, q ), 1 )
! 885: WORK( p ) = WORK( p )*CS
! 886: WORK( q ) = WORK( q ) / CS
! 887: IF( RSVEC ) THEN
! 888: CALL DAXPY( MVL, -T*AQOAP,
! 889: + V( 1, q ), 1,
! 890: + V( 1, p ), 1 )
! 891: CALL DAXPY( MVL,
! 892: + CS*SN*APOAQ,
! 893: + V( 1, p ), 1,
! 894: + V( 1, q ), 1 )
! 895: END IF
! 896: END IF
! 897: ELSE
! 898: IF( WORK( q ).GE.ONE ) THEN
! 899: CALL DAXPY( M, T*APOAQ,
! 900: + A( 1, p ), 1,
! 901: + A( 1, q ), 1 )
! 902: CALL DAXPY( M, -CS*SN*AQOAP,
! 903: + A( 1, q ), 1,
! 904: + A( 1, p ), 1 )
! 905: WORK( p ) = WORK( p ) / CS
! 906: WORK( q ) = WORK( q )*CS
! 907: IF( RSVEC ) THEN
! 908: CALL DAXPY( MVL, T*APOAQ,
! 909: + V( 1, p ), 1,
! 910: + V( 1, q ), 1 )
! 911: CALL DAXPY( MVL,
! 912: + -CS*SN*AQOAP,
! 913: + V( 1, q ), 1,
! 914: + V( 1, p ), 1 )
! 915: END IF
! 916: ELSE
! 917: IF( WORK( p ).GE.WORK( q ) )
! 918: + THEN
! 919: CALL DAXPY( M, -T*AQOAP,
! 920: + A( 1, q ), 1,
! 921: + A( 1, p ), 1 )
! 922: CALL DAXPY( M, CS*SN*APOAQ,
! 923: + A( 1, p ), 1,
! 924: + A( 1, q ), 1 )
! 925: WORK( p ) = WORK( p )*CS
! 926: WORK( q ) = WORK( q ) / CS
! 927: IF( RSVEC ) THEN
! 928: CALL DAXPY( MVL,
! 929: + -T*AQOAP,
! 930: + V( 1, q ), 1,
! 931: + V( 1, p ), 1 )
! 932: CALL DAXPY( MVL,
! 933: + CS*SN*APOAQ,
! 934: + V( 1, p ), 1,
! 935: + V( 1, q ), 1 )
! 936: END IF
! 937: ELSE
! 938: CALL DAXPY( M, T*APOAQ,
! 939: + A( 1, p ), 1,
! 940: + A( 1, q ), 1 )
! 941: CALL DAXPY( M,
! 942: + -CS*SN*AQOAP,
! 943: + A( 1, q ), 1,
! 944: + A( 1, p ), 1 )
! 945: WORK( p ) = WORK( p ) / CS
! 946: WORK( q ) = WORK( q )*CS
! 947: IF( RSVEC ) THEN
! 948: CALL DAXPY( MVL,
! 949: + T*APOAQ, V( 1, p ),
! 950: + 1, V( 1, q ), 1 )
! 951: CALL DAXPY( MVL,
! 952: + -CS*SN*AQOAP,
! 953: + V( 1, q ), 1,
! 954: + V( 1, p ), 1 )
! 955: END IF
! 956: END IF
! 957: END IF
! 958: END IF
! 959: END IF
! 960: *
! 961: ELSE
! 962: * .. have to use modified Gram-Schmidt like transformation
! 963: CALL DCOPY( M, A( 1, p ), 1,
! 964: + WORK( N+1 ), 1 )
! 965: CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
! 966: + 1, WORK( N+1 ), LDA,
! 967: + IERR )
! 968: CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
! 969: + 1, A( 1, q ), LDA, IERR )
! 970: TEMP1 = -AAPQ*WORK( p ) / WORK( q )
! 971: CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,
! 972: + A( 1, q ), 1 )
! 973: CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
! 974: + 1, A( 1, q ), LDA, IERR )
! 975: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 976: + ONE-AAPQ*AAPQ ) )
! 977: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 978: END IF
! 979: * END IF ROTOK THEN ... ELSE
! 980: *
! 981: * In the case of cancellation in updating SVA(q), SVA(p)
! 982: * recompute SVA(q), SVA(p).
! 983: *
! 984: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 985: + THEN
! 986: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 987: + ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 988: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
! 989: + WORK( q )
! 990: ELSE
! 991: T = ZERO
! 992: AAQQ = ZERO
! 993: CALL DLASSQ( M, A( 1, q ), 1, T,
! 994: + AAQQ )
! 995: SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
! 996: END IF
! 997: END IF
! 998: IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
! 999: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 1000: + ( AAPP.GT.ROOTSFMIN ) ) THEN
! 1001: AAPP = DNRM2( M, A( 1, p ), 1 )*
! 1002: + WORK( p )
! 1003: ELSE
! 1004: T = ZERO
! 1005: AAPP = ZERO
! 1006: CALL DLASSQ( M, A( 1, p ), 1, T,
! 1007: + AAPP )
! 1008: AAPP = T*DSQRT( AAPP )*WORK( p )
! 1009: END IF
! 1010: SVA( p ) = AAPP
! 1011: END IF
! 1012: *
! 1013: ELSE
! 1014: * A(:,p) and A(:,q) already numerically orthogonal
! 1015: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
! 1016: *[RTD] SKIPPED = SKIPPED + 1
! 1017: PSKIPPED = PSKIPPED + 1
! 1018: END IF
! 1019: ELSE
! 1020: * A(:,q) is zero column
! 1021: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
! 1022: PSKIPPED = PSKIPPED + 1
! 1023: END IF
! 1024: *
! 1025: IF( ( i.LE.SWBAND ) .AND.
! 1026: + ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 1027: IF( ir1.EQ.0 )AAPP = -AAPP
! 1028: NOTROT = 0
! 1029: GO TO 2103
! 1030: END IF
! 1031: *
! 1032: 2002 CONTINUE
! 1033: * END q-LOOP
! 1034: *
! 1035: 2103 CONTINUE
! 1036: * bailed out of q-loop
! 1037: *
! 1038: SVA( p ) = AAPP
! 1039: *
! 1040: ELSE
! 1041: SVA( p ) = AAPP
! 1042: IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
! 1043: + NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
! 1044: END IF
! 1045: *
! 1046: 2001 CONTINUE
! 1047: * end of the p-loop
! 1048: * end of doing the block ( ibr, ibr )
! 1049: 1002 CONTINUE
! 1050: * end of ir1-loop
! 1051: *
! 1052: * ... go to the off diagonal blocks
! 1053: *
! 1054: igl = ( ibr-1 )*KBL + 1
! 1055: *
! 1056: DO 2010 jbc = ibr + 1, NBL
! 1057: *
! 1058: jgl = ( jbc-1 )*KBL + 1
! 1059: *
! 1060: * doing the block at ( ibr, jbc )
! 1061: *
! 1062: IJBLSK = 0
! 1063: DO 2100 p = igl, MIN0( igl+KBL-1, N )
! 1064: *
! 1065: AAPP = SVA( p )
! 1066: IF( AAPP.GT.ZERO ) THEN
! 1067: *
! 1068: PSKIPPED = 0
! 1069: *
! 1070: DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
! 1071: *
! 1072: AAQQ = SVA( q )
! 1073: IF( AAQQ.GT.ZERO ) THEN
! 1074: AAPP0 = AAPP
! 1075: *
! 1076: * .. M x 2 Jacobi SVD ..
! 1077: *
! 1078: * Safe Gram matrix computation
! 1079: *
! 1080: IF( AAQQ.GE.ONE ) THEN
! 1081: IF( AAPP.GE.AAQQ ) THEN
! 1082: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 1083: ELSE
! 1084: ROTOK = ( SMALL*AAQQ ).LE.AAPP
! 1085: END IF
! 1086: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
! 1087: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 1088: + q ), 1 )*WORK( p )*WORK( q ) /
! 1089: + AAQQ ) / AAPP
! 1090: ELSE
! 1091: CALL DCOPY( M, A( 1, p ), 1,
! 1092: + WORK( N+1 ), 1 )
! 1093: CALL DLASCL( 'G', 0, 0, AAPP,
! 1094: + WORK( p ), M, 1,
! 1095: + WORK( N+1 ), LDA, IERR )
! 1096: AAPQ = DDOT( M, WORK( N+1 ), 1,
! 1097: + A( 1, q ), 1 )*WORK( q ) / AAQQ
! 1098: END IF
! 1099: ELSE
! 1100: IF( AAPP.GE.AAQQ ) THEN
! 1101: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 1102: ELSE
! 1103: ROTOK = AAQQ.LE.( AAPP / SMALL )
! 1104: END IF
! 1105: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 1106: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
! 1107: + q ), 1 )*WORK( p )*WORK( q ) /
! 1108: + AAQQ ) / AAPP
! 1109: ELSE
! 1110: CALL DCOPY( M, A( 1, q ), 1,
! 1111: + WORK( N+1 ), 1 )
! 1112: CALL DLASCL( 'G', 0, 0, AAQQ,
! 1113: + WORK( q ), M, 1,
! 1114: + WORK( N+1 ), LDA, IERR )
! 1115: AAPQ = DDOT( M, WORK( N+1 ), 1,
! 1116: + A( 1, p ), 1 )*WORK( p ) / AAPP
! 1117: END IF
! 1118: END IF
! 1119: *
! 1120: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
! 1121: *
! 1122: * TO rotate or NOT to rotate, THAT is the question ...
! 1123: *
! 1124: IF( DABS( AAPQ ).GT.TOL ) THEN
! 1125: NOTROT = 0
! 1126: *[RTD] ROTATED = ROTATED + 1
! 1127: PSKIPPED = 0
! 1128: ISWROT = ISWROT + 1
! 1129: *
! 1130: IF( ROTOK ) THEN
! 1131: *
! 1132: AQOAP = AAQQ / AAPP
! 1133: APOAQ = AAPP / AAQQ
! 1134: THETA = -HALF*DABS( AQOAP-APOAQ ) /
! 1135: + AAPQ
! 1136: IF( AAQQ.GT.AAPP0 )THETA = -THETA
! 1137: *
! 1138: IF( DABS( THETA ).GT.BIGTHETA ) THEN
! 1139: T = HALF / THETA
! 1140: FASTR( 3 ) = T*WORK( p ) / WORK( q )
! 1141: FASTR( 4 ) = -T*WORK( q ) /
! 1142: + WORK( p )
! 1143: CALL DROTM( M, A( 1, p ), 1,
! 1144: + A( 1, q ), 1, FASTR )
! 1145: IF( RSVEC )CALL DROTM( MVL,
! 1146: + V( 1, p ), 1,
! 1147: + V( 1, q ), 1,
! 1148: + FASTR )
! 1149: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 1150: + ONE+T*APOAQ*AAPQ ) )
! 1151: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 1152: + ONE-T*AQOAP*AAPQ ) )
! 1153: MXSINJ = DMAX1( MXSINJ, DABS( T ) )
! 1154: ELSE
! 1155: *
! 1156: * .. choose correct signum for THETA and rotate
! 1157: *
! 1158: THSIGN = -DSIGN( ONE, AAPQ )
! 1159: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
! 1160: T = ONE / ( THETA+THSIGN*
! 1161: + DSQRT( ONE+THETA*THETA ) )
! 1162: CS = DSQRT( ONE / ( ONE+T*T ) )
! 1163: SN = T*CS
! 1164: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
! 1165: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 1166: + ONE+T*APOAQ*AAPQ ) )
! 1167: AAPP = AAPP*DSQRT( ONE-T*AQOAP*
! 1168: + AAPQ )
! 1169: *
! 1170: APOAQ = WORK( p ) / WORK( q )
! 1171: AQOAP = WORK( q ) / WORK( p )
! 1172: IF( WORK( p ).GE.ONE ) THEN
! 1173: *
! 1174: IF( WORK( q ).GE.ONE ) THEN
! 1175: FASTR( 3 ) = T*APOAQ
! 1176: FASTR( 4 ) = -T*AQOAP
! 1177: WORK( p ) = WORK( p )*CS
! 1178: WORK( q ) = WORK( q )*CS
! 1179: CALL DROTM( M, A( 1, p ), 1,
! 1180: + A( 1, q ), 1,
! 1181: + FASTR )
! 1182: IF( RSVEC )CALL DROTM( MVL,
! 1183: + V( 1, p ), 1, V( 1, q ),
! 1184: + 1, FASTR )
! 1185: ELSE
! 1186: CALL DAXPY( M, -T*AQOAP,
! 1187: + A( 1, q ), 1,
! 1188: + A( 1, p ), 1 )
! 1189: CALL DAXPY( M, CS*SN*APOAQ,
! 1190: + A( 1, p ), 1,
! 1191: + A( 1, q ), 1 )
! 1192: IF( RSVEC ) THEN
! 1193: CALL DAXPY( MVL, -T*AQOAP,
! 1194: + V( 1, q ), 1,
! 1195: + V( 1, p ), 1 )
! 1196: CALL DAXPY( MVL,
! 1197: + CS*SN*APOAQ,
! 1198: + V( 1, p ), 1,
! 1199: + V( 1, q ), 1 )
! 1200: END IF
! 1201: WORK( p ) = WORK( p )*CS
! 1202: WORK( q ) = WORK( q ) / CS
! 1203: END IF
! 1204: ELSE
! 1205: IF( WORK( q ).GE.ONE ) THEN
! 1206: CALL DAXPY( M, T*APOAQ,
! 1207: + A( 1, p ), 1,
! 1208: + A( 1, q ), 1 )
! 1209: CALL DAXPY( M, -CS*SN*AQOAP,
! 1210: + A( 1, q ), 1,
! 1211: + A( 1, p ), 1 )
! 1212: IF( RSVEC ) THEN
! 1213: CALL DAXPY( MVL, T*APOAQ,
! 1214: + V( 1, p ), 1,
! 1215: + V( 1, q ), 1 )
! 1216: CALL DAXPY( MVL,
! 1217: + -CS*SN*AQOAP,
! 1218: + V( 1, q ), 1,
! 1219: + V( 1, p ), 1 )
! 1220: END IF
! 1221: WORK( p ) = WORK( p ) / CS
! 1222: WORK( q ) = WORK( q )*CS
! 1223: ELSE
! 1224: IF( WORK( p ).GE.WORK( q ) )
! 1225: + THEN
! 1226: CALL DAXPY( M, -T*AQOAP,
! 1227: + A( 1, q ), 1,
! 1228: + A( 1, p ), 1 )
! 1229: CALL DAXPY( M, CS*SN*APOAQ,
! 1230: + A( 1, p ), 1,
! 1231: + A( 1, q ), 1 )
! 1232: WORK( p ) = WORK( p )*CS
! 1233: WORK( q ) = WORK( q ) / CS
! 1234: IF( RSVEC ) THEN
! 1235: CALL DAXPY( MVL,
! 1236: + -T*AQOAP,
! 1237: + V( 1, q ), 1,
! 1238: + V( 1, p ), 1 )
! 1239: CALL DAXPY( MVL,
! 1240: + CS*SN*APOAQ,
! 1241: + V( 1, p ), 1,
! 1242: + V( 1, q ), 1 )
! 1243: END IF
! 1244: ELSE
! 1245: CALL DAXPY( M, T*APOAQ,
! 1246: + A( 1, p ), 1,
! 1247: + A( 1, q ), 1 )
! 1248: CALL DAXPY( M,
! 1249: + -CS*SN*AQOAP,
! 1250: + A( 1, q ), 1,
! 1251: + A( 1, p ), 1 )
! 1252: WORK( p ) = WORK( p ) / CS
! 1253: WORK( q ) = WORK( q )*CS
! 1254: IF( RSVEC ) THEN
! 1255: CALL DAXPY( MVL,
! 1256: + T*APOAQ, V( 1, p ),
! 1257: + 1, V( 1, q ), 1 )
! 1258: CALL DAXPY( MVL,
! 1259: + -CS*SN*AQOAP,
! 1260: + V( 1, q ), 1,
! 1261: + V( 1, p ), 1 )
! 1262: END IF
! 1263: END IF
! 1264: END IF
! 1265: END IF
! 1266: END IF
! 1267: *
! 1268: ELSE
! 1269: IF( AAPP.GT.AAQQ ) THEN
! 1270: CALL DCOPY( M, A( 1, p ), 1,
! 1271: + WORK( N+1 ), 1 )
! 1272: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
! 1273: + M, 1, WORK( N+1 ), LDA,
! 1274: + IERR )
! 1275: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
! 1276: + M, 1, A( 1, q ), LDA,
! 1277: + IERR )
! 1278: TEMP1 = -AAPQ*WORK( p ) / WORK( q )
! 1279: CALL DAXPY( M, TEMP1, WORK( N+1 ),
! 1280: + 1, A( 1, q ), 1 )
! 1281: CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
! 1282: + M, 1, A( 1, q ), LDA,
! 1283: + IERR )
! 1284: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 1285: + ONE-AAPQ*AAPQ ) )
! 1286: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 1287: ELSE
! 1288: CALL DCOPY( M, A( 1, q ), 1,
! 1289: + WORK( N+1 ), 1 )
! 1290: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
! 1291: + M, 1, WORK( N+1 ), LDA,
! 1292: + IERR )
! 1293: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
! 1294: + M, 1, A( 1, p ), LDA,
! 1295: + IERR )
! 1296: TEMP1 = -AAPQ*WORK( q ) / WORK( p )
! 1297: CALL DAXPY( M, TEMP1, WORK( N+1 ),
! 1298: + 1, A( 1, p ), 1 )
! 1299: CALL DLASCL( 'G', 0, 0, ONE, AAPP,
! 1300: + M, 1, A( 1, p ), LDA,
! 1301: + IERR )
! 1302: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
! 1303: + ONE-AAPQ*AAPQ ) )
! 1304: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 1305: END IF
! 1306: END IF
! 1307: * END IF ROTOK THEN ... ELSE
! 1308: *
! 1309: * In the case of cancellation in updating SVA(q)
! 1310: * .. recompute SVA(q)
! 1311: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 1312: + THEN
! 1313: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 1314: + ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 1315: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
! 1316: + WORK( q )
! 1317: ELSE
! 1318: T = ZERO
! 1319: AAQQ = ZERO
! 1320: CALL DLASSQ( M, A( 1, q ), 1, T,
! 1321: + AAQQ )
! 1322: SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
! 1323: END IF
! 1324: END IF
! 1325: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
! 1326: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 1327: + ( AAPP.GT.ROOTSFMIN ) ) THEN
! 1328: AAPP = DNRM2( M, A( 1, p ), 1 )*
! 1329: + WORK( p )
! 1330: ELSE
! 1331: T = ZERO
! 1332: AAPP = ZERO
! 1333: CALL DLASSQ( M, A( 1, p ), 1, T,
! 1334: + AAPP )
! 1335: AAPP = T*DSQRT( AAPP )*WORK( p )
! 1336: END IF
! 1337: SVA( p ) = AAPP
! 1338: END IF
! 1339: * end of OK rotation
! 1340: ELSE
! 1341: NOTROT = NOTROT + 1
! 1342: *[RTD] SKIPPED = SKIPPED + 1
! 1343: PSKIPPED = PSKIPPED + 1
! 1344: IJBLSK = IJBLSK + 1
! 1345: END IF
! 1346: ELSE
! 1347: NOTROT = NOTROT + 1
! 1348: PSKIPPED = PSKIPPED + 1
! 1349: IJBLSK = IJBLSK + 1
! 1350: END IF
! 1351: *
! 1352: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
! 1353: + THEN
! 1354: SVA( p ) = AAPP
! 1355: NOTROT = 0
! 1356: GO TO 2011
! 1357: END IF
! 1358: IF( ( i.LE.SWBAND ) .AND.
! 1359: + ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 1360: AAPP = -AAPP
! 1361: NOTROT = 0
! 1362: GO TO 2203
! 1363: END IF
! 1364: *
! 1365: 2200 CONTINUE
! 1366: * end of the q-loop
! 1367: 2203 CONTINUE
! 1368: *
! 1369: SVA( p ) = AAPP
! 1370: *
! 1371: ELSE
! 1372: *
! 1373: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
! 1374: + MIN0( jgl+KBL-1, N ) - jgl + 1
! 1375: IF( AAPP.LT.ZERO )NOTROT = 0
! 1376: *
! 1377: END IF
! 1378: *
! 1379: 2100 CONTINUE
! 1380: * end of the p-loop
! 1381: 2010 CONTINUE
! 1382: * end of the jbc-loop
! 1383: 2011 CONTINUE
! 1384: *2011 bailed out of the jbc-loop
! 1385: DO 2012 p = igl, MIN0( igl+KBL-1, N )
! 1386: SVA( p ) = DABS( SVA( p ) )
! 1387: 2012 CONTINUE
! 1388: ***
! 1389: 2000 CONTINUE
! 1390: *2000 :: end of the ibr-loop
! 1391: *
! 1392: * .. update SVA(N)
! 1393: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
! 1394: + THEN
! 1395: SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
! 1396: ELSE
! 1397: T = ZERO
! 1398: AAPP = ZERO
! 1399: CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
! 1400: SVA( N ) = T*DSQRT( AAPP )*WORK( N )
! 1401: END IF
! 1402: *
! 1403: * Additional steering devices
! 1404: *
! 1405: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
! 1406: + ( ISWROT.LE.N ) ) )SWBAND = i
! 1407: *
! 1408: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
! 1409: + TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
! 1410: GO TO 1994
! 1411: END IF
! 1412: *
! 1413: IF( NOTROT.GE.EMPTSW )GO TO 1994
! 1414: *
! 1415: 1993 CONTINUE
! 1416: * end i=1:NSWEEP loop
! 1417: *
! 1418: * #:( Reaching this point means that the procedure has not converged.
! 1419: INFO = NSWEEP - 1
! 1420: GO TO 1995
! 1421: *
! 1422: 1994 CONTINUE
! 1423: * #:) Reaching this point means numerical convergence after the i-th
! 1424: * sweep.
! 1425: *
! 1426: INFO = 0
! 1427: * #:) INFO = 0 confirms successful iterations.
! 1428: 1995 CONTINUE
! 1429: *
! 1430: * Sort the singular values and find how many are above
! 1431: * the underflow threshold.
! 1432: *
! 1433: N2 = 0
! 1434: N4 = 0
! 1435: DO 5991 p = 1, N - 1
! 1436: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 1437: IF( p.NE.q ) THEN
! 1438: TEMP1 = SVA( p )
! 1439: SVA( p ) = SVA( q )
! 1440: SVA( q ) = TEMP1
! 1441: TEMP1 = WORK( p )
! 1442: WORK( p ) = WORK( q )
! 1443: WORK( q ) = TEMP1
! 1444: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 1445: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
! 1446: END IF
! 1447: IF( SVA( p ).NE.ZERO ) THEN
! 1448: N4 = N4 + 1
! 1449: IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
! 1450: END IF
! 1451: 5991 CONTINUE
! 1452: IF( SVA( N ).NE.ZERO ) THEN
! 1453: N4 = N4 + 1
! 1454: IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
! 1455: END IF
! 1456: *
! 1457: * Normalize the left singular vectors.
! 1458: *
! 1459: IF( LSVEC .OR. UCTOL ) THEN
! 1460: DO 1998 p = 1, N2
! 1461: CALL DSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
! 1462: 1998 CONTINUE
! 1463: END IF
! 1464: *
! 1465: * Scale the product of Jacobi rotations (assemble the fast rotations).
! 1466: *
! 1467: IF( RSVEC ) THEN
! 1468: IF( APPLV ) THEN
! 1469: DO 2398 p = 1, N
! 1470: CALL DSCAL( MVL, WORK( p ), V( 1, p ), 1 )
! 1471: 2398 CONTINUE
! 1472: ELSE
! 1473: DO 2399 p = 1, N
! 1474: TEMP1 = ONE / DNRM2( MVL, V( 1, p ), 1 )
! 1475: CALL DSCAL( MVL, TEMP1, V( 1, p ), 1 )
! 1476: 2399 CONTINUE
! 1477: END IF
! 1478: END IF
! 1479: *
! 1480: * Undo scaling, if necessary (and possible).
! 1481: IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
! 1482: + SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
! 1483: + ( SFMIN / SCALE ) ) ) ) THEN
! 1484: DO 2400 p = 1, N
! 1485: SVA( p ) = SCALE*SVA( p )
! 1486: 2400 CONTINUE
! 1487: SCALE = ONE
! 1488: END IF
! 1489: *
! 1490: WORK( 1 ) = SCALE
! 1491: * The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
! 1492: * then some of the singular values may overflow or underflow and
! 1493: * the spectrum is given in this factored representation.
! 1494: *
! 1495: WORK( 2 ) = DBLE( N4 )
! 1496: * N4 is the number of computed nonzero singular values of A.
! 1497: *
! 1498: WORK( 3 ) = DBLE( N2 )
! 1499: * N2 is the number of singular values of A greater than SFMIN.
! 1500: * If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
! 1501: * that may carry some information.
! 1502: *
! 1503: WORK( 4 ) = DBLE( i )
! 1504: * i is the index of the last sweep before declaring convergence.
! 1505: *
! 1506: WORK( 5 ) = MXAAPQ
! 1507: * MXAAPQ is the largest absolute value of scaled pivots in the
! 1508: * last sweep
! 1509: *
! 1510: WORK( 6 ) = MXSINJ
! 1511: * MXSINJ is the largest absolute value of the sines of Jacobi angles
! 1512: * in the last sweep
! 1513: *
! 1514: RETURN
! 1515: * ..
! 1516: * .. END OF DGESVJ
! 1517: * ..
! 1518: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgesvj.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack