1STGSNA(1) LAPACK routine (version 3.2) STGSNA(1)
2
6 STGSNA - estimates reciprocal condition numbers for specified eigenval‐
7 ues and/or eigenvectors of a matrix pair (A, B) in generalized real
8 Schur canonical form (or of any matrix pair (Q*A*Z', Q*B*Z') with
9 orthogonal matrices Q and Z, where Z' denotes the transpose of Z
10
12 SUBROUTINE STGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL,
13 VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO )
14 CHARACTER HOWMNY, JOB
16 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
18 LOGICAL SELECT( * )
20 INTEGER IWORK( * )
22 REAL A( LDA, * ), B( LDB, * ), DIF( * ), S( * ), VL(
24 LDVL, * ), VR( LDVR, * ), WORK( * )
25
27 STGSNA estimates reciprocal condition numbers for specified eigenvalues
28 and/or eigenvectors of a matrix pair (A, B) in generalized real Schur
29 canonical form (or of any matrix pair (Q*A*Z', Q*B*Z') with orthogonal
30 matrices Q and Z, where Z' denotes the transpose of Z. (A, B) must be
31 in generalized real Schur form (as returned by SGGES), i.e. A is block
32 upper triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper
33 triangular.
34
36 JOB (input) CHARACTER*1
37 Specifies whether condition numbers are required for eigenval‐
38 ues (S) or eigenvectors (DIF):
39 = 'E': for eigenvalues only (S);
40 = 'V': for eigenvectors only (DIF);
41 = 'B': for both eigenvalues and eigenvectors (S and DIF).
42 HOWMNY (input) CHARACTER*1
44 = 'A': compute condition numbers for all eigenpairs;
45 = 'S': compute condition numbers for selected eigenpairs speci‐
46 fied by the array SELECT.
47 SELECT (input) LOGICAL array, dimension (N)
49 If HOWMNY = 'S', SELECT specifies the eigenpairs for which con‐
50 dition numbers are required. To select condition numbers for
51 the eigenpair corresponding to a real eigenvalue w(j),
52 SELECT(j) must be set to .TRUE.. To select condition numbers
53 corresponding to a complex conjugate pair of eigenvalues w(j)
54 and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
55 set to .TRUE.. If HOWMNY = 'A', SELECT is not referenced.
56 N (input) INTEGER
58 The order of the square matrix pair (A, B). N >= 0.
59 A (input) REAL array, dimension (LDA,N)
61 The upper quasi-triangular matrix A in the pair (A,B).
62 LDA (input) INTEGER
64 The leading dimension of the array A. LDA >= max(1,N).
65 B (input) REAL array, dimension (LDB,N)
67 The upper triangular matrix B in the pair (A,B).
68 LDB (input) INTEGER
70 The leading dimension of the array B. LDB >= max(1,N).
71 VL (input) REAL array, dimension (LDVL,M)
73 If JOB = 'E' or 'B', VL must contain left eigenvectors of (A,
74 B), corresponding to the eigenpairs specified by HOWMNY and
75 SELECT. The eigenvectors must be stored in consecutive columns
76 of VL, as returned by STGEVC. If JOB = 'V', VL is not refer‐
77 enced.
78 LDVL (input) INTEGER
80 The leading dimension of the array VL. LDVL >= 1. If JOB = 'E'
81 or 'B', LDVL >= N.
82 VR (input) REAL array, dimension (LDVR,M)
84 If JOB = 'E' or 'B', VR must contain right eigenvectors of (A,
85 B), corresponding to the eigenpairs specified by HOWMNY and
86 SELECT. The eigenvectors must be stored in consecutive columns
87 ov VR, as returned by STGEVC. If JOB = 'V', VR is not refer‐
88 enced.
89 LDVR (input) INTEGER
91 The leading dimension of the array VR. LDVR >= 1. If JOB = 'E'
92 or 'B', LDVR >= N.
93 S (output) REAL array, dimension (MM)
95 If JOB = 'E' or 'B', the reciprocal condition numbers of the
96 selected eigenvalues, stored in consecutive elements of the
97 array. For a complex conjugate pair of eigenvalues two consecu‐
98 tive elements of S are set to the same value. Thus S(j),
99 DIF(j), and the j-th columns of VL and VR all correspond to the
100 same eigenpair (but not in general the j-th eigenpair, unless
101 all eigenpairs are selected). If JOB = 'V', S is not refer‐
102 enced.
103 DIF (output) REAL array, dimension (MM)
105 If JOB = 'V' or 'B', the estimated reciprocal condition numbers
106 of the selected eigenvectors, stored in consecutive elements of
107 the array. For a complex eigenvector two consecutive elements
108 of DIF are set to the same value. If the eigenvalues cannot be
109 reordered to compute DIF(j), DIF(j) is set to 0; this can only
110 occur when the true value would be very small anyway. If JOB =
111 'E', DIF is not referenced.
112 MM (input) INTEGER
114 The number of elements in the arrays S and DIF. MM >= M.
115 M (output) INTEGER
117 The number of elements of the arrays S and DIF used to store
118 the specified condition numbers; for each selected real eigen‐
119 value one element is used, and for each selected complex conju‐
120 gate pair of eigenvalues, two elements are used. If HOWMNY =
121 'A', M is set to N.
122 WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
124 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
125 LWORK (input) INTEGER
127 The dimension of the array WORK. LWORK >= max(1,N). If JOB =
128 'V' or 'B' LWORK >= 2*N*(N+2)+16. If LWORK = -1, then a
129 workspace query is assumed; the routine only calculates the
130 optimal size of the WORK array, returns this value as the first
131 entry of the WORK array, and no error message related to LWORK
132 is issued by XERBLA.
133 IWORK (workspace) INTEGER array, dimension (N + 6)
135 If JOB = 'E', IWORK is not referenced.
136 INFO (output) INTEGER
138 =0: Successful exit
139 <0: If INFO = -i, the i-th argument had an illegal value
140
142 The reciprocal of the condition number of a generalized eigenvalue w =
143 (a, b) is defined as
144 S(w) = (|u'Av|**2 + |u'Bv|**2)**(1/2) / (norm(u)*norm(v)) where u
145 and v are the left and right eigenvectors of (A, B) corresponding to w;
146 |z| denotes the absolute value of the complex number, and norm(u)
147 denotes the 2-norm of the vector u.
148 The pair (a, b) corresponds to an eigenvalue w = a/b (= u'Av/u'Bv) of
149 the matrix pair (A, B). If both a and b equal zero, then (A B) is sin‐
150 gular and S(I) = -1 is returned.
151 An approximate error bound on the chordal distance between the i-th
152 computed generalized eigenvalue w and the corresponding exact eigenval‐
153 ue lambda is
154 chord(w, lambda) <= EPS * norm(A, B) / S(I)
155 where EPS is the machine precision.
156 The reciprocal of the condition number DIF(i) of right eigenvector u
157 and left eigenvector v corresponding to the generalized eigenvalue w is
158 defined as follows:
159 a) If the i-th eigenvalue w = (a,b) is real
160 Suppose U and V are orthogonal transformations such that
161 U'*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
162 ( 0 S22 ),( 0 T22 ) n-1
163 1 n-1 1 n-1
164 Then the reciprocal condition number DIF(i) is
165 Difl((a, b), (S22, T22)) = sigma-min( Zl ),
166 where sigma-min(Zl) denotes the smallest singular value of the
167 2(n-1)-by-2(n-1) matrix
168 Zl = [ kron(a, In-1) -kron(1, S22) ]
169 [ kron(b, In-1) -kron(1, T22) ] .
170 Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
171 Kronecker product between the matrices X and Y.
172 Note that if the default method for computing DIF(i) is wanted
173 (see SLATDF), then the parameter DIFDRI (see below) should be
174 changed from 3 to 4 (routine SLATDF(IJOB = 2 will be used)).
175 See STGSYL for more details.
176 b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
177 Suppose U and V are orthogonal transformations such that
178 U'*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
179 ( 0 S22 ),( 0 T22) n-2
180 2 n-2 2 n-2
181 and (S11, T11) corresponds to the complex conjugate eigenvalue
182 pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
183 that
184 U1'*S11*V1 = ( s11 s12 ) and U1'*T11*V1 = ( t11 t12 )
185 ( 0 s22 ) ( 0 t22 )
186 where the generalized eigenvalues w = s11/t11 and
187 conjg(w) = s22/t22.
188 Then the reciprocal condition number DIF(i) is bounded by
189 min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
190 where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
191 Z1 is the complex 2-by-2 matrix
192 Z1 = [ s11 -s22 ]
193 [ t11 -t22 ],
194 This is done by computing (using real arithmetic) the
195 roots of the characteristical polynomial det(Z1' * Z1 - lambda I),
196 where Z1' denotes the conjugate transpose of Z1 and det(X) denotes
197 the determinant of X.
198 and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
199 upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
200 Z2 = [ kron(S11', In-2) -kron(I2, S22) ]
201 [ kron(T11', In-2) -kron(I2, T22) ]
202 Note that if the default method for computing DIF is wanted (see
203 SLATDF), then the parameter DIFDRI (see below) should be changed
204 from 3 to 4 (routine SLATDF(IJOB = 2 will be used)). See STGSYL
205 for more details.
206 For each eigenvalue/vector specified by SELECT, DIF stores a Frobenius
207 norm-based estimate of Difl.
208 An approximate error bound for the i-th computed eigenvector VL(i) or
209 VR(i) is given by
210 EPS * norm(A, B) / DIF(i).
211 See ref. [2-3] for more details and further references.
212 Based on contributions by
213 Bo Kagstrom and Peter Poromaa, Department of Computing Science,
214 Umea University, S-901 87 Umea, Sweden.
215 References
216 ==========
217 [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
218 Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
219 M.S. Moonen et al (eds), Linear Algebra for Large Scale and
220 Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
221 [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
222 Eigenvalues of a Regular Matrix Pair (A, B) and Condition
223 Estimation: Theory, Algorithms and Software,
224 Report UMINF - 94.04, Department of Computing Science, Umea
225 University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
226 Note 87. To appear in Numerical Algorithms, 1996.
227 [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
228 for Solving the Generalized Sylvester Equation and Estimating the
229 Separation between Regular Matrix Pairs, Report UMINF - 93.23,
230 Department of Computing Science, Umea University, S-901 87 Umea,
231 Sweden, December 1993, Revised April 1994, Also as LAPACK Working
232 Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
233 No 1, 1996.
234 LAPACK routine (version 3.2) November 2008 STGSNA(1)