Annotation of rpl/lapack/lapack/ztrevc3.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZTREVC3
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZTREVC3 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc3.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc3.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc3.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
! 22: * VR, LDVR, MM, M, WORK, LWORK, RWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER HOWMNY, SIDE
! 26: * INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * LOGICAL SELECT( * )
! 30: * DOUBLE PRECISION RWORK( * )
! 31: * COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
! 32: * $ WORK( * )
! 33: * ..
! 34: *
! 35: *
! 36: *> \par Purpose:
! 37: * =============
! 38: *>
! 39: *> \verbatim
! 40: *>
! 41: *> ZTREVC3 computes some or all of the right and/or left eigenvectors of
! 42: *> a complex upper triangular matrix T.
! 43: *> Matrices of this type are produced by the Schur factorization of
! 44: *> a complex general matrix: A = Q*T*Q**H, as computed by ZHSEQR.
! 45: *>
! 46: *> The right eigenvector x and the left eigenvector y of T corresponding
! 47: *> to an eigenvalue w are defined by:
! 48: *>
! 49: *> T*x = w*x, (y**H)*T = w*(y**H)
! 50: *>
! 51: *> where y**H denotes the conjugate transpose of the vector y.
! 52: *> The eigenvalues are not input to this routine, but are read directly
! 53: *> from the diagonal of T.
! 54: *>
! 55: *> This routine returns the matrices X and/or Y of right and left
! 56: *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
! 57: *> input matrix. If Q is the unitary factor that reduces a matrix A to
! 58: *> Schur form T, then Q*X and Q*Y are the matrices of right and left
! 59: *> eigenvectors of A.
! 60: *>
! 61: *> This uses a Level 3 BLAS version of the back transformation.
! 62: *> \endverbatim
! 63: *
! 64: * Arguments:
! 65: * ==========
! 66: *
! 67: *> \param[in] SIDE
! 68: *> \verbatim
! 69: *> SIDE is CHARACTER*1
! 70: *> = 'R': compute right eigenvectors only;
! 71: *> = 'L': compute left eigenvectors only;
! 72: *> = 'B': compute both right and left eigenvectors.
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[in] HOWMNY
! 76: *> \verbatim
! 77: *> HOWMNY is CHARACTER*1
! 78: *> = 'A': compute all right and/or left eigenvectors;
! 79: *> = 'B': compute all right and/or left eigenvectors,
! 80: *> backtransformed using the matrices supplied in
! 81: *> VR and/or VL;
! 82: *> = 'S': compute selected right and/or left eigenvectors,
! 83: *> as indicated by the logical array SELECT.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[in] SELECT
! 87: *> \verbatim
! 88: *> SELECT is LOGICAL array, dimension (N)
! 89: *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
! 90: *> computed.
! 91: *> The eigenvector corresponding to the j-th eigenvalue is
! 92: *> computed if SELECT(j) = .TRUE..
! 93: *> Not referenced if HOWMNY = 'A' or 'B'.
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[in] N
! 97: *> \verbatim
! 98: *> N is INTEGER
! 99: *> The order of the matrix T. N >= 0.
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[in,out] T
! 103: *> \verbatim
! 104: *> T is COMPLEX*16 array, dimension (LDT,N)
! 105: *> The upper triangular matrix T. T is modified, but restored
! 106: *> on exit.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[in] LDT
! 110: *> \verbatim
! 111: *> LDT is INTEGER
! 112: *> The leading dimension of the array T. LDT >= max(1,N).
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in,out] VL
! 116: *> \verbatim
! 117: *> VL is COMPLEX*16 array, dimension (LDVL,MM)
! 118: *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
! 119: *> contain an N-by-N matrix Q (usually the unitary matrix Q of
! 120: *> Schur vectors returned by ZHSEQR).
! 121: *> On exit, if SIDE = 'L' or 'B', VL contains:
! 122: *> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
! 123: *> if HOWMNY = 'B', the matrix Q*Y;
! 124: *> if HOWMNY = 'S', the left eigenvectors of T specified by
! 125: *> SELECT, stored consecutively in the columns
! 126: *> of VL, in the same order as their
! 127: *> eigenvalues.
! 128: *> Not referenced if SIDE = 'R'.
! 129: *> \endverbatim
! 130: *>
! 131: *> \param[in] LDVL
! 132: *> \verbatim
! 133: *> LDVL is INTEGER
! 134: *> The leading dimension of the array VL.
! 135: *> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[in,out] VR
! 139: *> \verbatim
! 140: *> VR is COMPLEX*16 array, dimension (LDVR,MM)
! 141: *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
! 142: *> contain an N-by-N matrix Q (usually the unitary matrix Q of
! 143: *> Schur vectors returned by ZHSEQR).
! 144: *> On exit, if SIDE = 'R' or 'B', VR contains:
! 145: *> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
! 146: *> if HOWMNY = 'B', the matrix Q*X;
! 147: *> if HOWMNY = 'S', the right eigenvectors of T specified by
! 148: *> SELECT, stored consecutively in the columns
! 149: *> of VR, in the same order as their
! 150: *> eigenvalues.
! 151: *> Not referenced if SIDE = 'L'.
! 152: *> \endverbatim
! 153: *>
! 154: *> \param[in] LDVR
! 155: *> \verbatim
! 156: *> LDVR is INTEGER
! 157: *> The leading dimension of the array VR.
! 158: *> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
! 159: *> \endverbatim
! 160: *>
! 161: *> \param[in] MM
! 162: *> \verbatim
! 163: *> MM is INTEGER
! 164: *> The number of columns in the arrays VL and/or VR. MM >= M.
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[out] M
! 168: *> \verbatim
! 169: *> M is INTEGER
! 170: *> The number of columns in the arrays VL and/or VR actually
! 171: *> used to store the eigenvectors.
! 172: *> If HOWMNY = 'A' or 'B', M is set to N.
! 173: *> Each selected eigenvector occupies one column.
! 174: *> \endverbatim
! 175: *>
! 176: *> \param[out] WORK
! 177: *> \verbatim
! 178: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
! 179: *> \endverbatim
! 180: *>
! 181: *> \param[in] LWORK
! 182: *> \verbatim
! 183: *> LWORK is INTEGER
! 184: *> The dimension of array WORK. LWORK >= max(1,2*N).
! 185: *> For optimum performance, LWORK >= N + 2*N*NB, where NB is
! 186: *> the optimal blocksize.
! 187: *>
! 188: *> If LWORK = -1, then a workspace query is assumed; the routine
! 189: *> only calculates the optimal size of the WORK array, returns
! 190: *> this value as the first entry of the WORK array, and no error
! 191: *> message related to LWORK is issued by XERBLA.
! 192: *> \endverbatim
! 193: *>
! 194: *> \param[out] RWORK
! 195: *> \verbatim
! 196: *> RWORK is DOUBLE PRECISION array, dimension (LRWORK)
! 197: *> \endverbatim
! 198: *>
! 199: *> \param[in] LRWORK
! 200: *> \verbatim
! 201: *> LRWORK is INTEGER
! 202: *> The dimension of array RWORK. LRWORK >= max(1,N).
! 203: *>
! 204: *> If LRWORK = -1, then a workspace query is assumed; the routine
! 205: *> only calculates the optimal size of the RWORK array, returns
! 206: *> this value as the first entry of the RWORK array, and no error
! 207: *> message related to LRWORK is issued by XERBLA.
! 208: *> \endverbatim
! 209: *>
! 210: *> \param[out] INFO
! 211: *> \verbatim
! 212: *> INFO is INTEGER
! 213: *> = 0: successful exit
! 214: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 215: *> \endverbatim
! 216: *
! 217: * Authors:
! 218: * ========
! 219: *
! 220: *> \author Univ. of Tennessee
! 221: *> \author Univ. of California Berkeley
! 222: *> \author Univ. of Colorado Denver
! 223: *> \author NAG Ltd.
! 224: *
! 225: *> \date November 2011
! 226: *
! 227: * @precisions fortran z -> c
! 228: *
! 229: *> \ingroup complex16OTHERcomputational
! 230: *
! 231: *> \par Further Details:
! 232: * =====================
! 233: *>
! 234: *> \verbatim
! 235: *>
! 236: *> The algorithm used in this program is basically backward (forward)
! 237: *> substitution, with scaling to make the the code robust against
! 238: *> possible overflow.
! 239: *>
! 240: *> Each eigenvector is normalized so that the element of largest
! 241: *> magnitude has magnitude 1; here the magnitude of a complex number
! 242: *> (x,y) is taken to be |x| + |y|.
! 243: *> \endverbatim
! 244: *>
! 245: * =====================================================================
! 246: SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
! 247: $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
! 248: IMPLICIT NONE
! 249: *
! 250: * -- LAPACK computational routine (version 3.4.0) --
! 251: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 252: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 253: * November 2011
! 254: *
! 255: * .. Scalar Arguments ..
! 256: CHARACTER HOWMNY, SIDE
! 257: INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
! 258: * ..
! 259: * .. Array Arguments ..
! 260: LOGICAL SELECT( * )
! 261: DOUBLE PRECISION RWORK( * )
! 262: COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
! 263: $ WORK( * )
! 264: * ..
! 265: *
! 266: * =====================================================================
! 267: *
! 268: * .. Parameters ..
! 269: DOUBLE PRECISION ZERO, ONE
! 270: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 271: COMPLEX*16 CZERO, CONE
! 272: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 273: $ CONE = ( 1.0D+0, 0.0D+0 ) )
! 274: INTEGER NBMIN, NBMAX
! 275: PARAMETER ( NBMIN = 8, NBMAX = 128 )
! 276: * ..
! 277: * .. Local Scalars ..
! 278: LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
! 279: INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB
! 280: DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
! 281: COMPLEX*16 CDUM
! 282: * ..
! 283: * .. External Functions ..
! 284: LOGICAL LSAME
! 285: INTEGER ILAENV, IZAMAX
! 286: DOUBLE PRECISION DLAMCH, DZASUM
! 287: EXTERNAL LSAME, ILAENV, IZAMAX, DLAMCH, DZASUM
! 288: * ..
! 289: * .. External Subroutines ..
! 290: EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS
! 291: * ..
! 292: * .. Intrinsic Functions ..
! 293: INTRINSIC ABS, DBLE, DCMPLX, CONJG, AIMAG, MAX
! 294: * ..
! 295: * .. Statement Functions ..
! 296: DOUBLE PRECISION CABS1
! 297: * ..
! 298: * .. Statement Function definitions ..
! 299: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( AIMAG( CDUM ) )
! 300: * ..
! 301: * .. Executable Statements ..
! 302: *
! 303: * Decode and test the input parameters
! 304: *
! 305: BOTHV = LSAME( SIDE, 'B' )
! 306: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
! 307: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
! 308: *
! 309: ALLV = LSAME( HOWMNY, 'A' )
! 310: OVER = LSAME( HOWMNY, 'B' )
! 311: SOMEV = LSAME( HOWMNY, 'S' )
! 312: *
! 313: * Set M to the number of columns required to store the selected
! 314: * eigenvectors.
! 315: *
! 316: IF( SOMEV ) THEN
! 317: M = 0
! 318: DO 10 J = 1, N
! 319: IF( SELECT( J ) )
! 320: $ M = M + 1
! 321: 10 CONTINUE
! 322: ELSE
! 323: M = N
! 324: END IF
! 325: *
! 326: INFO = 0
! 327: NB = ILAENV( 1, 'ZTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
! 328: MAXWRK = N + 2*N*NB
! 329: WORK(1) = MAXWRK
! 330: RWORK(1) = N
! 331: LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
! 332: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
! 333: INFO = -1
! 334: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
! 335: INFO = -2
! 336: ELSE IF( N.LT.0 ) THEN
! 337: INFO = -4
! 338: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
! 339: INFO = -6
! 340: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
! 341: INFO = -8
! 342: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
! 343: INFO = -10
! 344: ELSE IF( MM.LT.M ) THEN
! 345: INFO = -11
! 346: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
! 347: INFO = -14
! 348: ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
! 349: INFO = -16
! 350: END IF
! 351: IF( INFO.NE.0 ) THEN
! 352: CALL XERBLA( 'ZTREVC3', -INFO )
! 353: RETURN
! 354: ELSE IF( LQUERY ) THEN
! 355: RETURN
! 356: END IF
! 357: *
! 358: * Quick return if possible.
! 359: *
! 360: IF( N.EQ.0 )
! 361: $ RETURN
! 362: *
! 363: * Use blocked version of back-transformation if sufficient workspace.
! 364: * Zero-out the workspace to avoid potential NaN propagation.
! 365: *
! 366: IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
! 367: NB = (LWORK - N) / (2*N)
! 368: NB = MIN( NB, NBMAX )
! 369: CALL ZLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N )
! 370: ELSE
! 371: NB = 1
! 372: END IF
! 373: *
! 374: * Set the constants to control overflow.
! 375: *
! 376: UNFL = DLAMCH( 'Safe minimum' )
! 377: OVFL = ONE / UNFL
! 378: CALL DLABAD( UNFL, OVFL )
! 379: ULP = DLAMCH( 'Precision' )
! 380: SMLNUM = UNFL*( N / ULP )
! 381: *
! 382: * Store the diagonal elements of T in working array WORK.
! 383: *
! 384: DO 20 I = 1, N
! 385: WORK( I ) = T( I, I )
! 386: 20 CONTINUE
! 387: *
! 388: * Compute 1-norm of each column of strictly upper triangular
! 389: * part of T to control overflow in triangular solver.
! 390: *
! 391: RWORK( 1 ) = ZERO
! 392: DO 30 J = 2, N
! 393: RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
! 394: 30 CONTINUE
! 395: *
! 396: IF( RIGHTV ) THEN
! 397: *
! 398: * ============================================================
! 399: * Compute right eigenvectors.
! 400: *
! 401: * IV is index of column in current block.
! 402: * Non-blocked version always uses IV=NB=1;
! 403: * blocked version starts with IV=NB, goes down to 1.
! 404: * (Note the "0-th" column is used to store the original diagonal.)
! 405: IV = NB
! 406: IS = M
! 407: DO 80 KI = N, 1, -1
! 408: IF( SOMEV ) THEN
! 409: IF( .NOT.SELECT( KI ) )
! 410: $ GO TO 80
! 411: END IF
! 412: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
! 413: *
! 414: * --------------------------------------------------------
! 415: * Complex right eigenvector
! 416: *
! 417: WORK( KI + IV*N ) = CONE
! 418: *
! 419: * Form right-hand side.
! 420: *
! 421: DO 40 K = 1, KI - 1
! 422: WORK( K + IV*N ) = -T( K, KI )
! 423: 40 CONTINUE
! 424: *
! 425: * Solve upper triangular system:
! 426: * [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
! 427: *
! 428: DO 50 K = 1, KI - 1
! 429: T( K, K ) = T( K, K ) - T( KI, KI )
! 430: IF( CABS1( T( K, K ) ).LT.SMIN )
! 431: $ T( K, K ) = SMIN
! 432: 50 CONTINUE
! 433: *
! 434: IF( KI.GT.1 ) THEN
! 435: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
! 436: $ KI-1, T, LDT, WORK( 1 + IV*N ), SCALE,
! 437: $ RWORK, INFO )
! 438: WORK( KI + IV*N ) = SCALE
! 439: END IF
! 440: *
! 441: * Copy the vector x or Q*x to VR and normalize.
! 442: *
! 443: IF( .NOT.OVER ) THEN
! 444: * ------------------------------
! 445: * no back-transform: copy x to VR and normalize.
! 446: CALL ZCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
! 447: *
! 448: II = IZAMAX( KI, VR( 1, IS ), 1 )
! 449: REMAX = ONE / CABS1( VR( II, IS ) )
! 450: CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
! 451: *
! 452: DO 60 K = KI + 1, N
! 453: VR( K, IS ) = CZERO
! 454: 60 CONTINUE
! 455: *
! 456: ELSE IF( NB.EQ.1 ) THEN
! 457: * ------------------------------
! 458: * version 1: back-transform each vector with GEMV, Q*x.
! 459: IF( KI.GT.1 )
! 460: $ CALL ZGEMV( 'N', N, KI-1, CONE, VR, LDVR,
! 461: $ WORK( 1 + IV*N ), 1, DCMPLX( SCALE ),
! 462: $ VR( 1, KI ), 1 )
! 463: *
! 464: II = IZAMAX( N, VR( 1, KI ), 1 )
! 465: REMAX = ONE / CABS1( VR( II, KI ) )
! 466: CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
! 467: *
! 468: ELSE
! 469: * ------------------------------
! 470: * version 2: back-transform block of vectors with GEMM
! 471: * zero out below vector
! 472: DO K = KI + 1, N
! 473: WORK( K + IV*N ) = CZERO
! 474: END DO
! 475: *
! 476: * Columns IV:NB of work are valid vectors.
! 477: * When the number of vectors stored reaches NB,
! 478: * or if this was last vector, do the GEMM
! 479: IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN
! 480: CALL ZGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE,
! 481: $ VR, LDVR,
! 482: $ WORK( 1 + (IV)*N ), N,
! 483: $ CZERO,
! 484: $ WORK( 1 + (NB+IV)*N ), N )
! 485: * normalize vectors
! 486: DO K = IV, NB
! 487: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
! 488: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
! 489: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
! 490: END DO
! 491: CALL ZLACPY( 'F', N, NB-IV+1,
! 492: $ WORK( 1 + (NB+IV)*N ), N,
! 493: $ VR( 1, KI ), LDVR )
! 494: IV = NB
! 495: ELSE
! 496: IV = IV - 1
! 497: END IF
! 498: END IF
! 499: *
! 500: * Restore the original diagonal elements of T.
! 501: *
! 502: DO 70 K = 1, KI - 1
! 503: T( K, K ) = WORK( K )
! 504: 70 CONTINUE
! 505: *
! 506: IS = IS - 1
! 507: 80 CONTINUE
! 508: END IF
! 509: *
! 510: IF( LEFTV ) THEN
! 511: *
! 512: * ============================================================
! 513: * Compute left eigenvectors.
! 514: *
! 515: * IV is index of column in current block.
! 516: * Non-blocked version always uses IV=1;
! 517: * blocked version starts with IV=1, goes up to NB.
! 518: * (Note the "0-th" column is used to store the original diagonal.)
! 519: IV = 1
! 520: IS = 1
! 521: DO 130 KI = 1, N
! 522: *
! 523: IF( SOMEV ) THEN
! 524: IF( .NOT.SELECT( KI ) )
! 525: $ GO TO 130
! 526: END IF
! 527: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
! 528: *
! 529: * --------------------------------------------------------
! 530: * Complex left eigenvector
! 531: *
! 532: WORK( KI + IV*N ) = CONE
! 533: *
! 534: * Form right-hand side.
! 535: *
! 536: DO 90 K = KI + 1, N
! 537: WORK( K + IV*N ) = -CONJG( T( KI, K ) )
! 538: 90 CONTINUE
! 539: *
! 540: * Solve conjugate-transposed triangular system:
! 541: * [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
! 542: *
! 543: DO 100 K = KI + 1, N
! 544: T( K, K ) = T( K, K ) - T( KI, KI )
! 545: IF( CABS1( T( K, K ) ).LT.SMIN )
! 546: $ T( K, K ) = SMIN
! 547: 100 CONTINUE
! 548: *
! 549: IF( KI.LT.N ) THEN
! 550: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
! 551: $ 'Y', N-KI, T( KI+1, KI+1 ), LDT,
! 552: $ WORK( KI+1 + IV*N ), SCALE, RWORK, INFO )
! 553: WORK( KI + IV*N ) = SCALE
! 554: END IF
! 555: *
! 556: * Copy the vector x or Q*x to VL and normalize.
! 557: *
! 558: IF( .NOT.OVER ) THEN
! 559: * ------------------------------
! 560: * no back-transform: copy x to VL and normalize.
! 561: CALL ZCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 )
! 562: *
! 563: II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
! 564: REMAX = ONE / CABS1( VL( II, IS ) )
! 565: CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
! 566: *
! 567: DO 110 K = 1, KI - 1
! 568: VL( K, IS ) = CZERO
! 569: 110 CONTINUE
! 570: *
! 571: ELSE IF( NB.EQ.1 ) THEN
! 572: * ------------------------------
! 573: * version 1: back-transform each vector with GEMV, Q*x.
! 574: IF( KI.LT.N )
! 575: $ CALL ZGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL,
! 576: $ WORK( KI+1 + IV*N ), 1, DCMPLX( SCALE ),
! 577: $ VL( 1, KI ), 1 )
! 578: *
! 579: II = IZAMAX( N, VL( 1, KI ), 1 )
! 580: REMAX = ONE / CABS1( VL( II, KI ) )
! 581: CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
! 582: *
! 583: ELSE
! 584: * ------------------------------
! 585: * version 2: back-transform block of vectors with GEMM
! 586: * zero out above vector
! 587: * could go from KI-NV+1 to KI-1
! 588: DO K = 1, KI - 1
! 589: WORK( K + IV*N ) = CZERO
! 590: END DO
! 591: *
! 592: * Columns 1:IV of work are valid vectors.
! 593: * When the number of vectors stored reaches NB,
! 594: * or if this was last vector, do the GEMM
! 595: IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN
! 596: CALL ZGEMM( 'N', 'N', N, IV, N-KI+IV, ONE,
! 597: $ VL( 1, KI-IV+1 ), LDVL,
! 598: $ WORK( KI-IV+1 + (1)*N ), N,
! 599: $ CZERO,
! 600: $ WORK( 1 + (NB+1)*N ), N )
! 601: * normalize vectors
! 602: DO K = 1, IV
! 603: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
! 604: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
! 605: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
! 606: END DO
! 607: CALL ZLACPY( 'F', N, IV,
! 608: $ WORK( 1 + (NB+1)*N ), N,
! 609: $ VL( 1, KI-IV+1 ), LDVL )
! 610: IV = 1
! 611: ELSE
! 612: IV = IV + 1
! 613: END IF
! 614: END IF
! 615: *
! 616: * Restore the original diagonal elements of T.
! 617: *
! 618: DO 120 K = KI + 1, N
! 619: T( K, K ) = WORK( K )
! 620: 120 CONTINUE
! 621: *
! 622: IS = IS + 1
! 623: 130 CONTINUE
! 624: END IF
! 625: *
! 626: RETURN
! 627: *
! 628: * End of ZTREVC3
! 629: *
! 630: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztrevc3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack