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