Annotation of rpl/lapack/lapack/zgsvj0.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZGSVJ0 pre-processor for the routine dgesvj.
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZGSVJ0 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgsvj0.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgsvj0.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgsvj0.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
! 22: * SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
! 26: * DOUBLE PRECISION EPS, SFMIN, TOL
! 27: * CHARACTER*1 JOBV
! 28: * ..
! 29: * .. Array Arguments ..
! 30: * COMPLEX*16 A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
! 31: * DOUBLE PRECISION SVA( N )
! 32: * ..
! 33: *
! 34: *
! 35: *> \par Purpose:
! 36: * =============
! 37: *>
! 38: *> \verbatim
! 39: *>
! 40: *> ZGSVJ0 is called from ZGESVJ as a pre-processor and that is its main
! 41: *> purpose. It applies Jacobi rotations in the same way as ZGESVJ does, but
! 42: *> it does not check convergence (stopping criterion). Few tuning
! 43: *> parameters (marked by [TP]) are available for the implementer.
! 44: *> \endverbatim
! 45: *
! 46: * Arguments:
! 47: * ==========
! 48: *
! 49: *> \param[in] JOBV
! 50: *> \verbatim
! 51: *> JOBV is CHARACTER*1
! 52: *> Specifies whether the output from this procedure is used
! 53: *> to compute the matrix V:
! 54: *> = 'V': the product of the Jacobi rotations is accumulated
! 55: *> by postmulyiplying the N-by-N array V.
! 56: *> (See the description of V.)
! 57: *> = 'A': the product of the Jacobi rotations is accumulated
! 58: *> by postmulyiplying the MV-by-N array V.
! 59: *> (See the descriptions of MV and V.)
! 60: *> = 'N': the Jacobi rotations are not accumulated.
! 61: *> \endverbatim
! 62: *>
! 63: *> \param[in] M
! 64: *> \verbatim
! 65: *> M is INTEGER
! 66: *> The number of rows of the input matrix A. M >= 0.
! 67: *> \endverbatim
! 68: *>
! 69: *> \param[in] N
! 70: *> \verbatim
! 71: *> N is INTEGER
! 72: *> The number of columns of the input matrix A.
! 73: *> M >= N >= 0.
! 74: *> \endverbatim
! 75: *>
! 76: *> \param[in,out] A
! 77: *> \verbatim
! 78: *> A is COMPLEX*16 array, dimension (LDA,N)
! 79: *> On entry, M-by-N matrix A, such that A*diag(D) represents
! 80: *> the input matrix.
! 81: *> On exit,
! 82: *> A_onexit * diag(D_onexit) represents the input matrix A*diag(D)
! 83: *> post-multiplied by a sequence of Jacobi rotations, where the
! 84: *> rotation threshold and the total number of sweeps are given in
! 85: *> TOL and NSWEEP, respectively.
! 86: *> (See the descriptions of D, TOL and NSWEEP.)
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] LDA
! 90: *> \verbatim
! 91: *> LDA is INTEGER
! 92: *> The leading dimension of the array A. LDA >= max(1,M).
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[in,out] D
! 96: *> \verbatim
! 97: *> D is COMPLEX*16 array, dimension (N)
! 98: *> The array D accumulates the scaling factors from the complex scaled
! 99: *> Jacobi rotations.
! 100: *> On entry, A*diag(D) represents the input matrix.
! 101: *> On exit, A_onexit*diag(D_onexit) represents the input matrix
! 102: *> post-multiplied by a sequence of Jacobi rotations, where the
! 103: *> rotation threshold and the total number of sweeps are given in
! 104: *> TOL and NSWEEP, respectively.
! 105: *> (See the descriptions of A, TOL and NSWEEP.)
! 106: *> \endverbatim
! 107: *>
! 108: *> \param[in,out] SVA
! 109: *> \verbatim
! 110: *> SVA is DOUBLE PRECISION array, dimension (N)
! 111: *> On entry, SVA contains the Euclidean norms of the columns of
! 112: *> the matrix A*diag(D).
! 113: *> On exit, SVA contains the Euclidean norms of the columns of
! 114: *> the matrix A_onexit*diag(D_onexit).
! 115: *>
! 116: *> \param[in] MV
! 117: *> \verbatim
! 118: *> MV is INTEGER
! 119: *> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
! 120: *> sequence of Jacobi rotations.
! 121: *> If JOBV = 'N', then MV is not referenced.
! 122: *> \endverbatim
! 123: *>
! 124: *> \param[in,out] V
! 125: *> \verbatim
! 126: *> V is COMPLEX*16 array, dimension (LDV,N)
! 127: *> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
! 128: *> sequence of Jacobi rotations.
! 129: *> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
! 130: *> sequence of Jacobi rotations.
! 131: *> If JOBV = 'N', then V is not referenced.
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[in] LDV
! 135: *> \verbatim
! 136: *> LDV is INTEGER
! 137: *> The leading dimension of the array V, LDV >= 1.
! 138: *> If JOBV = 'V', LDV .GE. N.
! 139: *> If JOBV = 'A', LDV .GE. MV.
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[in] EPS
! 143: *> \verbatim
! 144: *> EPS is DOUBLE PRECISION
! 145: *> EPS = DLAMCH('Epsilon')
! 146: *> \endverbatim
! 147: *>
! 148: *> \param[in] SFMIN
! 149: *> \verbatim
! 150: *> SFMIN is DOUBLE PRECISION
! 151: *> SFMIN = DLAMCH('Safe Minimum')
! 152: *> \endverbatim
! 153: *>
! 154: *> \param[in] TOL
! 155: *> \verbatim
! 156: *> TOL is DOUBLE PRECISION
! 157: *> TOL is the threshold for Jacobi rotations. For a pair
! 158: *> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
! 159: *> applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
! 160: *> \endverbatim
! 161: *>
! 162: *> \param[in] NSWEEP
! 163: *> \verbatim
! 164: *> NSWEEP is INTEGER
! 165: *> NSWEEP is the number of sweeps of Jacobi rotations to be
! 166: *> performed.
! 167: *> \endverbatim
! 168: *>
! 169: *> \param[out] WORK
! 170: *> \verbatim
! 171: *> WORK is COMPLEX*16 array, dimension LWORK.
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[in] LWORK
! 175: *> \verbatim
! 176: *> LWORK is INTEGER
! 177: *> LWORK is the dimension of WORK. LWORK .GE. M.
! 178: *> \endverbatim
! 179: *>
! 180: *> \param[out] INFO
! 181: *> \verbatim
! 182: *> INFO is INTEGER
! 183: *> = 0 : successful exit.
! 184: *> < 0 : if INFO = -i, then the i-th argument had an illegal value
! 185: *> \endverbatim
! 186: *
! 187: * Authors:
! 188: * ========
! 189: *
! 190: *> \author Univ. of Tennessee
! 191: *> \author Univ. of California Berkeley
! 192: *> \author Univ. of Colorado Denver
! 193: *> \author NAG Ltd.
! 194: *
! 195: *> \date November 2015
! 196: *
! 197: *> \ingroup complex16OTHERcomputational
! 198: *>
! 199: *> \par Further Details:
! 200: * =====================
! 201: *>
! 202: *> ZGSVJ0 is used just to enable ZGESVJ to call a simplified version of
! 203: *> itself to work on a submatrix of the original matrix.
! 204: *>
! 205: *> Contributors:
! 206: * =============
! 207: *>
! 208: *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
! 209: *>
! 210: *> Bugs, Examples and Comments:
! 211: * ============================
! 212: *>
! 213: *> Please report all bugs and send interesting test examples and comments to
! 214: *> drmac@math.hr. Thank you.
! 215: *
! 216: * =====================================================================
! 217: SUBROUTINE ZGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
! 218: $ SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
! 219: *
! 220: * -- LAPACK computational routine (version 3.6.0) --
! 221: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 222: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 223: * November 2015
! 224: *
! 225: IMPLICIT NONE
! 226: * .. Scalar Arguments ..
! 227: INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
! 228: DOUBLE PRECISION EPS, SFMIN, TOL
! 229: CHARACTER*1 JOBV
! 230: * ..
! 231: * .. Array Arguments ..
! 232: COMPLEX*16 A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
! 233: DOUBLE PRECISION SVA( N )
! 234: * ..
! 235: *
! 236: * =====================================================================
! 237: *
! 238: * .. Local Parameters ..
! 239: DOUBLE PRECISION ZERO, HALF, ONE
! 240: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0)
! 241: COMPLEX*16 CZERO, CONE
! 242: PARAMETER ( CZERO = (0.0D0, 0.0D0), CONE = (1.0D0, 0.0D0) )
! 243: * ..
! 244: * .. Local Scalars ..
! 245: COMPLEX*16 AAPQ, OMPQ
! 246: DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
! 247: $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
! 248: $ ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
! 249: $ THSIGN
! 250: INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
! 251: $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
! 252: $ NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
! 253: LOGICAL APPLV, ROTOK, RSVEC
! 254: * ..
! 255: * ..
! 256: * .. Intrinsic Functions ..
! 257: INTRINSIC ABS, DMAX1, DCONJG, DFLOAT, MIN0, DSIGN, DSQRT
! 258: * ..
! 259: * .. External Functions ..
! 260: DOUBLE PRECISION DZNRM2
! 261: COMPLEX*16 ZDOTC
! 262: INTEGER IDAMAX
! 263: LOGICAL LSAME
! 264: EXTERNAL IDAMAX, LSAME, ZDOTC, DZNRM2
! 265: * ..
! 266: * ..
! 267: * .. External Subroutines ..
! 268: * ..
! 269: * from BLAS
! 270: EXTERNAL ZCOPY, ZROT, ZSWAP
! 271: * from LAPACK
! 272: EXTERNAL ZLASCL, ZLASSQ, XERBLA
! 273: * ..
! 274: * .. Executable Statements ..
! 275: *
! 276: * Test the input parameters.
! 277: *
! 278: APPLV = LSAME( JOBV, 'A' )
! 279: RSVEC = LSAME( JOBV, 'V' )
! 280: IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
! 281: INFO = -1
! 282: ELSE IF( M.LT.0 ) THEN
! 283: INFO = -2
! 284: ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
! 285: INFO = -3
! 286: ELSE IF( LDA.LT.M ) THEN
! 287: INFO = -5
! 288: ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
! 289: INFO = -8
! 290: ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
! 291: $ ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
! 292: INFO = -10
! 293: ELSE IF( TOL.LE.EPS ) THEN
! 294: INFO = -13
! 295: ELSE IF( NSWEEP.LT.0 ) THEN
! 296: INFO = -14
! 297: ELSE IF( LWORK.LT.M ) THEN
! 298: INFO = -16
! 299: ELSE
! 300: INFO = 0
! 301: END IF
! 302: *
! 303: * #:(
! 304: IF( INFO.NE.0 ) THEN
! 305: CALL XERBLA( 'ZGSVJ0', -INFO )
! 306: RETURN
! 307: END IF
! 308: *
! 309: IF( RSVEC ) THEN
! 310: MVL = N
! 311: ELSE IF( APPLV ) THEN
! 312: MVL = MV
! 313: END IF
! 314: RSVEC = RSVEC .OR. APPLV
! 315:
! 316: ROOTEPS = DSQRT( EPS )
! 317: ROOTSFMIN = DSQRT( SFMIN )
! 318: SMALL = SFMIN / EPS
! 319: BIG = ONE / SFMIN
! 320: ROOTBIG = ONE / ROOTSFMIN
! 321: BIGTHETA = ONE / ROOTEPS
! 322: ROOTTOL = DSQRT( TOL )
! 323: *
! 324: * .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
! 325: *
! 326: EMPTSW = ( N*( N-1 ) ) / 2
! 327: NOTROT = 0
! 328: *
! 329: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
! 330: *
! 331:
! 332: SWBAND = 0
! 333: *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
! 334: * if ZGESVJ is used as a computational routine in the preconditioned
! 335: * Jacobi SVD algorithm ZGEJSV. For sweeps i=1:SWBAND the procedure
! 336: * works on pivots inside a band-like region around the diagonal.
! 337: * The boundaries are determined dynamically, based on the number of
! 338: * pivots above a threshold.
! 339: *
! 340: KBL = MIN0( 8, N )
! 341: *[TP] KBL is a tuning parameter that defines the tile size in the
! 342: * tiling of the p-q loops of pivot pairs. In general, an optimal
! 343: * value of KBL depends on the matrix dimensions and on the
! 344: * parameters of the computer's memory.
! 345: *
! 346: NBL = N / KBL
! 347: IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
! 348: *
! 349: BLSKIP = KBL**2
! 350: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
! 351: *
! 352: ROWSKIP = MIN0( 5, KBL )
! 353: *[TP] ROWSKIP is a tuning parameter.
! 354: *
! 355: LKAHEAD = 1
! 356: *[TP] LKAHEAD is a tuning parameter.
! 357: *
! 358: * Quasi block transformations, using the lower (upper) triangular
! 359: * structure of the input matrix. The quasi-block-cycling usually
! 360: * invokes cubic convergence. Big part of this cycle is done inside
! 361: * canonical subspaces of dimensions less than M.
! 362: *
! 363: *
! 364: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
! 365: *
! 366: DO 1993 i = 1, NSWEEP
! 367: *
! 368: * .. go go go ...
! 369: *
! 370: MXAAPQ = ZERO
! 371: MXSINJ = ZERO
! 372: ISWROT = 0
! 373: *
! 374: NOTROT = 0
! 375: PSKIPPED = 0
! 376: *
! 377: * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
! 378: * 1 <= p < q <= N. This is the first step toward a blocked implementation
! 379: * of the rotations. New implementation, based on block transformations,
! 380: * is under development.
! 381: *
! 382: DO 2000 ibr = 1, NBL
! 383: *
! 384: igl = ( ibr-1 )*KBL + 1
! 385: *
! 386: DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
! 387: *
! 388: igl = igl + ir1*KBL
! 389: *
! 390: DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
! 391: *
! 392: * .. de Rijk's pivoting
! 393: *
! 394: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 395: IF( p.NE.q ) THEN
! 396: CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 397: IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1,
! 398: $ V( 1, q ), 1 )
! 399: TEMP1 = SVA( p )
! 400: SVA( p ) = SVA( q )
! 401: SVA( q ) = TEMP1
! 402: AAPQ = D(p)
! 403: D(p) = D(q)
! 404: D(q) = AAPQ
! 405: END IF
! 406: *
! 407: IF( ir1.EQ.0 ) THEN
! 408: *
! 409: * Column norms are periodically updated by explicit
! 410: * norm computation.
! 411: * Caveat:
! 412: * Unfortunately, some BLAS implementations compute SNCRM2(M,A(1,p),1)
! 413: * as SQRT(S=ZDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to
! 414: * overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
! 415: * underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
! 416: * Hence, DZNRM2 cannot be trusted, not even in the case when
! 417: * the true norm is far from the under(over)flow boundaries.
! 418: * If properly implemented DZNRM2 is available, the IF-THEN-ELSE-END IF
! 419: * below should be replaced with "AAPP = DZNRM2( M, A(1,p), 1 )".
! 420: *
! 421: IF( ( SVA( p ).LT.ROOTBIG ) .AND.
! 422: $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN
! 423: SVA( p ) = DZNRM2( M, A( 1, p ), 1 )
! 424: ELSE
! 425: TEMP1 = ZERO
! 426: AAPP = ONE
! 427: CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
! 428: SVA( p ) = TEMP1*DSQRT( AAPP )
! 429: END IF
! 430: AAPP = SVA( p )
! 431: ELSE
! 432: AAPP = SVA( p )
! 433: END IF
! 434: *
! 435: IF( AAPP.GT.ZERO ) THEN
! 436: *
! 437: PSKIPPED = 0
! 438: *
! 439: DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
! 440: *
! 441: AAQQ = SVA( q )
! 442: *
! 443: IF( AAQQ.GT.ZERO ) THEN
! 444: *
! 445: AAPP0 = AAPP
! 446: IF( AAQQ.GE.ONE ) THEN
! 447: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 448: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
! 449: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
! 450: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
! 451: ELSE
! 452: CALL ZCOPY( M, A( 1, p ), 1,
! 453: $ WORK, 1 )
! 454: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
! 455: $ M, 1, WORK, LDA, IERR )
! 456: AAPQ = ZDOTC( M, WORK, 1,
! 457: $ A( 1, q ), 1 ) / AAQQ
! 458: END IF
! 459: ELSE
! 460: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 461: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 462: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
! 463: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
! 464: ELSE
! 465: CALL ZCOPY( M, A( 1, q ), 1,
! 466: $ WORK, 1 )
! 467: CALL ZLASCL( 'G', 0, 0, AAQQ,
! 468: $ ONE, M, 1,
! 469: $ WORK, LDA, IERR )
! 470: AAPQ = ZDOTC( M, A( 1, p ), 1,
! 471: $ WORK, 1 ) / AAPP
! 472: END IF
! 473: END IF
! 474: *
! 475: OMPQ = AAPQ / ABS(AAPQ)
! 476: * AAPQ = AAPQ * DCONJG( CWORK(p) ) * CWORK(q)
! 477: AAPQ1 = -ABS(AAPQ)
! 478: MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
! 479: *
! 480: * TO rotate or NOT to rotate, THAT is the question ...
! 481: *
! 482: IF( ABS( AAPQ1 ).GT.TOL ) THEN
! 483: *
! 484: * .. rotate
! 485: *[RTD] ROTATED = ROTATED + ONE
! 486: *
! 487: IF( ir1.EQ.0 ) THEN
! 488: NOTROT = 0
! 489: PSKIPPED = 0
! 490: ISWROT = ISWROT + 1
! 491: END IF
! 492: *
! 493: IF( ROTOK ) THEN
! 494: *
! 495: AQOAP = AAQQ / AAPP
! 496: APOAQ = AAPP / AAQQ
! 497: THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1
! 498: *
! 499: IF( ABS( THETA ).GT.BIGTHETA ) THEN
! 500: *
! 501: T = HALF / THETA
! 502: CS = ONE
! 503:
! 504: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
! 505: $ CS, DCONJG(OMPQ)*T )
! 506: IF ( RSVEC ) THEN
! 507: CALL ZROT( MVL, V(1,p), 1,
! 508: $ V(1,q), 1, CS, DCONJG(OMPQ)*T )
! 509: END IF
! 510:
! 511: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 512: $ ONE+T*APOAQ*AAPQ1 ) )
! 513: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 514: $ ONE-T*AQOAP*AAPQ1 ) )
! 515: MXSINJ = DMAX1( MXSINJ, ABS( T ) )
! 516: *
! 517: ELSE
! 518: *
! 519: * .. choose correct signum for THETA and rotate
! 520: *
! 521: THSIGN = -DSIGN( ONE, AAPQ1 )
! 522: T = ONE / ( THETA+THSIGN*
! 523: $ DSQRT( ONE+THETA*THETA ) )
! 524: CS = DSQRT( ONE / ( ONE+T*T ) )
! 525: SN = T*CS
! 526: *
! 527: MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
! 528: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 529: $ ONE+T*APOAQ*AAPQ1 ) )
! 530: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 531: $ ONE-T*AQOAP*AAPQ1 ) )
! 532: *
! 533: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
! 534: $ CS, DCONJG(OMPQ)*SN )
! 535: IF ( RSVEC ) THEN
! 536: CALL ZROT( MVL, V(1,p), 1,
! 537: $ V(1,q), 1, CS, DCONJG(OMPQ)*SN )
! 538: END IF
! 539: END IF
! 540: D(p) = -D(q) * OMPQ
! 541: *
! 542: ELSE
! 543: * .. have to use modified Gram-Schmidt like transformation
! 544: CALL ZCOPY( M, A( 1, p ), 1,
! 545: $ WORK, 1 )
! 546: CALL ZLASCL( 'G', 0, 0, AAPP, ONE, M,
! 547: $ 1, WORK, LDA,
! 548: $ IERR )
! 549: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE, M,
! 550: $ 1, A( 1, q ), LDA, IERR )
! 551: CALL ZAXPY( M, -AAPQ, WORK, 1,
! 552: $ A( 1, q ), 1 )
! 553: CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M,
! 554: $ 1, A( 1, q ), LDA, IERR )
! 555: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 556: $ ONE-AAPQ1*AAPQ1 ) )
! 557: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 558: END IF
! 559: * END IF ROTOK THEN ... ELSE
! 560: *
! 561: * In the case of cancellation in updating SVA(q), SVA(p)
! 562: * recompute SVA(q), SVA(p).
! 563: *
! 564: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 565: $ THEN
! 566: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 567: $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 568: SVA( q ) = DZNRM2( M, A( 1, q ), 1 )
! 569: ELSE
! 570: T = ZERO
! 571: AAQQ = ONE
! 572: CALL ZLASSQ( M, A( 1, q ), 1, T,
! 573: $ AAQQ )
! 574: SVA( q ) = T*DSQRT( AAQQ )
! 575: END IF
! 576: END IF
! 577: IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
! 578: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 579: $ ( AAPP.GT.ROOTSFMIN ) ) THEN
! 580: AAPP = DZNRM2( M, A( 1, p ), 1 )
! 581: ELSE
! 582: T = ZERO
! 583: AAPP = ONE
! 584: CALL ZLASSQ( M, A( 1, p ), 1, T,
! 585: $ AAPP )
! 586: AAPP = T*DSQRT( AAPP )
! 587: END IF
! 588: SVA( p ) = AAPP
! 589: END IF
! 590: *
! 591: ELSE
! 592: * A(:,p) and A(:,q) already numerically orthogonal
! 593: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
! 594: *[RTD] SKIPPED = SKIPPED + 1
! 595: PSKIPPED = PSKIPPED + 1
! 596: END IF
! 597: ELSE
! 598: * A(:,q) is zero column
! 599: IF( ir1.EQ.0 )NOTROT = NOTROT + 1
! 600: PSKIPPED = PSKIPPED + 1
! 601: END IF
! 602: *
! 603: IF( ( i.LE.SWBAND ) .AND.
! 604: $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 605: IF( ir1.EQ.0 )AAPP = -AAPP
! 606: NOTROT = 0
! 607: GO TO 2103
! 608: END IF
! 609: *
! 610: 2002 CONTINUE
! 611: * END q-LOOP
! 612: *
! 613: 2103 CONTINUE
! 614: * bailed out of q-loop
! 615: *
! 616: SVA( p ) = AAPP
! 617: *
! 618: ELSE
! 619: SVA( p ) = AAPP
! 620: IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
! 621: $ NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
! 622: END IF
! 623: *
! 624: 2001 CONTINUE
! 625: * end of the p-loop
! 626: * end of doing the block ( ibr, ibr )
! 627: 1002 CONTINUE
! 628: * end of ir1-loop
! 629: *
! 630: * ... go to the off diagonal blocks
! 631: *
! 632: igl = ( ibr-1 )*KBL + 1
! 633: *
! 634: DO 2010 jbc = ibr + 1, NBL
! 635: *
! 636: jgl = ( jbc-1 )*KBL + 1
! 637: *
! 638: * doing the block at ( ibr, jbc )
! 639: *
! 640: IJBLSK = 0
! 641: DO 2100 p = igl, MIN0( igl+KBL-1, N )
! 642: *
! 643: AAPP = SVA( p )
! 644: IF( AAPP.GT.ZERO ) THEN
! 645: *
! 646: PSKIPPED = 0
! 647: *
! 648: DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
! 649: *
! 650: AAQQ = SVA( q )
! 651: IF( AAQQ.GT.ZERO ) THEN
! 652: AAPP0 = AAPP
! 653: *
! 654: * .. M x 2 Jacobi SVD ..
! 655: *
! 656: * Safe Gram matrix computation
! 657: *
! 658: IF( AAQQ.GE.ONE ) THEN
! 659: IF( AAPP.GE.AAQQ ) THEN
! 660: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 661: ELSE
! 662: ROTOK = ( SMALL*AAQQ ).LE.AAPP
! 663: END IF
! 664: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
! 665: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
! 666: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
! 667: ELSE
! 668: CALL ZCOPY( M, A( 1, p ), 1,
! 669: $ WORK, 1 )
! 670: CALL ZLASCL( 'G', 0, 0, AAPP,
! 671: $ ONE, M, 1,
! 672: $ WORK, LDA, IERR )
! 673: AAPQ = ZDOTC( M, WORK, 1,
! 674: $ A( 1, q ), 1 ) / AAQQ
! 675: END IF
! 676: ELSE
! 677: IF( AAPP.GE.AAQQ ) THEN
! 678: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 679: ELSE
! 680: ROTOK = AAQQ.LE.( AAPP / SMALL )
! 681: END IF
! 682: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 683: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
! 684: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
! 685: ELSE
! 686: CALL ZCOPY( M, A( 1, q ), 1,
! 687: $ WORK, 1 )
! 688: CALL ZLASCL( 'G', 0, 0, AAQQ,
! 689: $ ONE, M, 1,
! 690: $ WORK, LDA, IERR )
! 691: AAPQ = ZDOTC( M, A( 1, p ), 1,
! 692: $ WORK, 1 ) / AAPP
! 693: END IF
! 694: END IF
! 695: *
! 696: OMPQ = AAPQ / ABS(AAPQ)
! 697: * AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q)
! 698: AAPQ1 = -ABS(AAPQ)
! 699: MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
! 700: *
! 701: * TO rotate or NOT to rotate, THAT is the question ...
! 702: *
! 703: IF( ABS( AAPQ1 ).GT.TOL ) THEN
! 704: NOTROT = 0
! 705: *[RTD] ROTATED = ROTATED + 1
! 706: PSKIPPED = 0
! 707: ISWROT = ISWROT + 1
! 708: *
! 709: IF( ROTOK ) THEN
! 710: *
! 711: AQOAP = AAQQ / AAPP
! 712: APOAQ = AAPP / AAQQ
! 713: THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1
! 714: IF( AAQQ.GT.AAPP0 )THETA = -THETA
! 715: *
! 716: IF( ABS( THETA ).GT.BIGTHETA ) THEN
! 717: T = HALF / THETA
! 718: CS = ONE
! 719: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
! 720: $ CS, DCONJG(OMPQ)*T )
! 721: IF( RSVEC ) THEN
! 722: CALL ZROT( MVL, V(1,p), 1,
! 723: $ V(1,q), 1, CS, DCONJG(OMPQ)*T )
! 724: END IF
! 725: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 726: $ ONE+T*APOAQ*AAPQ1 ) )
! 727: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 728: $ ONE-T*AQOAP*AAPQ1 ) )
! 729: MXSINJ = DMAX1( MXSINJ, ABS( T ) )
! 730: ELSE
! 731: *
! 732: * .. choose correct signum for THETA and rotate
! 733: *
! 734: THSIGN = -DSIGN( ONE, AAPQ1 )
! 735: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
! 736: T = ONE / ( THETA+THSIGN*
! 737: $ DSQRT( ONE+THETA*THETA ) )
! 738: CS = DSQRT( ONE / ( ONE+T*T ) )
! 739: SN = T*CS
! 740: MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
! 741: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 742: $ ONE+T*APOAQ*AAPQ1 ) )
! 743: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 744: $ ONE-T*AQOAP*AAPQ1 ) )
! 745: *
! 746: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
! 747: $ CS, DCONJG(OMPQ)*SN )
! 748: IF( RSVEC ) THEN
! 749: CALL ZROT( MVL, V(1,p), 1,
! 750: $ V(1,q), 1, CS, DCONJG(OMPQ)*SN )
! 751: END IF
! 752: END IF
! 753: D(p) = -D(q) * OMPQ
! 754: *
! 755: ELSE
! 756: * .. have to use modified Gram-Schmidt like transformation
! 757: IF( AAPP.GT.AAQQ ) THEN
! 758: CALL ZCOPY( M, A( 1, p ), 1,
! 759: $ WORK, 1 )
! 760: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
! 761: $ M, 1, WORK,LDA,
! 762: $ IERR )
! 763: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
! 764: $ M, 1, A( 1, q ), LDA,
! 765: $ IERR )
! 766: CALL ZAXPY( M, -AAPQ, WORK,
! 767: $ 1, A( 1, q ), 1 )
! 768: CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
! 769: $ M, 1, A( 1, q ), LDA,
! 770: $ IERR )
! 771: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 772: $ ONE-AAPQ1*AAPQ1 ) )
! 773: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 774: ELSE
! 775: CALL ZCOPY( M, A( 1, q ), 1,
! 776: $ WORK, 1 )
! 777: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
! 778: $ M, 1, WORK,LDA,
! 779: $ IERR )
! 780: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
! 781: $ M, 1, A( 1, p ), LDA,
! 782: $ IERR )
! 783: CALL ZAXPY( M, -DCONJG(AAPQ),
! 784: $ WORK, 1, A( 1, p ), 1 )
! 785: CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
! 786: $ M, 1, A( 1, p ), LDA,
! 787: $ IERR )
! 788: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
! 789: $ ONE-AAPQ1*AAPQ1 ) )
! 790: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 791: END IF
! 792: END IF
! 793: * END IF ROTOK THEN ... ELSE
! 794: *
! 795: * In the case of cancellation in updating SVA(q), SVA(p)
! 796: * .. recompute SVA(q), SVA(p)
! 797: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 798: $ THEN
! 799: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 800: $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 801: SVA( q ) = DZNRM2( M, A( 1, q ), 1)
! 802: ELSE
! 803: T = ZERO
! 804: AAQQ = ONE
! 805: CALL ZLASSQ( M, A( 1, q ), 1, T,
! 806: $ AAQQ )
! 807: SVA( q ) = T*DSQRT( AAQQ )
! 808: END IF
! 809: END IF
! 810: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
! 811: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 812: $ ( AAPP.GT.ROOTSFMIN ) ) THEN
! 813: AAPP = DZNRM2( M, A( 1, p ), 1 )
! 814: ELSE
! 815: T = ZERO
! 816: AAPP = ONE
! 817: CALL ZLASSQ( M, A( 1, p ), 1, T,
! 818: $ AAPP )
! 819: AAPP = T*DSQRT( AAPP )
! 820: END IF
! 821: SVA( p ) = AAPP
! 822: END IF
! 823: * end of OK rotation
! 824: ELSE
! 825: NOTROT = NOTROT + 1
! 826: *[RTD] SKIPPED = SKIPPED + 1
! 827: PSKIPPED = PSKIPPED + 1
! 828: IJBLSK = IJBLSK + 1
! 829: END IF
! 830: ELSE
! 831: NOTROT = NOTROT + 1
! 832: PSKIPPED = PSKIPPED + 1
! 833: IJBLSK = IJBLSK + 1
! 834: END IF
! 835: *
! 836: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
! 837: $ THEN
! 838: SVA( p ) = AAPP
! 839: NOTROT = 0
! 840: GO TO 2011
! 841: END IF
! 842: IF( ( i.LE.SWBAND ) .AND.
! 843: $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 844: AAPP = -AAPP
! 845: NOTROT = 0
! 846: GO TO 2203
! 847: END IF
! 848: *
! 849: 2200 CONTINUE
! 850: * end of the q-loop
! 851: 2203 CONTINUE
! 852: *
! 853: SVA( p ) = AAPP
! 854: *
! 855: ELSE
! 856: *
! 857: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
! 858: $ MIN0( jgl+KBL-1, N ) - jgl + 1
! 859: IF( AAPP.LT.ZERO )NOTROT = 0
! 860: *
! 861: END IF
! 862: *
! 863: 2100 CONTINUE
! 864: * end of the p-loop
! 865: 2010 CONTINUE
! 866: * end of the jbc-loop
! 867: 2011 CONTINUE
! 868: *2011 bailed out of the jbc-loop
! 869: DO 2012 p = igl, MIN0( igl+KBL-1, N )
! 870: SVA( p ) = ABS( SVA( p ) )
! 871: 2012 CONTINUE
! 872: ***
! 873: 2000 CONTINUE
! 874: *2000 :: end of the ibr-loop
! 875: *
! 876: * .. update SVA(N)
! 877: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
! 878: $ THEN
! 879: SVA( N ) = DZNRM2( M, A( 1, N ), 1 )
! 880: ELSE
! 881: T = ZERO
! 882: AAPP = ONE
! 883: CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
! 884: SVA( N ) = T*DSQRT( AAPP )
! 885: END IF
! 886: *
! 887: * Additional steering devices
! 888: *
! 889: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
! 890: $ ( ISWROT.LE.N ) ) )SWBAND = i
! 891: *
! 892: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DFLOAT( N ) )*
! 893: $ TOL ) .AND. ( DFLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
! 894: GO TO 1994
! 895: END IF
! 896: *
! 897: IF( NOTROT.GE.EMPTSW )GO TO 1994
! 898: *
! 899: 1993 CONTINUE
! 900: * end i=1:NSWEEP loop
! 901: *
! 902: * #:( Reaching this point means that the procedure has not converged.
! 903: INFO = NSWEEP - 1
! 904: GO TO 1995
! 905: *
! 906: 1994 CONTINUE
! 907: * #:) Reaching this point means numerical convergence after the i-th
! 908: * sweep.
! 909: *
! 910: INFO = 0
! 911: * #:) INFO = 0 confirms successful iterations.
! 912: 1995 CONTINUE
! 913: *
! 914: * Sort the vector SVA() of column norms.
! 915: DO 5991 p = 1, N - 1
! 916: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 917: IF( p.NE.q ) THEN
! 918: TEMP1 = SVA( p )
! 919: SVA( p ) = SVA( q )
! 920: SVA( q ) = TEMP1
! 921: AAPQ = D( p )
! 922: D( p ) = D( q )
! 923: D( q ) = AAPQ
! 924: CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 925: IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
! 926: END IF
! 927: 5991 CONTINUE
! 928: *
! 929: RETURN
! 930: * ..
! 931: * .. END OF ZGSVJ0
! 932: * ..
! 933: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgsvj0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack