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