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