Annotation of rpl/lapack/lapack/zgeev.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
! 2: $ WORK, LWORK, RWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER JOBVL, JOBVR
! 11: INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION RWORK( * )
! 15: COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
! 16: $ W( * ), WORK( * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
! 23: * eigenvalues and, optionally, the left and/or right eigenvectors.
! 24: *
! 25: * The right eigenvector v(j) of A satisfies
! 26: * A * v(j) = lambda(j) * v(j)
! 27: * where lambda(j) is its eigenvalue.
! 28: * The left eigenvector u(j) of A satisfies
! 29: * u(j)**H * A = lambda(j) * u(j)**H
! 30: * where u(j)**H denotes the conjugate transpose of u(j).
! 31: *
! 32: * The computed eigenvectors are normalized to have Euclidean norm
! 33: * equal to 1 and largest component real.
! 34: *
! 35: * Arguments
! 36: * =========
! 37: *
! 38: * JOBVL (input) CHARACTER*1
! 39: * = 'N': left eigenvectors of A are not computed;
! 40: * = 'V': left eigenvectors of are computed.
! 41: *
! 42: * JOBVR (input) CHARACTER*1
! 43: * = 'N': right eigenvectors of A are not computed;
! 44: * = 'V': right eigenvectors of A are computed.
! 45: *
! 46: * N (input) INTEGER
! 47: * The order of the matrix A. N >= 0.
! 48: *
! 49: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 50: * On entry, the N-by-N matrix A.
! 51: * On exit, A has been overwritten.
! 52: *
! 53: * LDA (input) INTEGER
! 54: * The leading dimension of the array A. LDA >= max(1,N).
! 55: *
! 56: * W (output) COMPLEX*16 array, dimension (N)
! 57: * W contains the computed eigenvalues.
! 58: *
! 59: * VL (output) COMPLEX*16 array, dimension (LDVL,N)
! 60: * If JOBVL = 'V', the left eigenvectors u(j) are stored one
! 61: * after another in the columns of VL, in the same order
! 62: * as their eigenvalues.
! 63: * If JOBVL = 'N', VL is not referenced.
! 64: * u(j) = VL(:,j), the j-th column of VL.
! 65: *
! 66: * LDVL (input) INTEGER
! 67: * The leading dimension of the array VL. LDVL >= 1; if
! 68: * JOBVL = 'V', LDVL >= N.
! 69: *
! 70: * VR (output) COMPLEX*16 array, dimension (LDVR,N)
! 71: * If JOBVR = 'V', the right eigenvectors v(j) are stored one
! 72: * after another in the columns of VR, in the same order
! 73: * as their eigenvalues.
! 74: * If JOBVR = 'N', VR is not referenced.
! 75: * v(j) = VR(:,j), the j-th column of VR.
! 76: *
! 77: * LDVR (input) INTEGER
! 78: * The leading dimension of the array VR. LDVR >= 1; if
! 79: * JOBVR = 'V', LDVR >= N.
! 80: *
! 81: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 82: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 83: *
! 84: * LWORK (input) INTEGER
! 85: * The dimension of the array WORK. LWORK >= max(1,2*N).
! 86: * For good performance, LWORK must generally be larger.
! 87: *
! 88: * If LWORK = -1, then a workspace query is assumed; the routine
! 89: * only calculates the optimal size of the WORK array, returns
! 90: * this value as the first entry of the WORK array, and no error
! 91: * message related to LWORK is issued by XERBLA.
! 92: *
! 93: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 94: *
! 95: * INFO (output) INTEGER
! 96: * = 0: successful exit
! 97: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 98: * > 0: if INFO = i, the QR algorithm failed to compute all the
! 99: * eigenvalues, and no eigenvectors have been computed;
! 100: * elements and i+1:N of W contain eigenvalues which have
! 101: * converged.
! 102: *
! 103: * =====================================================================
! 104: *
! 105: * .. Parameters ..
! 106: DOUBLE PRECISION ZERO, ONE
! 107: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 108: * ..
! 109: * .. Local Scalars ..
! 110: LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
! 111: CHARACTER SIDE
! 112: INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
! 113: $ IWRK, K, MAXWRK, MINWRK, NOUT
! 114: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
! 115: COMPLEX*16 TMP
! 116: * ..
! 117: * .. Local Arrays ..
! 118: LOGICAL SELECT( 1 )
! 119: DOUBLE PRECISION DUM( 1 )
! 120: * ..
! 121: * .. External Subroutines ..
! 122: EXTERNAL DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,
! 123: $ ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR
! 124: * ..
! 125: * .. External Functions ..
! 126: LOGICAL LSAME
! 127: INTEGER IDAMAX, ILAENV
! 128: DOUBLE PRECISION DLAMCH, DZNRM2, ZLANGE
! 129: EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
! 130: * ..
! 131: * .. Intrinsic Functions ..
! 132: INTRINSIC DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
! 133: * ..
! 134: * .. Executable Statements ..
! 135: *
! 136: * Test the input arguments
! 137: *
! 138: INFO = 0
! 139: LQUERY = ( LWORK.EQ.-1 )
! 140: WANTVL = LSAME( JOBVL, 'V' )
! 141: WANTVR = LSAME( JOBVR, 'V' )
! 142: IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
! 143: INFO = -1
! 144: ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
! 145: INFO = -2
! 146: ELSE IF( N.LT.0 ) THEN
! 147: INFO = -3
! 148: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 149: INFO = -5
! 150: ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
! 151: INFO = -8
! 152: ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
! 153: INFO = -10
! 154: END IF
! 155: *
! 156: * Compute workspace
! 157: * (Note: Comments in the code beginning "Workspace:" describe the
! 158: * minimal amount of workspace needed at that point in the code,
! 159: * as well as the preferred amount for good performance.
! 160: * CWorkspace refers to complex workspace, and RWorkspace to real
! 161: * workspace. NB refers to the optimal block size for the
! 162: * immediately following subroutine, as returned by ILAENV.
! 163: * HSWORK refers to the workspace preferred by ZHSEQR, as
! 164: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
! 165: * the worst case.)
! 166: *
! 167: IF( INFO.EQ.0 ) THEN
! 168: IF( N.EQ.0 ) THEN
! 169: MINWRK = 1
! 170: MAXWRK = 1
! 171: ELSE
! 172: MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
! 173: MINWRK = 2*N
! 174: IF( WANTVL ) THEN
! 175: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
! 176: $ ' ', N, 1, N, -1 ) )
! 177: CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
! 178: $ WORK, -1, INFO )
! 179: ELSE IF( WANTVR ) THEN
! 180: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
! 181: $ ' ', N, 1, N, -1 ) )
! 182: CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
! 183: $ WORK, -1, INFO )
! 184: ELSE
! 185: CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
! 186: $ WORK, -1, INFO )
! 187: END IF
! 188: HSWORK = WORK( 1 )
! 189: MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
! 190: END IF
! 191: WORK( 1 ) = MAXWRK
! 192: *
! 193: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
! 194: INFO = -12
! 195: END IF
! 196: END IF
! 197: *
! 198: IF( INFO.NE.0 ) THEN
! 199: CALL XERBLA( 'ZGEEV ', -INFO )
! 200: RETURN
! 201: ELSE IF( LQUERY ) THEN
! 202: RETURN
! 203: END IF
! 204: *
! 205: * Quick return if possible
! 206: *
! 207: IF( N.EQ.0 )
! 208: $ RETURN
! 209: *
! 210: * Get machine constants
! 211: *
! 212: EPS = DLAMCH( 'P' )
! 213: SMLNUM = DLAMCH( 'S' )
! 214: BIGNUM = ONE / SMLNUM
! 215: CALL DLABAD( SMLNUM, BIGNUM )
! 216: SMLNUM = SQRT( SMLNUM ) / EPS
! 217: BIGNUM = ONE / SMLNUM
! 218: *
! 219: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 220: *
! 221: ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
! 222: SCALEA = .FALSE.
! 223: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 224: SCALEA = .TRUE.
! 225: CSCALE = SMLNUM
! 226: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 227: SCALEA = .TRUE.
! 228: CSCALE = BIGNUM
! 229: END IF
! 230: IF( SCALEA )
! 231: $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
! 232: *
! 233: * Balance the matrix
! 234: * (CWorkspace: none)
! 235: * (RWorkspace: need N)
! 236: *
! 237: IBAL = 1
! 238: CALL ZGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
! 239: *
! 240: * Reduce to upper Hessenberg form
! 241: * (CWorkspace: need 2*N, prefer N+N*NB)
! 242: * (RWorkspace: none)
! 243: *
! 244: ITAU = 1
! 245: IWRK = ITAU + N
! 246: CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
! 247: $ LWORK-IWRK+1, IERR )
! 248: *
! 249: IF( WANTVL ) THEN
! 250: *
! 251: * Want left eigenvectors
! 252: * Copy Householder vectors to VL
! 253: *
! 254: SIDE = 'L'
! 255: CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
! 256: *
! 257: * Generate unitary matrix in VL
! 258: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
! 259: * (RWorkspace: none)
! 260: *
! 261: CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
! 262: $ LWORK-IWRK+1, IERR )
! 263: *
! 264: * Perform QR iteration, accumulating Schur vectors in VL
! 265: * (CWorkspace: need 1, prefer HSWORK (see comments) )
! 266: * (RWorkspace: none)
! 267: *
! 268: IWRK = ITAU
! 269: CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
! 270: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
! 271: *
! 272: IF( WANTVR ) THEN
! 273: *
! 274: * Want left and right eigenvectors
! 275: * Copy Schur vectors to VR
! 276: *
! 277: SIDE = 'B'
! 278: CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
! 279: END IF
! 280: *
! 281: ELSE IF( WANTVR ) THEN
! 282: *
! 283: * Want right eigenvectors
! 284: * Copy Householder vectors to VR
! 285: *
! 286: SIDE = 'R'
! 287: CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
! 288: *
! 289: * Generate unitary matrix in VR
! 290: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
! 291: * (RWorkspace: none)
! 292: *
! 293: CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
! 294: $ LWORK-IWRK+1, IERR )
! 295: *
! 296: * Perform QR iteration, accumulating Schur vectors in VR
! 297: * (CWorkspace: need 1, prefer HSWORK (see comments) )
! 298: * (RWorkspace: none)
! 299: *
! 300: IWRK = ITAU
! 301: CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
! 302: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
! 303: *
! 304: ELSE
! 305: *
! 306: * Compute eigenvalues only
! 307: * (CWorkspace: need 1, prefer HSWORK (see comments) )
! 308: * (RWorkspace: none)
! 309: *
! 310: IWRK = ITAU
! 311: CALL ZHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
! 312: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
! 313: END IF
! 314: *
! 315: * If INFO > 0 from ZHSEQR, then quit
! 316: *
! 317: IF( INFO.GT.0 )
! 318: $ GO TO 50
! 319: *
! 320: IF( WANTVL .OR. WANTVR ) THEN
! 321: *
! 322: * Compute left and/or right eigenvectors
! 323: * (CWorkspace: need 2*N)
! 324: * (RWorkspace: need 2*N)
! 325: *
! 326: IRWORK = IBAL + N
! 327: CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
! 328: $ N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
! 329: END IF
! 330: *
! 331: IF( WANTVL ) THEN
! 332: *
! 333: * Undo balancing of left eigenvectors
! 334: * (CWorkspace: none)
! 335: * (RWorkspace: need N)
! 336: *
! 337: CALL ZGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
! 338: $ IERR )
! 339: *
! 340: * Normalize left eigenvectors and make largest component real
! 341: *
! 342: DO 20 I = 1, N
! 343: SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
! 344: CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
! 345: DO 10 K = 1, N
! 346: RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 +
! 347: $ DIMAG( VL( K, I ) )**2
! 348: 10 CONTINUE
! 349: K = IDAMAX( N, RWORK( IRWORK ), 1 )
! 350: TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
! 351: CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
! 352: VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
! 353: 20 CONTINUE
! 354: END IF
! 355: *
! 356: IF( WANTVR ) THEN
! 357: *
! 358: * Undo balancing of right eigenvectors
! 359: * (CWorkspace: none)
! 360: * (RWorkspace: need N)
! 361: *
! 362: CALL ZGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
! 363: $ IERR )
! 364: *
! 365: * Normalize right eigenvectors and make largest component real
! 366: *
! 367: DO 40 I = 1, N
! 368: SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
! 369: CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
! 370: DO 30 K = 1, N
! 371: RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 +
! 372: $ DIMAG( VR( K, I ) )**2
! 373: 30 CONTINUE
! 374: K = IDAMAX( N, RWORK( IRWORK ), 1 )
! 375: TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
! 376: CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
! 377: VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
! 378: 40 CONTINUE
! 379: END IF
! 380: *
! 381: * Undo scaling if necessary
! 382: *
! 383: 50 CONTINUE
! 384: IF( SCALEA ) THEN
! 385: CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
! 386: $ MAX( N-INFO, 1 ), IERR )
! 387: IF( INFO.GT.0 ) THEN
! 388: CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
! 389: END IF
! 390: END IF
! 391: *
! 392: WORK( 1 ) = MAXWRK
! 393: RETURN
! 394: *
! 395: * End of ZGEEV
! 396: *
! 397: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgeev.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack