--- rpl/lapack/lapack/dgeevx.f 2010/12/21 13:53:25 1.7 +++ rpl/lapack/lapack/dgeevx.f 2023/08/07 08:38:48 1.20 @@ -1,11 +1,313 @@ +*> \brief DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGEEVX + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, +* VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, +* RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER BALANC, JOBVL, JOBVR, SENSE +* INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N +* DOUBLE PRECISION ABNRM +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ), +* $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ), +* $ WI( * ), WORK( * ), WR( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGEEVX computes for an N-by-N real nonsymmetric matrix A, the +*> eigenvalues and, optionally, the left and/or right eigenvectors. +*> +*> Optionally also, it computes a balancing transformation to improve +*> the conditioning of the eigenvalues and eigenvectors (ILO, IHI, +*> SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues +*> (RCONDE), and reciprocal condition numbers for the right +*> eigenvectors (RCONDV). +*> +*> 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. +*> +*> Balancing a matrix means permuting the rows and columns to make it +*> more nearly upper triangular, and applying a diagonal similarity +*> transformation D * A * D**(-1), where D is a diagonal matrix, to +*> make its rows and columns closer in norm and the condition numbers +*> of its eigenvalues and eigenvectors smaller. The computed +*> reciprocal condition numbers correspond to the balanced matrix. +*> Permuting rows and columns will not change the condition numbers +*> (in exact arithmetic) but diagonal scaling will. For further +*> explanation of balancing, see section 4.10.2 of the LAPACK +*> Users' Guide. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] BALANC +*> \verbatim +*> BALANC is CHARACTER*1 +*> Indicates how the input matrix should be diagonally scaled +*> and/or permuted to improve the conditioning of its +*> eigenvalues. +*> = 'N': Do not diagonally scale or permute; +*> = 'P': Perform permutations to make the matrix more nearly +*> upper triangular. Do not diagonally scale; +*> = 'S': Diagonally scale the matrix, i.e. replace A by +*> D*A*D**(-1), where D is a diagonal matrix chosen +*> to make the rows and columns of A more equal in +*> norm. Do not permute; +*> = 'B': Both diagonally scale and permute A. +*> +*> Computed reciprocal condition numbers will be for the matrix +*> after balancing and/or permuting. Permuting does not change +*> condition numbers (in exact arithmetic), but balancing does. +*> \endverbatim +*> +*> \param[in] JOBVL +*> \verbatim +*> JOBVL is CHARACTER*1 +*> = 'N': left eigenvectors of A are not computed; +*> = 'V': left eigenvectors of A are computed. +*> If SENSE = 'E' or 'B', JOBVL must = 'V'. +*> \endverbatim +*> +*> \param[in] JOBVR +*> \verbatim +*> JOBVR is CHARACTER*1 +*> = 'N': right eigenvectors of A are not computed; +*> = 'V': right eigenvectors of A are computed. +*> If SENSE = 'E' or 'B', JOBVR must = 'V'. +*> \endverbatim +*> +*> \param[in] SENSE +*> \verbatim +*> SENSE is CHARACTER*1 +*> Determines which reciprocal condition numbers are computed. +*> = 'N': None are computed; +*> = 'E': Computed for eigenvalues only; +*> = 'V': Computed for right eigenvectors only; +*> = 'B': Computed for eigenvalues and right eigenvectors. +*> +*> If SENSE = 'E' or 'B', both left and right eigenvectors +*> must also be computed (JOBVL = 'V' and JOBVR = 'V'). +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the N-by-N matrix A. +*> On exit, A has been overwritten. If JOBVL = 'V' or +*> JOBVR = 'V', A contains the real Schur form of the balanced +*> version of the input matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] WR +*> \verbatim +*> WR is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] WI +*> \verbatim +*> WI is DOUBLE PRECISION array, dimension (N) +*> WR and WI contain the real and imaginary parts, +*> respectively, of the computed eigenvalues. Complex +*> conjugate pairs of eigenvalues will appear consecutively +*> with the eigenvalue having the positive imaginary part +*> first. +*> \endverbatim +*> +*> \param[out] VL +*> \verbatim +*> VL is 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). +*> \endverbatim +*> +*> \param[in] LDVL +*> \verbatim +*> LDVL is INTEGER +*> The leading dimension of the array VL. LDVL >= 1; if +*> JOBVL = 'V', LDVL >= N. +*> \endverbatim +*> +*> \param[out] VR +*> \verbatim +*> VR is 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). +*> \endverbatim +*> +*> \param[in] LDVR +*> \verbatim +*> LDVR is INTEGER +*> The leading dimension of the array VR. LDVR >= 1, and if +*> JOBVR = 'V', LDVR >= N. +*> \endverbatim +*> +*> \param[out] ILO +*> \verbatim +*> ILO is INTEGER +*> \endverbatim +*> +*> \param[out] IHI +*> \verbatim +*> IHI is INTEGER +*> ILO and IHI are integer values determined when A was +*> balanced. The balanced A(i,j) = 0 if I > J and +*> J = 1,...,ILO-1 or I = IHI+1,...,N. +*> \endverbatim +*> +*> \param[out] SCALE +*> \verbatim +*> SCALE is DOUBLE PRECISION array, dimension (N) +*> Details of the permutations and scaling factors applied +*> when balancing A. If P(j) is the index of the row and column +*> interchanged with row and column j, and D(j) is the scaling +*> factor applied to row and column j, then +*> SCALE(J) = P(J), for J = 1,...,ILO-1 +*> = D(J), for J = ILO,...,IHI +*> = P(J) for J = IHI+1,...,N. +*> The order in which the interchanges are made is N to IHI+1, +*> then 1 to ILO-1. +*> \endverbatim +*> +*> \param[out] ABNRM +*> \verbatim +*> ABNRM is DOUBLE PRECISION +*> The one-norm of the balanced matrix (the maximum +*> of the sum of absolute values of elements of any column). +*> \endverbatim +*> +*> \param[out] RCONDE +*> \verbatim +*> RCONDE is DOUBLE PRECISION array, dimension (N) +*> RCONDE(j) is the reciprocal condition number of the j-th +*> eigenvalue. +*> \endverbatim +*> +*> \param[out] RCONDV +*> \verbatim +*> RCONDV is DOUBLE PRECISION array, dimension (N) +*> RCONDV(j) is the reciprocal condition number of the j-th +*> right eigenvector. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK. If SENSE = 'N' or 'E', +*> LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V', +*> LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6). +*> 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. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*N-2) +*> If SENSE = 'N' or 'E', not referenced. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is 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 or condition numbers +*> have been computed; elements 1:ILO-1 and i+1:N of WR +*> and WI contain eigenvalues which have converged. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +* +* @precisions fortran d -> s +* +*> \ingroup doubleGEeigen +* +* ===================================================================== SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) + implicit none * -* -- LAPACK driver routine (version 3.2) -- +* -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 * * .. Scalar Arguments .. CHARACTER BALANC, JOBVL, JOBVR, SENSE @@ -19,183 +321,6 @@ $ WI( * ), WORK( * ), WR( * ) * .. * -* Purpose -* ======= -* -* DGEEVX computes for an N-by-N real nonsymmetric matrix A, the -* eigenvalues and, optionally, the left and/or right eigenvectors. -* -* Optionally also, it computes a balancing transformation to improve -* the conditioning of the eigenvalues and eigenvectors (ILO, IHI, -* SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues -* (RCONDE), and reciprocal condition numbers for the right -* eigenvectors (RCONDV). -* -* 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. -* -* Balancing a matrix means permuting the rows and columns to make it -* more nearly upper triangular, and applying a diagonal similarity -* transformation D * A * D**(-1), where D is a diagonal matrix, to -* make its rows and columns closer in norm and the condition numbers -* of its eigenvalues and eigenvectors smaller. The computed -* reciprocal condition numbers correspond to the balanced matrix. -* Permuting rows and columns will not change the condition numbers -* (in exact arithmetic) but diagonal scaling will. For further -* explanation of balancing, see section 4.10.2 of the LAPACK -* Users' Guide. -* -* Arguments -* ========= -* -* BALANC (input) CHARACTER*1 -* Indicates how the input matrix should be diagonally scaled -* and/or permuted to improve the conditioning of its -* eigenvalues. -* = 'N': Do not diagonally scale or permute; -* = 'P': Perform permutations to make the matrix more nearly -* upper triangular. Do not diagonally scale; -* = 'S': Diagonally scale the matrix, i.e. replace A by -* D*A*D**(-1), where D is a diagonal matrix chosen -* to make the rows and columns of A more equal in -* norm. Do not permute; -* = 'B': Both diagonally scale and permute A. -* -* Computed reciprocal condition numbers will be for the matrix -* after balancing and/or permuting. Permuting does not change -* condition numbers (in exact arithmetic), but balancing does. -* -* JOBVL (input) CHARACTER*1 -* = 'N': left eigenvectors of A are not computed; -* = 'V': left eigenvectors of A are computed. -* If SENSE = 'E' or 'B', JOBVL must = 'V'. -* -* JOBVR (input) CHARACTER*1 -* = 'N': right eigenvectors of A are not computed; -* = 'V': right eigenvectors of A are computed. -* If SENSE = 'E' or 'B', JOBVR must = 'V'. -* -* SENSE (input) CHARACTER*1 -* Determines which reciprocal condition numbers are computed. -* = 'N': None are computed; -* = 'E': Computed for eigenvalues only; -* = 'V': Computed for right eigenvectors only; -* = 'B': Computed for eigenvalues and right eigenvectors. -* -* If SENSE = 'E' or 'B', both left and right eigenvectors -* must also be computed (JOBVL = 'V' and JOBVR = 'V'). -* -* 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. If JOBVL = 'V' or -* JOBVR = 'V', A contains the real Schur form of the balanced -* version of the input matrix A. -* -* 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 will 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, and if -* JOBVR = 'V', LDVR >= N. -* -* ILO (output) INTEGER -* IHI (output) INTEGER -* ILO and IHI are integer values determined when A was -* balanced. The balanced A(i,j) = 0 if I > J and -* J = 1,...,ILO-1 or I = IHI+1,...,N. -* -* SCALE (output) DOUBLE PRECISION array, dimension (N) -* Details of the permutations and scaling factors applied -* when balancing A. If P(j) is the index of the row and column -* interchanged with row and column j, and D(j) is the scaling -* factor applied to row and column j, then -* SCALE(J) = P(J), for J = 1,...,ILO-1 -* = D(J), for J = ILO,...,IHI -* = P(J) for J = IHI+1,...,N. -* The order in which the interchanges are made is N to IHI+1, -* then 1 to ILO-1. -* -* ABNRM (output) DOUBLE PRECISION -* The one-norm of the balanced matrix (the maximum -* of the sum of absolute values of elements of any column). -* -* RCONDE (output) DOUBLE PRECISION array, dimension (N) -* RCONDE(j) is the reciprocal condition number of the j-th -* eigenvalue. -* -* RCONDV (output) DOUBLE PRECISION array, dimension (N) -* RCONDV(j) is the reciprocal condition number of the j-th -* right eigenvector. -* -* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) -* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. -* -* LWORK (input) INTEGER -* The dimension of the array WORK. If SENSE = 'N' or 'E', -* LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V', -* LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6). -* 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. -* -* IWORK (workspace) INTEGER array, dimension (2*N-2) -* If SENSE = 'N' or 'E', not referenced. -* -* 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 or condition numbers -* have been computed; elements 1:ILO-1 and i+1:N of WR -* and WI contain eigenvalues which have converged. -* * ===================================================================== * * .. Parameters .. @@ -206,8 +331,8 @@ LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE, $ WNTSNN, WNTSNV CHARACTER JOB, SIDE - INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK, - $ MINWRK, NOUT + INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, + $ LWORK_TREVC, MAXWRK, MINWRK, NOUT DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, $ SN * .. @@ -217,7 +342,7 @@ * .. * .. External Subroutines .. EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY, - $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC, + $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3, $ DTRSNA, XERBLA * .. * .. External Functions .. @@ -242,8 +367,8 @@ WNTSNE = LSAME( SENSE, 'E' ) WNTSNV = LSAME( SENSE, 'V' ) WNTSNB = LSAME( SENSE, 'B' ) - IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, - $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) + IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) + $ .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) $ THEN INFO = -1 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN @@ -282,9 +407,19 @@ MAXWRK = N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 ) * IF( WANTVL ) THEN + CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA, + $ VL, LDVL, VR, LDVR, + $ N, NOUT, WORK, -1, IERR ) + LWORK_TREVC = INT( WORK(1) ) + MAXWRK = MAX( MAXWRK, N + LWORK_TREVC ) CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL, $ WORK, -1, INFO ) ELSE IF( WANTVR ) THEN + CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA, + $ VL, LDVL, VR, LDVR, + $ N, NOUT, WORK, -1, IERR ) + LWORK_TREVC = INT( WORK(1) ) + MAXWRK = MAX( MAXWRK, N + LWORK_TREVC ) CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR, $ WORK, -1, INFO ) ELSE @@ -296,7 +431,7 @@ $ LDVR, WORK, -1, INFO ) END IF END IF - HSWORK = WORK( 1 ) + HSWORK = INT( WORK(1) ) * IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN MINWRK = 2*N @@ -448,18 +583,18 @@ $ WORK( IWRK ), LWORK-IWRK+1, INFO ) END IF * -* If INFO > 0 from DHSEQR, then quit +* If INFO .NE. 0 from DHSEQR, then quit * - IF( INFO.GT.0 ) + IF( INFO.NE.0 ) $ GO TO 50 * IF( WANTVL .OR. WANTVR ) THEN * * Compute left and/or right eigenvectors -* (Workspace: need 3*N) +* (Workspace: need 3*N, prefer N + 2*N*NB) * - CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, - $ N, NOUT, WORK( IWRK ), IERR ) + CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, + $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR ) END IF * * Compute condition numbers if desired