1DTGSNA(1) LAPACK routine (version 3.1) DTGSNA(1)
2
6 DTGSNA - 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 DTGSNA( 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DIF( * ), S( *
24 ), VL( LDVL, * ), VR( LDVR, * ), WORK( * )
25
27 DTGSNA 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 DGGES),
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DTGEVC. 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) DOUBLE PRECISION 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 DTGEVC. 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension
127 (MAX(1,LWORK))
128 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
129 LWORK (input) INTEGER
131 The dimension of the array WORK. LWORK >= max(1,N). If JOB =
132 'V' or 'B' LWORK >= 2*N*(N+2)+16.
133 If LWORK = -1, then a workspace query is assumed; the routine
135 only calculates the optimal size of the WORK array, returns
136 this value as the first entry of the WORK array, and no error
137 message related to LWORK is issued by XERBLA.
138 IWORK (workspace) INTEGER array, dimension (N + 6)
140 If JOB = 'E', IWORK is not referenced.
141 INFO (output) INTEGER
143 =0: Successful exit
144 <0: If INFO = -i, the i-th argument had an illegal value
145
147 The reciprocal of the condition number of a generalized eigenvalue w =
148 (a, b) is defined as
149 S(w) = (|u'Av|**2 + |u'Bv|**2)**(1/2) / (norm(u)*norm(v))
151 where u and v are the left and right eigenvectors of (A, B) correspond‐
153 ing to w; |z| denotes the absolute value of the complex number, and
154 norm(u) denotes the 2-norm of the vector u.
155 The pair (a, b) corresponds to an eigenvalue w = a/b (= u'Av/u'Bv) of
156 the matrix pair (A, B). If both a and b equal zero, then (A B) is sin‐
157 gular and S(I) = -1 is returned.
158 An approximate error bound on the chordal distance between the i-th
160 computed generalized eigenvalue w and the corresponding exact eigenval‐
161 ue lambda is
162 chord(w, lambda) <= EPS * norm(A, B) / S(I)
164 where EPS is the machine precision.
166 The reciprocal of the condition number DIF(i) of right eigenvector u
168 and left eigenvector v corresponding to the generalized eigenvalue w is
169 defined as follows:
170 a) If the i-th eigenvalue w = (a,b) is real
172 Suppose U and V are orthogonal transformations such that
174 U'*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
176 ( 0 S22 ),( 0 T22 ) n-1
177 1 n-1 1 n-1
178 Then the reciprocal condition number DIF(i) is
180 Difl((a, b), (S22, T22)) = sigma-min( Zl ),
182 where sigma-min(Zl) denotes the smallest singular value of the
184 2(n-1)-by-2(n-1) matrix
185 Zl = [ kron(a, In-1) -kron(1, S22) ]
187 [ kron(b, In-1) -kron(1, T22) ] .
188 Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
190 Kronecker product between the matrices X and Y.
191 Note that if the default method for computing DIF(i) is wanted
193 (see DLATDF), then the parameter DIFDRI (see below) should be
194 changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
195 See DTGSYL for more details.
196 b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
198 Suppose U and V are orthogonal transformations such that
200 U'*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
202 ( 0 S22 ),( 0 T22) n-2
203 2 n-2 2 n-2
204 and (S11, T11) corresponds to the complex conjugate eigenvalue
206 pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
207 that
208 U1'*S11*V1 = ( s11 s12 ) and U1'*T11*V1 = ( t11 t12 )
210 ( 0 s22 ) ( 0 t22 )
211 where the generalized eigenvalues w = s11/t11 and
213 conjg(w) = s22/t22.
214 Then the reciprocal condition number DIF(i) is bounded by
216 min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
218 where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
220 Z1 is the complex 2-by-2 matrix
221 Z1 = [ s11 -s22 ]
223 [ t11 -t22 ],
224 This is done by computing (using real arithmetic) the
226 roots of the characteristical polynomial det(Z1' * Z1 - lambda I),
227 where Z1' denotes the conjugate transpose of Z1 and det(X) denotes
228 the determinant of X.
229 and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
231 upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
232 Z2 = [ kron(S11', In-2) -kron(I2, S22) ]
234 [ kron(T11', In-2) -kron(I2, T22) ]
235 Note that if the default method for computing DIF is wanted (see
237 DLATDF), then the parameter DIFDRI (see below) should be changed
238 from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
239 for more details.
240 For each eigenvalue/vector specified by SELECT, DIF stores a Frobenius
242 norm-based estimate of Difl.
243 An approximate error bound for the i-th computed eigenvector VL(i) or
245 VR(i) is given by
246 EPS * norm(A, B) / DIF(i).
248 See ref. [2-3] for more details and further references.
250 Based on contributions by
252 Bo Kagstrom and Peter Poromaa, Department of Computing Science,
253 Umea University, S-901 87 Umea, Sweden.
254 References
256 ==========
257 [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
259 Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
260 M.S. Moonen et al (eds), Linear Algebra for Large Scale and
261 Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
262 [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
264 Eigenvalues of a Regular Matrix Pair (A, B) and Condition
265 Estimation: Theory, Algorithms and Software,
266 Report UMINF - 94.04, Department of Computing Science, Umea
267 University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
268 Note 87. To appear in Numerical Algorithms, 1996.
269 [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
271 for Solving the Generalized Sylvester Equation and Estimating the
272 Separation between Regular Matrix Pairs, Report UMINF - 93.23,
273 Department of Computing Science, Umea University, S-901 87 Umea,
274 Sweden, December 1993, Revised April 1994, Also as LAPACK Working
275 Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
276 No 1, 1996.
277 LAPACK routine (version 3.1) November 2006 DTGSNA(1)