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