Annotation of rpl/lapack/lapack/zgsvj1.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivots.
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZGSVJ1 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgsvj1.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgsvj1.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgsvj1.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
! 22: * EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * DOUBLE PRECISION EPS, SFMIN, TOL
! 26: * INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
! 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: *> ZGSVJ1 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 targets only particular pivots and it does not check convergence
! 43: *> (stopping criterion). Few tunning parameters (marked by [TP]) are
! 44: *> available for the implementer.
! 45: *>
! 46: *> Further Details
! 47: *> ~~~~~~~~~~~~~~~
! 48: *> ZGSVJ1 applies few sweeps of Jacobi rotations in the column space of
! 49: *> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
! 50: *> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
! 51: *> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
! 52: *> [x]'s in the following scheme:
! 53: *>
! 54: *> | * * * [x] [x] [x]|
! 55: *> | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
! 56: *> | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
! 57: *> |[x] [x] [x] * * * |
! 58: *> |[x] [x] [x] * * * |
! 59: *> |[x] [x] [x] * * * |
! 60: *>
! 61: *> In terms of the columns of A, the first N1 columns are rotated 'against'
! 62: *> the remaining N-N1 columns, trying to increase the angle between the
! 63: *> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
! 64: *> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
! 65: *> The number of sweeps is given in NSWEEP and the orthogonality threshold
! 66: *> is given in TOL.
! 67: *> \endverbatim
! 68: *
! 69: * Arguments:
! 70: * ==========
! 71: *
! 72: *> \param[in] JOBV
! 73: *> \verbatim
! 74: *> JOBV is CHARACTER*1
! 75: *> Specifies whether the output from this procedure is used
! 76: *> to compute the matrix V:
! 77: *> = 'V': the product of the Jacobi rotations is accumulated
! 78: *> by postmulyiplying the N-by-N array V.
! 79: *> (See the description of V.)
! 80: *> = 'A': the product of the Jacobi rotations is accumulated
! 81: *> by postmulyiplying the MV-by-N array V.
! 82: *> (See the descriptions of MV and V.)
! 83: *> = 'N': the Jacobi rotations are not accumulated.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[in] M
! 87: *> \verbatim
! 88: *> M is INTEGER
! 89: *> The number of rows of the input matrix A. M >= 0.
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] N
! 93: *> \verbatim
! 94: *> N is INTEGER
! 95: *> The number of columns of the input matrix A.
! 96: *> M >= N >= 0.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] N1
! 100: *> \verbatim
! 101: *> N1 is INTEGER
! 102: *> N1 specifies the 2 x 2 block partition, the first N1 columns are
! 103: *> rotated 'against' the remaining N-N1 columns of A.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[in,out] A
! 107: *> \verbatim
! 108: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 109: *> On entry, M-by-N matrix A, such that A*diag(D) represents
! 110: *> the input matrix.
! 111: *> On exit,
! 112: *> A_onexit * D_onexit represents the input matrix A*diag(D)
! 113: *> post-multiplied by a sequence of Jacobi rotations, where the
! 114: *> rotation threshold and the total number of sweeps are given in
! 115: *> TOL and NSWEEP, respectively.
! 116: *> (See the descriptions of N1, D, TOL and NSWEEP.)
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in] LDA
! 120: *> \verbatim
! 121: *> LDA is INTEGER
! 122: *> The leading dimension of the array A. LDA >= max(1,M).
! 123: *> \endverbatim
! 124: *>
! 125: *> \param[in,out] D
! 126: *> \verbatim
! 127: *> D is DOUBLE PRECISION array, dimension (N)
! 128: *> The array D accumulates the scaling factors from the fast scaled
! 129: *> Jacobi rotations.
! 130: *> On entry, A*diag(D) represents the input matrix.
! 131: *> On exit, A_onexit*diag(D_onexit) represents the input matrix
! 132: *> post-multiplied by a sequence of Jacobi rotations, where the
! 133: *> rotation threshold and the total number of sweeps are given in
! 134: *> TOL and NSWEEP, respectively.
! 135: *> (See the descriptions of N1, A, TOL and NSWEEP.)
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[in,out] SVA
! 139: *> \verbatim
! 140: *> SVA is DOUBLE PRECISION array, dimension (N)
! 141: *> On entry, SVA contains the Euclidean norms of the columns of
! 142: *> the matrix A*diag(D).
! 143: *> On exit, SVA contains the Euclidean norms of the columns of
! 144: *> the matrix onexit*diag(D_onexit).
! 145: *> \endverbatim
! 146: *>
! 147: *> \param[in] MV
! 148: *> \verbatim
! 149: *> MV is INTEGER
! 150: *> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
! 151: *> sequence of Jacobi rotations.
! 152: *> If JOBV = 'N', then MV is not referenced.
! 153: *> \endverbatim
! 154: *>
! 155: *> \param[in,out] V
! 156: *> \verbatim
! 157: *> V is DOUBLE PRECISION array, dimension (LDV,N)
! 158: *> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
! 159: *> sequence of Jacobi rotations.
! 160: *> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
! 161: *> sequence of Jacobi rotations.
! 162: *> If JOBV = 'N', then V is not referenced.
! 163: *> \endverbatim
! 164: *>
! 165: *> \param[in] LDV
! 166: *> \verbatim
! 167: *> LDV is INTEGER
! 168: *> The leading dimension of the array V, LDV >= 1.
! 169: *> If JOBV = 'V', LDV .GE. N.
! 170: *> If JOBV = 'A', LDV .GE. MV.
! 171: *> \endverbatim
! 172: *>
! 173: *> \param[in] EPS
! 174: *> \verbatim
! 175: *> EPS is DOUBLE PRECISION
! 176: *> EPS = DLAMCH('Epsilon')
! 177: *> \endverbatim
! 178: *>
! 179: *> \param[in] SFMIN
! 180: *> \verbatim
! 181: *> SFMIN is DOUBLE PRECISION
! 182: *> SFMIN = DLAMCH('Safe Minimum')
! 183: *> \endverbatim
! 184: *>
! 185: *> \param[in] TOL
! 186: *> \verbatim
! 187: *> TOL is DOUBLE PRECISION
! 188: *> TOL is the threshold for Jacobi rotations. For a pair
! 189: *> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
! 190: *> applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
! 191: *> \endverbatim
! 192: *>
! 193: *> \param[in] NSWEEP
! 194: *> \verbatim
! 195: *> NSWEEP is INTEGER
! 196: *> NSWEEP is the number of sweeps of Jacobi rotations to be
! 197: *> performed.
! 198: *> \endverbatim
! 199: *>
! 200: *> \param[out] WORK
! 201: *> \verbatim
! 202: *> WORK is DOUBLE PRECISION array, dimension (LWORK)
! 203: *> \endverbatim
! 204: *>
! 205: *> \param[in] LWORK
! 206: *> \verbatim
! 207: *> LWORK is INTEGER
! 208: *> LWORK is the dimension of WORK. LWORK .GE. M.
! 209: *> \endverbatim
! 210: *>
! 211: *> \param[out] INFO
! 212: *> \verbatim
! 213: *> INFO is INTEGER
! 214: *> = 0 : successful exit.
! 215: *> < 0 : if INFO = -i, then the i-th argument had an illegal value
! 216: *> \endverbatim
! 217: *
! 218: * Authors:
! 219: * ========
! 220: *
! 221: *> \author Univ. of Tennessee
! 222: *> \author Univ. of California Berkeley
! 223: *> \author Univ. of Colorado Denver
! 224: *> \author NAG Ltd.
! 225: *
! 226: *> \date November 2015
! 227: *
! 228: *> \ingroup complex16OTHERcomputational
! 229: *
! 230: *> \par Contributors:
! 231: * ==================
! 232: *>
! 233: *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
! 234: *
! 235: * =====================================================================
! 236: SUBROUTINE ZGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
! 237: $ EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
! 238: *
! 239: * -- LAPACK computational routine (version 3.6.0) --
! 240: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 241: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 242: * November 2015
! 243: *
! 244: IMPLICIT NONE
! 245: * .. Scalar Arguments ..
! 246: DOUBLE PRECISION EPS, SFMIN, TOL
! 247: INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
! 248: CHARACTER*1 JOBV
! 249: * ..
! 250: * .. Array Arguments ..
! 251: COMPLEX*16 A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
! 252: DOUBLE PRECISION SVA( N )
! 253: * ..
! 254: *
! 255: * =====================================================================
! 256: *
! 257: * .. Local Parameters ..
! 258: DOUBLE PRECISION ZERO, HALF, ONE
! 259: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0)
! 260: * ..
! 261: * .. Local Scalars ..
! 262: COMPLEX*16 AAPQ, OMPQ
! 263: DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
! 264: $ BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
! 265: $ ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
! 266: $ TEMP1, THETA, THSIGN
! 267: INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
! 268: $ ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
! 269: $ p, PSKIPPED, q, ROWSKIP, SWBAND
! 270: LOGICAL APPLV, ROTOK, RSVEC
! 271: * ..
! 272: * ..
! 273: * .. Intrinsic Functions ..
! 274: INTRINSIC ABS, DCONJG, DMAX1, DFLOAT, MIN0, DSIGN, DSQRT
! 275: * ..
! 276: * .. External Functions ..
! 277: DOUBLE PRECISION DZNRM2
! 278: COMPLEX*16 ZDOTC
! 279: INTEGER IDAMAX
! 280: LOGICAL LSAME
! 281: EXTERNAL IDAMAX, LSAME, ZDOTC, DZNRM2
! 282: * ..
! 283: * .. External Subroutines ..
! 284: * .. from BLAS
! 285: EXTERNAL ZCOPY, ZROT, ZSWAP
! 286: * .. from LAPACK
! 287: EXTERNAL ZLASCL, ZLASSQ, XERBLA
! 288: * ..
! 289: * .. Executable Statements ..
! 290: *
! 291: * Test the input parameters.
! 292: *
! 293: APPLV = LSAME( JOBV, 'A' )
! 294: RSVEC = LSAME( JOBV, 'V' )
! 295: IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
! 296: INFO = -1
! 297: ELSE IF( M.LT.0 ) THEN
! 298: INFO = -2
! 299: ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
! 300: INFO = -3
! 301: ELSE IF( N1.LT.0 ) THEN
! 302: INFO = -4
! 303: ELSE IF( LDA.LT.M ) THEN
! 304: INFO = -6
! 305: ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
! 306: INFO = -9
! 307: ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
! 308: $ ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
! 309: INFO = -11
! 310: ELSE IF( TOL.LE.EPS ) THEN
! 311: INFO = -14
! 312: ELSE IF( NSWEEP.LT.0 ) THEN
! 313: INFO = -15
! 314: ELSE IF( LWORK.LT.M ) THEN
! 315: INFO = -17
! 316: ELSE
! 317: INFO = 0
! 318: END IF
! 319: *
! 320: * #:(
! 321: IF( INFO.NE.0 ) THEN
! 322: CALL XERBLA( 'ZGSVJ1', -INFO )
! 323: RETURN
! 324: END IF
! 325: *
! 326: IF( RSVEC ) THEN
! 327: MVL = N
! 328: ELSE IF( APPLV ) THEN
! 329: MVL = MV
! 330: END IF
! 331: RSVEC = RSVEC .OR. APPLV
! 332:
! 333: ROOTEPS = DSQRT( EPS )
! 334: ROOTSFMIN = DSQRT( SFMIN )
! 335: SMALL = SFMIN / EPS
! 336: BIG = ONE / SFMIN
! 337: ROOTBIG = ONE / ROOTSFMIN
! 338: LARGE = BIG / DSQRT( DFLOAT( M*N ) )
! 339: BIGTHETA = ONE / ROOTEPS
! 340: ROOTTOL = DSQRT( TOL )
! 341: *
! 342: * .. Initialize the right singular vector matrix ..
! 343: *
! 344: * RSVEC = LSAME( JOBV, 'Y' )
! 345: *
! 346: EMPTSW = N1*( N-N1 )
! 347: NOTROT = 0
! 348: *
! 349: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
! 350: *
! 351: KBL = MIN0( 8, N )
! 352: NBLR = N1 / KBL
! 353: IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
! 354:
! 355: * .. the tiling is nblr-by-nblc [tiles]
! 356:
! 357: NBLC = ( N-N1 ) / KBL
! 358: IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
! 359: BLSKIP = ( KBL**2 ) + 1
! 360: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
! 361:
! 362: ROWSKIP = MIN0( 5, KBL )
! 363: *[TP] ROWSKIP is a tuning parameter.
! 364: SWBAND = 0
! 365: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
! 366: * if ZGESVJ is used as a computational routine in the preconditioned
! 367: * Jacobi SVD algorithm ZGEJSV.
! 368: *
! 369: *
! 370: * | * * * [x] [x] [x]|
! 371: * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
! 372: * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
! 373: * |[x] [x] [x] * * * |
! 374: * |[x] [x] [x] * * * |
! 375: * |[x] [x] [x] * * * |
! 376: *
! 377: *
! 378: DO 1993 i = 1, NSWEEP
! 379: *
! 380: * .. go go go ...
! 381: *
! 382: MXAAPQ = ZERO
! 383: MXSINJ = ZERO
! 384: ISWROT = 0
! 385: *
! 386: NOTROT = 0
! 387: PSKIPPED = 0
! 388: *
! 389: * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
! 390: * 1 <= p < q <= N. This is the first step toward a blocked implementation
! 391: * of the rotations. New implementation, based on block transformations,
! 392: * is under development.
! 393: *
! 394: DO 2000 ibr = 1, NBLR
! 395: *
! 396: igl = ( ibr-1 )*KBL + 1
! 397: *
! 398:
! 399: *
! 400: * ... go to the off diagonal blocks
! 401: *
! 402: igl = ( ibr-1 )*KBL + 1
! 403: *
! 404: * DO 2010 jbc = ibr + 1, NBL
! 405: DO 2010 jbc = 1, NBLC
! 406: *
! 407: jgl = ( jbc-1 )*KBL + N1 + 1
! 408: *
! 409: * doing the block at ( ibr, jbc )
! 410: *
! 411: IJBLSK = 0
! 412: DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
! 413: *
! 414: AAPP = SVA( p )
! 415: IF( AAPP.GT.ZERO ) THEN
! 416: *
! 417: PSKIPPED = 0
! 418: *
! 419: DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
! 420: *
! 421: AAQQ = SVA( q )
! 422: IF( AAQQ.GT.ZERO ) THEN
! 423: AAPP0 = AAPP
! 424: *
! 425: * .. M x 2 Jacobi SVD ..
! 426: *
! 427: * Safe Gram matrix computation
! 428: *
! 429: IF( AAQQ.GE.ONE ) THEN
! 430: IF( AAPP.GE.AAQQ ) THEN
! 431: ROTOK = ( SMALL*AAPP ).LE.AAQQ
! 432: ELSE
! 433: ROTOK = ( SMALL*AAQQ ).LE.AAPP
! 434: END IF
! 435: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
! 436: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
! 437: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
! 438: ELSE
! 439: CALL ZCOPY( M, A( 1, p ), 1,
! 440: $ WORK, 1 )
! 441: CALL ZLASCL( 'G', 0, 0, AAPP,
! 442: $ ONE, M, 1,
! 443: $ WORK, LDA, IERR )
! 444: AAPQ = ZDOTC( M, WORK, 1,
! 445: $ A( 1, q ), 1 ) / AAQQ
! 446: END IF
! 447: ELSE
! 448: IF( AAPP.GE.AAQQ ) THEN
! 449: ROTOK = AAPP.LE.( AAQQ / SMALL )
! 450: ELSE
! 451: ROTOK = AAQQ.LE.( AAPP / SMALL )
! 452: END IF
! 453: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
! 454: AAPQ = ( ZDOTC( M, A( 1, p ), 1,
! 455: $ A( 1, q ), 1 ) / AAQQ ) / AAPP
! 456: ELSE
! 457: CALL ZCOPY( M, A( 1, q ), 1,
! 458: $ WORK, 1 )
! 459: CALL ZLASCL( 'G', 0, 0, AAQQ,
! 460: $ ONE, M, 1,
! 461: $ WORK, LDA, IERR )
! 462: AAPQ = ZDOTC( M, A( 1, p ), 1,
! 463: $ WORK, 1 ) / AAPP
! 464: END IF
! 465: END IF
! 466: *
! 467: OMPQ = AAPQ / ABS(AAPQ)
! 468: * AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q)
! 469: AAPQ1 = -ABS(AAPQ)
! 470: MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
! 471: *
! 472: * TO rotate or NOT to rotate, THAT is the question ...
! 473: *
! 474: IF( ABS( AAPQ1 ).GT.TOL ) THEN
! 475: NOTROT = 0
! 476: *[RTD] ROTATED = ROTATED + 1
! 477: PSKIPPED = 0
! 478: ISWROT = ISWROT + 1
! 479: *
! 480: IF( ROTOK ) THEN
! 481: *
! 482: AQOAP = AAQQ / AAPP
! 483: APOAQ = AAPP / AAQQ
! 484: THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1
! 485: IF( AAQQ.GT.AAPP0 )THETA = -THETA
! 486: *
! 487: IF( ABS( THETA ).GT.BIGTHETA ) THEN
! 488: T = HALF / THETA
! 489: CS = ONE
! 490: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
! 491: $ CS, DCONJG(OMPQ)*T )
! 492: IF( RSVEC ) THEN
! 493: CALL ZROT( MVL, V(1,p), 1,
! 494: $ V(1,q), 1, CS, DCONJG(OMPQ)*T )
! 495: END IF
! 496: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 497: $ ONE+T*APOAQ*AAPQ1 ) )
! 498: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 499: $ ONE-T*AQOAP*AAPQ1 ) )
! 500: MXSINJ = DMAX1( MXSINJ, ABS( T ) )
! 501: ELSE
! 502: *
! 503: * .. choose correct signum for THETA and rotate
! 504: *
! 505: THSIGN = -DSIGN( ONE, AAPQ1 )
! 506: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
! 507: T = ONE / ( THETA+THSIGN*
! 508: $ DSQRT( ONE+THETA*THETA ) )
! 509: CS = DSQRT( ONE / ( ONE+T*T ) )
! 510: SN = T*CS
! 511: MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
! 512: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 513: $ ONE+T*APOAQ*AAPQ1 ) )
! 514: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 515: $ ONE-T*AQOAP*AAPQ1 ) )
! 516: *
! 517: CALL ZROT( M, A(1,p), 1, A(1,q), 1,
! 518: $ CS, DCONJG(OMPQ)*SN )
! 519: IF( RSVEC ) THEN
! 520: CALL ZROT( MVL, V(1,p), 1,
! 521: $ V(1,q), 1, CS, DCONJG(OMPQ)*SN )
! 522: END IF
! 523: END IF
! 524: D(p) = -D(q) * OMPQ
! 525: *
! 526: ELSE
! 527: * .. have to use modified Gram-Schmidt like transformation
! 528: IF( AAPP.GT.AAQQ ) THEN
! 529: CALL ZCOPY( M, A( 1, p ), 1,
! 530: $ WORK, 1 )
! 531: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
! 532: $ M, 1, WORK,LDA,
! 533: $ IERR )
! 534: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
! 535: $ M, 1, A( 1, q ), LDA,
! 536: $ IERR )
! 537: CALL ZAXPY( M, -AAPQ, WORK,
! 538: $ 1, A( 1, q ), 1 )
! 539: CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
! 540: $ M, 1, A( 1, q ), LDA,
! 541: $ IERR )
! 542: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
! 543: $ ONE-AAPQ1*AAPQ1 ) )
! 544: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 545: ELSE
! 546: CALL ZCOPY( M, A( 1, q ), 1,
! 547: $ WORK, 1 )
! 548: CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
! 549: $ M, 1, WORK,LDA,
! 550: $ IERR )
! 551: CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
! 552: $ M, 1, A( 1, p ), LDA,
! 553: $ IERR )
! 554: CALL ZAXPY( M, -DCONJG(AAPQ),
! 555: $ WORK, 1, A( 1, p ), 1 )
! 556: CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
! 557: $ M, 1, A( 1, p ), LDA,
! 558: $ IERR )
! 559: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
! 560: $ ONE-AAPQ1*AAPQ1 ) )
! 561: MXSINJ = DMAX1( MXSINJ, SFMIN )
! 562: END IF
! 563: END IF
! 564: * END IF ROTOK THEN ... ELSE
! 565: *
! 566: * In the case of cancellation in updating SVA(q), SVA(p)
! 567: * .. recompute SVA(q), SVA(p)
! 568: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
! 569: $ THEN
! 570: IF( ( AAQQ.LT.ROOTBIG ) .AND.
! 571: $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
! 572: SVA( q ) = DZNRM2( M, A( 1, q ), 1)
! 573: ELSE
! 574: T = ZERO
! 575: AAQQ = ONE
! 576: CALL ZLASSQ( M, A( 1, q ), 1, T,
! 577: $ AAQQ )
! 578: SVA( q ) = T*DSQRT( AAQQ )
! 579: END IF
! 580: END IF
! 581: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
! 582: IF( ( AAPP.LT.ROOTBIG ) .AND.
! 583: $ ( AAPP.GT.ROOTSFMIN ) ) THEN
! 584: AAPP = DZNRM2( M, A( 1, p ), 1 )
! 585: ELSE
! 586: T = ZERO
! 587: AAPP = ONE
! 588: CALL ZLASSQ( M, A( 1, p ), 1, T,
! 589: $ AAPP )
! 590: AAPP = T*DSQRT( AAPP )
! 591: END IF
! 592: SVA( p ) = AAPP
! 593: END IF
! 594: * end of OK rotation
! 595: ELSE
! 596: NOTROT = NOTROT + 1
! 597: *[RTD] SKIPPED = SKIPPED + 1
! 598: PSKIPPED = PSKIPPED + 1
! 599: IJBLSK = IJBLSK + 1
! 600: END IF
! 601: ELSE
! 602: NOTROT = NOTROT + 1
! 603: PSKIPPED = PSKIPPED + 1
! 604: IJBLSK = IJBLSK + 1
! 605: END IF
! 606: *
! 607: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
! 608: $ THEN
! 609: SVA( p ) = AAPP
! 610: NOTROT = 0
! 611: GO TO 2011
! 612: END IF
! 613: IF( ( i.LE.SWBAND ) .AND.
! 614: $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
! 615: AAPP = -AAPP
! 616: NOTROT = 0
! 617: GO TO 2203
! 618: END IF
! 619: *
! 620: 2200 CONTINUE
! 621: * end of the q-loop
! 622: 2203 CONTINUE
! 623: *
! 624: SVA( p ) = AAPP
! 625: *
! 626: ELSE
! 627: *
! 628: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
! 629: $ MIN0( jgl+KBL-1, N ) - jgl + 1
! 630: IF( AAPP.LT.ZERO )NOTROT = 0
! 631: *
! 632: END IF
! 633: *
! 634: 2100 CONTINUE
! 635: * end of the p-loop
! 636: 2010 CONTINUE
! 637: * end of the jbc-loop
! 638: 2011 CONTINUE
! 639: *2011 bailed out of the jbc-loop
! 640: DO 2012 p = igl, MIN0( igl+KBL-1, N )
! 641: SVA( p ) = ABS( SVA( p ) )
! 642: 2012 CONTINUE
! 643: ***
! 644: 2000 CONTINUE
! 645: *2000 :: end of the ibr-loop
! 646: *
! 647: * .. update SVA(N)
! 648: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
! 649: $ THEN
! 650: SVA( N ) = DZNRM2( M, A( 1, N ), 1 )
! 651: ELSE
! 652: T = ZERO
! 653: AAPP = ONE
! 654: CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
! 655: SVA( N ) = T*DSQRT( AAPP )
! 656: END IF
! 657: *
! 658: * Additional steering devices
! 659: *
! 660: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
! 661: $ ( ISWROT.LE.N ) ) )SWBAND = i
! 662: *
! 663: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DFLOAT( N ) )*
! 664: $ TOL ) .AND. ( DFLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
! 665: GO TO 1994
! 666: END IF
! 667: *
! 668: IF( NOTROT.GE.EMPTSW )GO TO 1994
! 669: *
! 670: 1993 CONTINUE
! 671: * end i=1:NSWEEP loop
! 672: *
! 673: * #:( Reaching this point means that the procedure has not converged.
! 674: INFO = NSWEEP - 1
! 675: GO TO 1995
! 676: *
! 677: 1994 CONTINUE
! 678: * #:) Reaching this point means numerical convergence after the i-th
! 679: * sweep.
! 680: *
! 681: INFO = 0
! 682: * #:) INFO = 0 confirms successful iterations.
! 683: 1995 CONTINUE
! 684: *
! 685: * Sort the vector SVA() of column norms.
! 686: DO 5991 p = 1, N - 1
! 687: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
! 688: IF( p.NE.q ) THEN
! 689: TEMP1 = SVA( p )
! 690: SVA( p ) = SVA( q )
! 691: SVA( q ) = TEMP1
! 692: AAPQ = D( p )
! 693: D( p ) = D( q )
! 694: D( q ) = AAPQ
! 695: CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
! 696: IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
! 697: END IF
! 698: 5991 CONTINUE
! 699: *
! 700: *
! 701: RETURN
! 702: * ..
! 703: * .. END OF ZGSVJ1
! 704: * ..
! 705: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgsvj1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack