Rev 22555 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
SUBROUTINE RGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,$ LDVR, WORK, LWORK, INFO )** -- LAPACK driver routine (version 3.0) --* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,* Courant Institute, Argonne National Lab, and Rice University* December 8, 1999** .. Scalar Arguments ..CHARACTER JOBVL, JOBVRINTEGER INFO, LDA, LDVL, LDVR, LWORK, N* ..* .. Array Arguments ..DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),$ WI( * ), WORK( * ), WR( * )* ..** Purpose* =======** DGEEV computes for an N-by-N real nonsymmetric matrix A, the* eigenvalues and, optionally, the left and/or right eigenvectors.** The right eigenvector v(j) of A satisfies* A * v(j) = lambda(j) * v(j)* where lambda(j) is its eigenvalue.* The left eigenvector u(j) of A satisfies* u(j)**H * A = lambda(j) * u(j)**H* where u(j)**H denotes the conjugate transpose of u(j).** The computed eigenvectors are normalized to have Euclidean norm* equal to 1 and largest component real.** Arguments* =========** JOBVL (input) CHARACTER*1* = 'N': left eigenvectors of A are not computed;* = 'V': left eigenvectors of A are computed.** JOBVR (input) CHARACTER*1* = 'N': right eigenvectors of A are not computed;* = 'V': right eigenvectors of A are computed.** N (input) INTEGER* The order of the matrix A. N >= 0.** A (input/output) DOUBLE PRECISION array, dimension (LDA,N)* On entry, the N-by-N matrix A.* On exit, A has been overwritten.** LDA (input) INTEGER* The leading dimension of the array A. LDA >= max(1,N).** WR (output) DOUBLE PRECISION array, dimension (N)* WI (output) DOUBLE PRECISION array, dimension (N)* WR and WI contain the real and imaginary parts,* respectively, of the computed eigenvalues. Complex* conjugate pairs of eigenvalues appear consecutively* with the eigenvalue having the positive imaginary part* first.** VL (output) DOUBLE PRECISION array, dimension (LDVL,N)* If JOBVL = 'V', the left eigenvectors u(j) are stored one* after another in the columns of VL, in the same order* as their eigenvalues.* If JOBVL = 'N', VL is not referenced.* If the j-th eigenvalue is real, then u(j) = VL(:,j),* the j-th column of VL.* If the j-th and (j+1)-st eigenvalues form a complex* conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and* u(j+1) = VL(:,j) - i*VL(:,j+1).** LDVL (input) INTEGER* The leading dimension of the array VL. LDVL >= 1; if* JOBVL = 'V', LDVL >= N.** VR (output) DOUBLE PRECISION array, dimension (LDVR,N)* If JOBVR = 'V', the right eigenvectors v(j) are stored one* after another in the columns of VR, in the same order* as their eigenvalues.* If JOBVR = 'N', VR is not referenced.* If the j-th eigenvalue is real, then v(j) = VR(:,j),* the j-th column of VR.* If the j-th and (j+1)-st eigenvalues form a complex* conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and* v(j+1) = VR(:,j) - i*VR(:,j+1).** LDVR (input) INTEGER* The leading dimension of the array VR. LDVR >= 1; if* JOBVR = 'V', LDVR >= N.** WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.** LWORK (input) INTEGER* The dimension of the array WORK. LWORK >= max(1,3*N), and* if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good* performance, LWORK must generally be larger.** If LWORK = -1, then a workspace query is assumed; the routine* only calculates the optimal size of the WORK array, returns* this value as the first entry of the WORK array, and no error* message related to LWORK is issued by XERBLA.** INFO (output) INTEGER* = 0: successful exit* < 0: if INFO = -i, the i-th argument had an illegal value.* > 0: if INFO = i, the QR algorithm failed to compute all the* eigenvalues, and no eigenvectors have been computed;* elements i+1:N of WR and WI contain eigenvalues which* have converged.** =====================================================================** .. Parameters ..DOUBLE PRECISION ZERO, ONEPARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )* ..* .. Local Scalars ..LOGICAL LQUERY, SCALEA, WANTVL, WANTVRCHARACTER SIDEINTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,$ MAXB, MAXWRK, MINWRK, NOUTDOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,$ SN* ..* .. Local Arrays ..LOGICAL SELECT( 1 )DOUBLE PRECISION DUM( 1 )* ..* .. External Subroutines ..EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY, DLARTG,$ DLASCL, DORGHR, DROT, DSCAL, DTREVC, XERBLA* ..* .. External Functions ..LOGICAL LSAMERINTEGER IDAMAX, ILAENVDOUBLE PRECISION DLAMCH, RLANGE, DLAPY2, DNRM2EXTERNAL LSAMER, IDAMAX, ILAENV, DLAMCH, RLANGE, DLAPY2,$ DNRM2* ..* .. Intrinsic Functions ..INTRINSIC MAX, MIN, SQRT* ..* .. Executable Statements ..** Test the input arguments*INFO = 0LQUERY = ( LWORK.EQ.-1 )WANTVL = LSAMER( JOBVL, 'V' )WANTVR = LSAMER( JOBVR, 'V' )IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAMER( JOBVL, 'N' ) ) ) THENINFO = -1ELSE IF( ( .NOT.WANTVR ) .AND. (.NOT.LSAMER( JOBVR, 'N' ) ) ) THENINFO = -2ELSE IF( N.LT.0 ) THENINFO = -3ELSE IF( LDA.LT.MAX( 1, N ) ) THENINFO = -5ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THENINFO = -9ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THENINFO = -11END IF** Compute workspace* (Note: Comments in the code beginning "Workspace:" describe the* minimal amount of workspace needed at that point in the code,* as well as the preferred amount for good performance.* NB refers to the optimal block size for the immediately* following subroutine, as returned by ILAENV.* HSWORK refers to the workspace preferred by DHSEQR, as* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,* the worst case.)*MINWRK = 1IF( INFO.EQ.0 .AND. ( LWORK.GE.1 .OR. LQUERY ) ) THENMAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THENMINWRK = MAX( 1, 3*N )MAXB = MAX( ILAENV( 8, 'DHSEQR', 'EN', N, 1, N, -1 ), 2 )K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'EN', N, 1,$ N, -1 ) ) )HSWORK = MAX( K*( K+2 ), 2*N )MAXWRK = MAX( MAXWRK, N+1, N+HSWORK )ELSEMINWRK = MAX( 1, 4*N )MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*$ ILAENV( 1, 'DORGHR', ' ', N, 1, N, -1 ) )MAXB = MAX( ILAENV( 8, 'DHSEQR', 'SV', N, 1, N, -1 ), 2 )K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'SV', N, 1,$ N, -1 ) ) )HSWORK = MAX( K*( K+2 ), 2*N )MAXWRK = MAX( MAXWRK, N+1, N+HSWORK )MAXWRK = MAX( MAXWRK, 4*N )END IFWORK( 1 ) = MAXWRKEND IFIF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THENINFO = -13END IFIF( INFO.NE.0 ) THENCALL XERBLA( 'DGEEV ', -INFO )RETURNELSE IF( LQUERY ) THENRETURNEND IF** Quick return if possible*IF( N.EQ.0 )$ RETURN** Get machine constants*EPS = DLAMCH( 'P' )SMLNUM = DLAMCH( 'S' )BIGNUM = ONE / SMLNUMCALL DLABAD( SMLNUM, BIGNUM )SMLNUM = SQRT( SMLNUM ) / EPSBIGNUM = ONE / SMLNUM** Scale A if max element outside range [SMLNUM,BIGNUM]*ANRM = RLANGE( 'M', N, N, A, LDA, DUM )SCALEA = .FALSE.IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THENSCALEA = .TRUE.CSCALE = SMLNUMELSE IF( ANRM.GT.BIGNUM ) THENSCALEA = .TRUE.CSCALE = BIGNUMEND IFIF( SCALEA )$ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )** Balance the matrix* (Workspace: need N)*IBAL = 1CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )** Reduce to upper Hessenberg form* (Workspace: need 3*N, prefer 2*N+N*NB)*ITAU = IBAL + NIWRK = ITAU + NCALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),$ LWORK-IWRK+1, IERR )*IF( WANTVL ) THEN** Want left eigenvectors* Copy Householder vectors to VL*SIDE = 'L'CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )** Generate orthogonal matrix in VL* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)*CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),$ LWORK-IWRK+1, IERR )** Perform QR iteration, accumulating Schur vectors in VL* (Workspace: need N+1, prefer N+HSWORK (see comments) )*IWRK = ITAUCALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,$ WORK( IWRK ), LWORK-IWRK+1, INFO )*IF( WANTVR ) THEN** Want left and right eigenvectors* Copy Schur vectors to VR*SIDE = 'B'CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )END IF*ELSE IF( WANTVR ) THEN** Want right eigenvectors* Copy Householder vectors to VR*SIDE = 'R'CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )** Generate orthogonal matrix in VR* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)*CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),$ LWORK-IWRK+1, IERR )** Perform QR iteration, accumulating Schur vectors in VR* (Workspace: need N+1, prefer N+HSWORK (see comments) )*IWRK = ITAUCALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,$ WORK( IWRK ), LWORK-IWRK+1, INFO )*ELSE** Compute eigenvalues only* (Workspace: need N+1, prefer N+HSWORK (see comments) )*IWRK = ITAUCALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,$ WORK( IWRK ), LWORK-IWRK+1, INFO )END IF** If INFO > 0 from DHSEQR, then quit*IF( INFO.GT.0 )$ GO TO 50*IF( WANTVL .OR. WANTVR ) THEN** Compute left and/or right eigenvectors* (Workspace: need 4*N)*CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,$ N, NOUT, WORK( IWRK ), IERR )END IF*IF( WANTVL ) THEN** Undo balancing of left eigenvectors* (Workspace: need N)*CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,$ IERR )** Normalize left eigenvectors and make largest component real*DO 20 I = 1, NIF( WI( I ).EQ.ZERO ) THENSCL = ONE / DNRM2( N, VL( 1, I ), 1 )CALL DSCAL( N, SCL, VL( 1, I ), 1 )ELSE IF( WI( I ).GT.ZERO ) THENSCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),$ DNRM2( N, VL( 1, I+1 ), 1 ) )CALL DSCAL( N, SCL, VL( 1, I ), 1 )CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )DO 10 K = 1, NWORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**210 CONTINUEK = IDAMAX( N, WORK( IWRK ), 1 )CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )VL( K, I+1 ) = ZEROEND IF20 CONTINUEEND IF*IF( WANTVR ) THEN** Undo balancing of right eigenvectors* (Workspace: need N)*CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,$ IERR )** Normalize right eigenvectors and make largest component real*DO 40 I = 1, NIF( WI( I ).EQ.ZERO ) THENSCL = ONE / DNRM2( N, VR( 1, I ), 1 )CALL DSCAL( N, SCL, VR( 1, I ), 1 )ELSE IF( WI( I ).GT.ZERO ) THENSCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),$ DNRM2( N, VR( 1, I+1 ), 1 ) )CALL DSCAL( N, SCL, VR( 1, I ), 1 )CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )DO 30 K = 1, NWORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**230 CONTINUEK = IDAMAX( N, WORK( IWRK ), 1 )CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )VR( K, I+1 ) = ZEROEND IF40 CONTINUEEND IF** Undo scaling if necessary*50 CONTINUEIF( SCALEA ) THENCALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),$ MAX( N-INFO, 1 ), IERR )CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),$ MAX( N-INFO, 1 ), IERR )IF( INFO.GT.0 ) THENCALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,$ IERR )CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,$ IERR )END IFEND IF*WORK( 1 ) = MAXWRKRETURN** End of DGEEV*ENDDOUBLE PRECISION FUNCTION RLANGE( NORM, M, N, A, LDA, WORK )** -- LAPACK auxiliary routine (version 3.0) --* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,* Courant Institute, Argonne National Lab, and Rice University* October 31, 1992** .. Scalar Arguments ..CHARACTER NORMINTEGER LDA, M, N* ..* .. Array Arguments ..DOUBLE PRECISION A( LDA, * ), WORK( * )* ..** Purpose* =======** DLANGE returns the value of the one norm, or the Frobenius norm, or* the infinity norm, or the element of largest absolute value of a* real matrix A.** Description* ===========** DLANGE returns the value** DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm'* (* ( norm1(A), NORM = '1', 'O' or 'o'* (* ( normI(A), NORM = 'I' or 'i'* (* ( normF(A), NORM = 'F', 'f', 'E' or 'e'** where norm1 denotes the one norm of a matrix (maximum column sum),* normI denotes the infinity norm of a matrix (maximum row sum) and* normF denotes the Frobenius norm of a matrix (square root of sum of* squares). Note that max(abs(A(i,j))) is not a matrix norm.** Arguments* =========** NORM (input) CHARACTER*1* Specifies the value to be returned in DLANGE as described* above.** M (input) INTEGER* The number of rows of the matrix A. M >= 0. When M = 0,* DLANGE is set to zero.** N (input) INTEGER* The number of columns of the matrix A. N >= 0. When N = 0,* DLANGE is set to zero.** A (input) DOUBLE PRECISION array, dimension (LDA,N)* The m by n matrix A.** LDA (input) INTEGER* The leading dimension of the array A. LDA >= max(M,1).** WORK (workspace) DOUBLE PRECISION array, dimension (LWORK),* where LWORK >= M when NORM = 'I'; otherwise, WORK is not* referenced.** =====================================================================** .. Parameters ..DOUBLE PRECISION ONE, ZEROPARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )* ..* .. Local Scalars ..INTEGER I, JDOUBLE PRECISION SCALE, SUM, VALUE* ..* .. External Subroutines ..EXTERNAL DLASSQ* ..* .. External Functions ..LOGICAL LSAMEREXTERNAL LSAMER* ..* .. Intrinsic Functions ..INTRINSIC ABS, MAX, MIN, SQRT* ..* .. Executable Statements ..*IF( MIN( M, N ).EQ.0 ) THENVALUE = ZEROELSE IF( LSAMER( NORM, 'M' ) ) THEN** Find max(abs(A(i,j))).*VALUE = ZERODO 20 J = 1, NDO 10 I = 1, MVALUE = MAX( VALUE, ABS( A( I, J ) ) )10 CONTINUE20 CONTINUEELSE IF( ( LSAMER( NORM, 'O' ) ) .OR. ( NORM.EQ.'1' ) ) THEN** Find norm1(A).*VALUE = ZERODO 40 J = 1, NSUM = ZERODO 30 I = 1, MSUM = SUM + ABS( A( I, J ) )30 CONTINUEVALUE = MAX( VALUE, SUM )40 CONTINUEELSE IF( LSAMER( NORM, 'I' ) ) THEN** Find normI(A).*DO 50 I = 1, MWORK( I ) = ZERO50 CONTINUEDO 70 J = 1, NDO 60 I = 1, MWORK( I ) = WORK( I ) + ABS( A( I, J ) )60 CONTINUE70 CONTINUEVALUE = ZERODO 80 I = 1, MVALUE = MAX( VALUE, WORK( I ) )80 CONTINUEELSE IF( (LSAMER( NORM, 'F' ) ) .OR. (LSAMER( NORM, 'E' ) ) ) THEN** Find normF(A).*SCALE = ZEROSUM = ONEDO 90 J = 1, NCALL DLASSQ( M, A( 1, J ), 1, SCALE, SUM )90 CONTINUEVALUE = SCALE*SQRT( SUM )END IF*RLANGE = VALUERETURN** End of DLANGEENDLOGICAL FUNCTION LSAMER( CA, CB )** -- LAPACK auxiliary routine (version 3.0) --* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,* Courant Institute, Argonne National Lab, and Rice University* September 30, 1994** .. Scalar Arguments ..CHARACTER CA, CB* ..** Purpose* =======** LSAME returns .TRUE. if CA is the same letter as CB regardless of* case.** Arguments* =========** CA (input) CHARACTER*1* CB (input) CHARACTER*1* CA and CB specify the single characters to be compared.** =====================================================================** .. Intrinsic Functions ..INTRINSIC ICHAR* ..* .. Local Scalars ..INTEGER INTA, INTB, ZCODE* ..* .. Executable Statements ..** Test if the characters are equal*LSAMER = CA.EQ.CBIF( LSAMER )$ RETURN** Now test for equivalence if both characters are alphabetic.*ZCODE = ICHAR( 'Z' )** Use 'Z' rather than 'A' so that ASCII can be detected on Prime* machines, on which ICHAR returns a value with bit 8 set.* ICHAR('A') on Prime machines returns 193 which is the same as* ICHAR('A') on an EBCDIC machine.*INTA = ICHAR( CA )INTB = ICHAR( CB )*IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN** ASCII is assumed - ZCODE is the ASCII code of either lower or* upper case 'Z'.*IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32*ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN** EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or* upper case 'Z'.*IF( INTA.GE.129 .AND. INTA.LE.137 .OR.$ INTA.GE.145 .AND. INTA.LE.153 .OR.$ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64IF( INTB.GE.129 .AND. INTB.LE.137 .OR.$ INTB.GE.145 .AND. INTB.LE.153 .OR.$ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64*ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN** ASCII is assumed, on Prime machines - ZCODE is the ASCII code* plus 128 of either lower or upper case 'Z'.*IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32END IFLSAMER = INTA.EQ.INTB** RETURN** End of LSAME*END