dtgevc.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 1,038 行 · 第 1/5 页
HTML
1,038 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dtgevc.f</title>
<meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!--
body { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default a { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.string { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.string a { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.comment { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.comment a { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
--></style>
</head>
<body>
<pre>
SUBROUTINE <a name="DTGEVC.1"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a>( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
$ LDVL, VR, LDVR, MM, M, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment"> Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment"> November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> CHARACTER HOWMNY, SIDE
INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> LOGICAL SELECT( * )
DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
$ VR( LDVR, * ), WORK( * )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Purpose
</span><span class="comment">*</span><span class="comment"> =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DTGEVC.22"></a><a href="dtgevc.f.html#DTGEVC.1">DTGEVC</a> computes some or all of the right and/or left eigenvectors of
</span><span class="comment">*</span><span class="comment"> a pair of real matrices (S,P), where S is a quasi-triangular matrix
</span><span class="comment">*</span><span class="comment"> and P is upper triangular. Matrix pairs of this type are produced by
</span><span class="comment">*</span><span class="comment"> the generalized Schur factorization of a matrix pair (A,B):
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> A = Q*S*Z**T, B = Q*P*Z**T
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> as computed by <a name="DGGHRD.29"></a><a href="dgghrd.f.html#DGGHRD.1">DGGHRD</a> + <a name="DHGEQZ.29"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The right eigenvector x and the left eigenvector y of (S,P)
</span><span class="comment">*</span><span class="comment"> corresponding to an eigenvalue w are defined by:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> S*x = w*P*x, (y**H)*S = w*(y**H)*P,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> where y**H denotes the conjugate tranpose of y.
</span><span class="comment">*</span><span class="comment"> The eigenvalues are not input to this routine, but are computed
</span><span class="comment">*</span><span class="comment"> directly from the diagonal blocks of S and P.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> This routine returns the matrices X and/or Y of right and left
</span><span class="comment">*</span><span class="comment"> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
</span><span class="comment">*</span><span class="comment"> where Z and Q are input matrices.
</span><span class="comment">*</span><span class="comment"> If Q and Z are the orthogonal factors from the generalized Schur
</span><span class="comment">*</span><span class="comment"> factorization of a matrix pair (A,B), then Z*X and Q*Y
</span><span class="comment">*</span><span class="comment"> are the matrices of right and left eigenvectors of (A,B).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Arguments
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SIDE (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment"> = 'R': compute right eigenvectors only;
</span><span class="comment">*</span><span class="comment"> = 'L': compute left eigenvectors only;
</span><span class="comment">*</span><span class="comment"> = 'B': compute both right and left eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> HOWMNY (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment"> = 'A': compute all right and/or left eigenvectors;
</span><span class="comment">*</span><span class="comment"> = 'B': compute all right and/or left eigenvectors,
</span><span class="comment">*</span><span class="comment"> backtransformed by the matrices in VR and/or VL;
</span><span class="comment">*</span><span class="comment"> = 'S': compute selected right and/or left eigenvectors,
</span><span class="comment">*</span><span class="comment"> specified by the logical array SELECT.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SELECT (input) LOGICAL array, dimension (N)
</span><span class="comment">*</span><span class="comment"> If HOWMNY='S', SELECT specifies the eigenvectors to be
</span><span class="comment">*</span><span class="comment"> computed. If w(j) is a real eigenvalue, the corresponding
</span><span class="comment">*</span><span class="comment"> real eigenvector is computed if SELECT(j) is .TRUE..
</span><span class="comment">*</span><span class="comment"> If w(j) and w(j+1) are the real and imaginary parts of a
</span><span class="comment">*</span><span class="comment"> complex eigenvalue, the corresponding complex eigenvector
</span><span class="comment">*</span><span class="comment"> is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
</span><span class="comment">*</span><span class="comment"> and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
</span><span class="comment">*</span><span class="comment"> set to .FALSE..
</span><span class="comment">*</span><span class="comment"> Not referenced if HOWMNY = 'A' or 'B'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The order of the matrices S and P. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> S (input) DOUBLE PRECISION array, dimension (LDS,N)
</span><span class="comment">*</span><span class="comment"> The upper quasi-triangular matrix S from a generalized Schur
</span><span class="comment">*</span><span class="comment"> factorization, as computed by <a name="DHGEQZ.78"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDS (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of array S. LDS >= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> P (input) DOUBLE PRECISION array, dimension (LDP,N)
</span><span class="comment">*</span><span class="comment"> The upper triangular matrix P from a generalized Schur
</span><span class="comment">*</span><span class="comment"> factorization, as computed by <a name="DHGEQZ.85"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>.
</span><span class="comment">*</span><span class="comment"> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
</span><span class="comment">*</span><span class="comment"> of S must be in positive diagonal form.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDP (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of array P. LDP >= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
</span><span class="comment">*</span><span class="comment"> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
</span><span class="comment">*</span><span class="comment"> contain an N-by-N matrix Q (usually the orthogonal matrix Q
</span><span class="comment">*</span><span class="comment"> of left Schur vectors returned by <a name="DHGEQZ.95"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>).
</span><span class="comment">*</span><span class="comment"> On exit, if SIDE = 'L' or 'B', VL contains:
</span><span class="comment">*</span><span class="comment"> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
</span><span class="comment">*</span><span class="comment"> if HOWMNY = 'B', the matrix Q*Y;
</span><span class="comment">*</span><span class="comment"> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
</span><span class="comment">*</span><span class="comment"> SELECT, stored consecutively in the columns of
</span><span class="comment">*</span><span class="comment"> VL, in the same order as their eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> A complex eigenvector corresponding to a complex eigenvalue
</span><span class="comment">*</span><span class="comment"> is stored in two consecutive columns, the first holding the
</span><span class="comment">*</span><span class="comment"> real part, and the second the imaginary part.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Not referenced if SIDE = 'R'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDVL (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of array VL. LDVL >= 1, and if
</span><span class="comment">*</span><span class="comment"> SIDE = 'L' or 'B', LDVL >= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
</span><span class="comment">*</span><span class="comment"> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
</span><span class="comment">*</span><span class="comment"> contain an N-by-N matrix Z (usually the orthogonal matrix Z
</span><span class="comment">*</span><span class="comment"> of right Schur vectors returned by <a name="DHGEQZ.116"></a><a href="dhgeqz.f.html#DHGEQZ.1">DHGEQZ</a>).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> On exit, if SIDE = 'R' or 'B', VR contains:
</span><span class="comment">*</span><span class="comment"> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
</span><span class="comment">*</span><span class="comment"> if HOWMNY = 'B' or 'b', the matrix Z*X;
</span><span class="comment">*</span><span class="comment"> if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
</span><span class="comment">*</span><span class="comment"> specified by SELECT, stored consecutively in the
</span><span class="comment">*</span><span class="comment"> columns of VR, in the same order as their
</span><span class="comment">*</span><span class="comment"> eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> A complex eigenvector corresponding to a complex eigenvalue
</span><span class="comment">*</span><span class="comment"> is stored in two consecutive columns, the first holding the
</span><span class="comment">*</span><span class="comment"> real part and the second the imaginary part.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Not referenced if SIDE = 'L'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> LDVR (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The leading dimension of the array VR. LDVR >= 1, and if
</span><span class="comment">*</span><span class="comment"> SIDE = 'R' or 'B', LDVR >= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> MM (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of columns in the arrays VL and/or VR. MM >= M.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> M (output) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of columns in the arrays VL and/or VR actually
</span><span class="comment">*</span><span class="comment"> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
</span><span class="comment">*</span><span class="comment"> is set to N. Each selected real eigenvector occupies one
</span><span class="comment">*</span><span class="comment"> column and each selected complex eigenvector occupies two
</span><span class="comment">*</span><span class="comment"> columns.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) DOUBLE PRECISION array, dimension (6*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> = 0: successful exit.
</span><span class="comment">*</span><span class="comment"> < 0: if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment"> > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex
</span><span class="comment">*</span><span class="comment"> eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Further Details
</span><span class="comment">*</span><span class="comment"> ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Allocation of workspace:
</span><span class="comment">*</span><span class="comment"> ---------- -- ---------
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK( j ) = 1-norm of j-th column of A, above the diagonal
</span><span class="comment">*</span><span class="comment"> WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
</span><span class="comment">*</span><span class="comment"> WORK( 2*N+1:3*N ) = real part of eigenvector
</span><span class="comment">*</span><span class="comment"> WORK( 3*N+1:4*N ) = imaginary part of eigenvector
</span><span class="comment">*</span><span class="comment"> WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
</span><span class="comment">*</span><span class="comment"> WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Rowwise vs. columnwise solution methods:
</span><span class="comment">*</span><span class="comment"> ------- -- ---------- -------- -------
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Finding a generalized eigenvector consists basically of solving the
</span><span class="comment">*</span><span class="comment"> singular triangular system
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Consider finding the i-th right eigenvector (assume all eigenvalues
</span><span class="comment">*</span><span class="comment"> are real). The equation to be solved is:
</span><span class="comment">*</span><span class="comment"> n i
</span><span class="comment">*</span><span class="comment"> 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1
</span><span class="comment">*</span><span class="comment"> k=j k=j
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> where C = (A - w B) (The components v(i+1:n) are 0.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The "rowwise" method is:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> (1) v(i) := 1
</span><span class="comment">*</span><span class="comment"> for j = i-1,. . .,1:
</span><span class="comment">*</span><span class="comment"> i
</span><span class="comment">*</span><span class="comment"> (2) compute s = - sum C(j,k) v(k) and
</span><span class="comment">*</span><span class="comment"> k=j+1
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?