version 1.3, 2010/08/06 15:28:49
|
version 1.15, 2017/06/17 10:54:06
|
Line 1
|
Line 1
|
|
*> \brief \b DTREVC |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DTREVC + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, |
|
* LDVR, MM, M, WORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER HOWMNY, SIDE |
|
* INTEGER INFO, LDT, LDVL, LDVR, M, MM, N |
|
* .. |
|
* .. Array Arguments .. |
|
* LOGICAL SELECT( * ) |
|
* DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), |
|
* $ WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DTREVC computes some or all of the right and/or left eigenvectors of |
|
*> a real upper quasi-triangular matrix T. |
|
*> Matrices of this type are produced by the Schur factorization of |
|
*> a real general matrix: A = Q*T*Q**T, as computed by DHSEQR. |
|
*> |
|
*> The right eigenvector x and the left eigenvector y of T corresponding |
|
*> to an eigenvalue w are defined by: |
|
*> |
|
*> T*x = w*x, (y**H)*T = w*(y**H) |
|
*> |
|
*> where y**H denotes the conjugate transpose of y. |
|
*> The eigenvalues are not input to this routine, but are read directly |
|
*> from the diagonal blocks of T. |
|
*> |
|
*> This routine returns the matrices X and/or Y of right and left |
|
*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an |
|
*> input matrix. If Q is the orthogonal factor that reduces a matrix |
|
*> A to Schur form T, then Q*X and Q*Y are the matrices of right and |
|
*> left eigenvectors of A. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] SIDE |
|
*> \verbatim |
|
*> SIDE is CHARACTER*1 |
|
*> = 'R': compute right eigenvectors only; |
|
*> = 'L': compute left eigenvectors only; |
|
*> = 'B': compute both right and left eigenvectors. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] HOWMNY |
|
*> \verbatim |
|
*> HOWMNY is CHARACTER*1 |
|
*> = 'A': compute all right and/or left eigenvectors; |
|
*> = 'B': compute all right and/or left eigenvectors, |
|
*> backtransformed by the matrices in VR and/or VL; |
|
*> = 'S': compute selected right and/or left eigenvectors, |
|
*> as indicated by the logical array SELECT. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] SELECT |
|
*> \verbatim |
|
*> SELECT is LOGICAL array, dimension (N) |
|
*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be |
|
*> computed. |
|
*> If w(j) is a real eigenvalue, the corresponding real |
|
*> eigenvector is computed if SELECT(j) is .TRUE.. |
|
*> If w(j) and w(j+1) are the real and imaginary parts of a |
|
*> complex eigenvalue, the corresponding complex eigenvector is |
|
*> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and |
|
*> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to |
|
*> .FALSE.. |
|
*> Not referenced if HOWMNY = 'A' or 'B'. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The order of the matrix T. N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] T |
|
*> \verbatim |
|
*> T is DOUBLE PRECISION array, dimension (LDT,N) |
|
*> The upper quasi-triangular matrix T in Schur canonical form. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDT |
|
*> \verbatim |
|
*> LDT is INTEGER |
|
*> The leading dimension of the array T. LDT >= max(1,N). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] VL |
|
*> \verbatim |
|
*> VL is DOUBLE PRECISION array, dimension (LDVL,MM) |
|
*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must |
|
*> contain an N-by-N matrix Q (usually the orthogonal matrix Q |
|
*> of Schur vectors returned by DHSEQR). |
|
*> On exit, if SIDE = 'L' or 'B', VL contains: |
|
*> if HOWMNY = 'A', the matrix Y of left eigenvectors of T; |
|
*> if HOWMNY = 'B', the matrix Q*Y; |
|
*> if HOWMNY = 'S', the left eigenvectors of T specified by |
|
*> SELECT, stored consecutively in the columns |
|
*> of VL, in the same order as their |
|
*> eigenvalues. |
|
*> A complex eigenvector corresponding to a complex eigenvalue |
|
*> is stored in two consecutive columns, the first holding the |
|
*> real part, and the second the imaginary part. |
|
*> Not referenced if SIDE = 'R'. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDVL |
|
*> \verbatim |
|
*> LDVL is INTEGER |
|
*> The leading dimension of the array VL. LDVL >= 1, and if |
|
*> SIDE = 'L' or 'B', LDVL >= N. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] VR |
|
*> \verbatim |
|
*> VR is DOUBLE PRECISION array, dimension (LDVR,MM) |
|
*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must |
|
*> contain an N-by-N matrix Q (usually the orthogonal matrix Q |
|
*> of Schur vectors returned by DHSEQR). |
|
*> On exit, if SIDE = 'R' or 'B', VR contains: |
|
*> if HOWMNY = 'A', the matrix X of right eigenvectors of T; |
|
*> if HOWMNY = 'B', the matrix Q*X; |
|
*> if HOWMNY = 'S', the right eigenvectors of T specified by |
|
*> SELECT, stored consecutively in the columns |
|
*> of VR, in the same order as their |
|
*> eigenvalues. |
|
*> A complex eigenvector corresponding to a complex eigenvalue |
|
*> is stored in two consecutive columns, the first holding the |
|
*> real part and the second the imaginary part. |
|
*> Not referenced if SIDE = 'L'. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDVR |
|
*> \verbatim |
|
*> LDVR is INTEGER |
|
*> The leading dimension of the array VR. LDVR >= 1, and if |
|
*> SIDE = 'R' or 'B', LDVR >= N. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] MM |
|
*> \verbatim |
|
*> MM is INTEGER |
|
*> The number of columns in the arrays VL and/or VR. MM >= M. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The number of columns in the arrays VL and/or VR actually |
|
*> used to store the eigenvectors. |
|
*> If HOWMNY = 'A' or 'B', M is set to N. |
|
*> Each selected real eigenvector occupies one column and each |
|
*> selected complex eigenvector occupies two columns. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension (3*N) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date December 2016 |
|
* |
|
*> \ingroup doubleOTHERcomputational |
|
* |
|
*> \par Further Details: |
|
* ===================== |
|
*> |
|
*> \verbatim |
|
*> |
|
*> The algorithm used in this program is basically backward (forward) |
|
*> substitution, with scaling to make the the code robust against |
|
*> possible overflow. |
|
*> |
|
*> Each eigenvector is normalized so that the element of largest |
|
*> magnitude has magnitude 1; here the magnitude of a complex number |
|
*> (x,y) is taken to be |x| + |y|. |
|
*> \endverbatim |
|
*> |
|
* ===================================================================== |
SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, |
SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, |
$ LDVR, MM, M, WORK, INFO ) |
$ LDVR, MM, M, WORK, INFO ) |
* |
* |
* -- LAPACK routine (version 3.2) -- |
* -- LAPACK computational routine (version 3.7.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2006 |
* December 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER HOWMNY, SIDE |
CHARACTER HOWMNY, SIDE |
Line 16
|
Line 237
|
$ WORK( * ) |
$ WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DTREVC computes some or all of the right and/or left eigenvectors of |
|
* a real upper quasi-triangular matrix T. |
|
* Matrices of this type are produced by the Schur factorization of |
|
* a real general matrix: A = Q*T*Q**T, as computed by DHSEQR. |
|
* |
|
* The right eigenvector x and the left eigenvector y of T corresponding |
|
* to an eigenvalue w are defined by: |
|
* |
|
* T*x = w*x, (y**H)*T = w*(y**H) |
|
* |
|
* where y**H denotes the conjugate transpose of y. |
|
* The eigenvalues are not input to this routine, but are read directly |
|
* from the diagonal blocks of T. |
|
* |
|
* This routine returns the matrices X and/or Y of right and left |
|
* eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an |
|
* input matrix. If Q is the orthogonal factor that reduces a matrix |
|
* A to Schur form T, then Q*X and Q*Y are the matrices of right and |
|
* left eigenvectors of A. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* SIDE (input) CHARACTER*1 |
|
* = 'R': compute right eigenvectors only; |
|
* = 'L': compute left eigenvectors only; |
|
* = 'B': compute both right and left eigenvectors. |
|
* |
|
* HOWMNY (input) CHARACTER*1 |
|
* = 'A': compute all right and/or left eigenvectors; |
|
* = 'B': compute all right and/or left eigenvectors, |
|
* backtransformed by the matrices in VR and/or VL; |
|
* = 'S': compute selected right and/or left eigenvectors, |
|
* as indicated by the logical array SELECT. |
|
* |
|
* SELECT (input/output) LOGICAL array, dimension (N) |
|
* If HOWMNY = 'S', SELECT specifies the eigenvectors to be |
|
* computed. |
|
* If w(j) is a real eigenvalue, the corresponding real |
|
* eigenvector is computed if SELECT(j) is .TRUE.. |
|
* If w(j) and w(j+1) are the real and imaginary parts of a |
|
* complex eigenvalue, the corresponding complex eigenvector is |
|
* computed if either SELECT(j) or SELECT(j+1) is .TRUE., and |
|
* on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to |
|
* .FALSE.. |
|
* Not referenced if HOWMNY = 'A' or 'B'. |
|
* |
|
* N (input) INTEGER |
|
* The order of the matrix T. N >= 0. |
|
* |
|
* T (input) DOUBLE PRECISION array, dimension (LDT,N) |
|
* The upper quasi-triangular matrix T in Schur canonical form. |
|
* |
|
* LDT (input) INTEGER |
|
* The leading dimension of the array T. LDT >= max(1,N). |
|
* |
|
* VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) |
|
* On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must |
|
* contain an N-by-N matrix Q (usually the orthogonal matrix Q |
|
* of Schur vectors returned by DHSEQR). |
|
* On exit, if SIDE = 'L' or 'B', VL contains: |
|
* if HOWMNY = 'A', the matrix Y of left eigenvectors of T; |
|
* if HOWMNY = 'B', the matrix Q*Y; |
|
* if HOWMNY = 'S', the left eigenvectors of T specified by |
|
* SELECT, stored consecutively in the columns |
|
* of VL, in the same order as their |
|
* eigenvalues. |
|
* A complex eigenvector corresponding to a complex eigenvalue |
|
* is stored in two consecutive columns, the first holding the |
|
* real part, and the second the imaginary part. |
|
* Not referenced if SIDE = 'R'. |
|
* |
|
* LDVL (input) INTEGER |
|
* The leading dimension of the array VL. LDVL >= 1, and if |
|
* SIDE = 'L' or 'B', LDVL >= N. |
|
* |
|
* VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) |
|
* On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must |
|
* contain an N-by-N matrix Q (usually the orthogonal matrix Q |
|
* of Schur vectors returned by DHSEQR). |
|
* On exit, if SIDE = 'R' or 'B', VR contains: |
|
* if HOWMNY = 'A', the matrix X of right eigenvectors of T; |
|
* if HOWMNY = 'B', the matrix Q*X; |
|
* if HOWMNY = 'S', the right eigenvectors of T specified by |
|
* SELECT, stored consecutively in the columns |
|
* of VR, in the same order as their |
|
* eigenvalues. |
|
* A complex eigenvector corresponding to a complex eigenvalue |
|
* is stored in two consecutive columns, the first holding the |
|
* real part and the second the imaginary part. |
|
* Not referenced if SIDE = 'L'. |
|
* |
|
* LDVR (input) INTEGER |
|
* The leading dimension of the array VR. LDVR >= 1, and if |
|
* SIDE = 'R' or 'B', LDVR >= N. |
|
* |
|
* MM (input) INTEGER |
|
* The number of columns in the arrays VL and/or VR. MM >= M. |
|
* |
|
* M (output) INTEGER |
|
* The number of columns in the arrays VL and/or VR actually |
|
* used to store the eigenvectors. |
|
* If HOWMNY = 'A' or 'B', M is set to N. |
|
* Each selected real eigenvector occupies one column and each |
|
* selected complex eigenvector occupies two columns. |
|
* |
|
* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit |
|
* < 0: if INFO = -i, the i-th argument had an illegal value |
|
* |
|
* Further Details |
|
* =============== |
|
* |
|
* The algorithm used in this program is basically backward (forward) |
|
* substitution, with scaling to make the the code robust against |
|
* possible overflow. |
|
* |
|
* Each eigenvector is normalized so that the element of largest |
|
* magnitude has magnitude 1; here the magnitude of a complex number |
|
* (x,y) is taken to be |x| + |y|. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 651
|
Line 746
|
160 CONTINUE |
160 CONTINUE |
* |
* |
* Solve the quasi-triangular system: |
* Solve the quasi-triangular system: |
* (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK |
* (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK |
* |
* |
VMAX = ONE |
VMAX = ONE |
VCRIT = BIGNUM |
VCRIT = BIGNUM |
Line 688
|
Line 783
|
$ DDOT( J-KI-1, T( KI+1, J ), 1, |
$ DDOT( J-KI-1, T( KI+1, J ), 1, |
$ WORK( KI+1+N ), 1 ) |
$ WORK( KI+1+N ), 1 ) |
* |
* |
* Solve (T(J,J)-WR)'*X = WORK |
* Solve (T(J,J)-WR)**T*X = WORK |
* |
* |
CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), |
CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), |
$ LDT, ONE, ONE, WORK( J+N ), N, WR, |
$ LDT, ONE, ONE, WORK( J+N ), N, WR, |
Line 726
|
Line 821
|
$ WORK( KI+1+N ), 1 ) |
$ WORK( KI+1+N ), 1 ) |
* |
* |
* Solve |
* Solve |
* [T(J,J)-WR T(J,J+1) ]'* X = SCALE*( WORK1 ) |
* [T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 ) |
* [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) |
* [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) |
* |
* |
CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), |
CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), |
$ LDT, ONE, ONE, WORK( J+N ), N, WR, |
$ LDT, ONE, ONE, WORK( J+N ), N, WR, |
Line 778
|
Line 873
|
* Complex left eigenvector. |
* Complex left eigenvector. |
* |
* |
* Initial solve: |
* Initial solve: |
* ((T(KI,KI) T(KI,KI+1) )' - (WR - I* WI))*X = 0. |
* ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0. |
* ((T(KI+1,KI) T(KI+1,KI+1)) ) |
* ((T(KI+1,KI) T(KI+1,KI+1)) ) |
* |
* |
IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN |
IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN |
Line 891
|
Line 986
|
$ WORK( KI+2+N2 ), 1 ) |
$ WORK( KI+2+N2 ), 1 ) |
* |
* |
* Solve 2-by-2 complex linear equation |
* Solve 2-by-2 complex linear equation |
* ([T(j,j) T(j,j+1) ]'-(wr-i*wi)*I)*X = SCALE*B |
* ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B |
* ([T(j+1,j) T(j+1,j+1)] ) |
* ([T(j+1,j) T(j+1,j+1)] ) |
* |
* |
CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), |
CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), |
$ LDT, ONE, ONE, WORK( J+N ), N, WR, |
$ LDT, ONE, ONE, WORK( J+N ), N, WR, |