1DTGSEN(1) LAPACK routine (version 3.1.1) DTGSEN(1)
2
6 DTGSEN - the generalized real Schur decomposition of a real matrix pair
7 (A, B) (in terms of an orthonormal equivalence trans- formation Q' *
8 (A, B) * Z), so that a selected cluster of eigenvalues appears in the
9 leading diagonal blocks of the upper quasi-triangular matrix A and the
10 upper triangular B
11
13 SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
14 ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR,
15 DIF, WORK, LWORK, IWORK, LIWORK, INFO )
16 LOGICAL WANTQ, WANTZ
18 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, M, N
20 DOUBLE PRECISION PL, PR
22 LOGICAL SELECT( * )
24 INTEGER IWORK( * )
26 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), B(
28 LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), WORK( *
29 ), Z( LDZ, * )
30
32 DTGSEN reorders the generalized real Schur decomposition of a real
33 matrix pair (A, B) (in terms of an orthonormal equivalence trans- for‐
34 mation Q' * (A, B) * Z), so that a selected cluster of eigenvalues
35 appears in the leading diagonal blocks of the upper quasi-triangular
36 matrix A and the upper triangular B. The leading columns of Q and Z
37 form orthonormal bases of the corresponding left and right eigen- spa‐
38 ces (deflating subspaces). (A, B) must be in generalized real Schur
39 canonical form (as returned by DGGES), i.e. A is block upper triangular
40 with 1-by-1 and 2-by-2 diagonal blocks. B is upper triangular.
41 DTGSEN also computes the generalized eigenvalues
43 w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
45 of the reordered matrix pair (A, B).
47 Optionally, DTGSEN computes the estimates of reciprocal condition num‐
49 bers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
50 (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
51 between the matrix pairs (A11, B11) and (A22,B22) that correspond to
52 the selected cluster and the eigenvalues outside the cluster, resp.,
53 and norms of "projections" onto left and right eigenspaces w.r.t. the
54 selected cluster in the (1,1)-block.
55
58 IJOB (input) INTEGER
59 Specifies whether condition numbers are required for the clus‐
60 ter of eigenvalues (PL and PR) or the deflating subspaces (Difu
61 and Difl):
62 =0: Only reorder w.r.t. SELECT. No extras.
63 =1: Reciprocal of norms of "projections" onto left and right
64 eigenspaces w.r.t. the selected cluster (PL and PR). =2: Upper
65 bounds on Difu and Difl. F-norm-based estimate
66 (DIF(1:2)).
67 =3: Estimate of Difu and Difl. 1-norm-based estimate
68 (DIF(1:2)). About 5 times as expensive as IJOB = 2. =4: Com‐
69 pute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic version
70 to get it all. =5: Compute PL, PR and DIF (i.e. 0, 1 and 3
71 above)
72 WANTQ (input) LOGICAL
74 WANTZ (input) LOGICAL
76 SELECT (input) LOGICAL array, dimension (N)
78 SELECT specifies the eigenvalues in the selected cluster. To
79 select a real eigenvalue w(j), SELECT(j) must be set to w(j)
80 and w(j+1), corresponding to a 2-by-2 diagonal block, either
81 SELECT(j) or SELECT(j+1) or both must be set to either both
82 included in the cluster or both excluded.
83 N (input) INTEGER
85 The order of the matrices A and B. N >= 0.
86 A (input/output) DOUBLE PRECISION array, dimension(LDA,N)
88 On entry, the upper quasi-triangular matrix A, with (A, B) in
89 generalized real Schur canonical form. On exit, A is overwrit‐
90 ten by the reordered matrix A.
91 LDA (input) INTEGER
93 The leading dimension of the array A. LDA >= max(1,N).
94 B (input/output) DOUBLE PRECISION array, dimension(LDB,N)
96 On entry, the upper triangular matrix B, with (A, B) in gener‐
97 alized real Schur canonical form. On exit, B is overwritten by
98 the reordered matrix B.
99 LDB (input) INTEGER
101 The leading dimension of the array B. LDB >= max(1,N).
102 ALPHAR (output) DOUBLE PRECISION array, dimension (N)
104 ALPHAI (output) DOUBLE PRECISION array, dimension (N) BETA
105 (output) DOUBLE PRECISION array, dimension (N) On exit,
106 (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will be the gen‐
107 eralized eigenvalues. ALPHAR(j) + ALPHAI(j)*i and
108 BETA(j),j=1,...,N are the diagonals of the complex Schur form
109 (S,T) that would result if the 2-by-2 diagonal blocks of the
110 real generalized Schur form of (A,B) were further reduced to
111 triangular form using complex unitary transformations. If
112 ALPHAI(j) is zero, then the j-th eigenvalue is real; if posi‐
113 tive, then the j-th and (j+1)-st eigenvalues are a complex con‐
114 jugate pair, with ALPHAI(j+1) negative.
115 Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
117 On entry, if WANTQ = .TRUE., Q is an N-by-N matrix. On exit, Q
118 has been postmultiplied by the left orthogonal transformation
119 matrix which reorder (A, B); The leading M columns of Q form
120 orthonormal bases for the specified pair of left eigenspaces
121 (deflating subspaces). If WANTQ = .FALSE., Q is not refer‐
122 enced.
123 LDQ (input) INTEGER
125 The leading dimension of the array Q. LDQ >= 1; and if WANTQ =
126 .TRUE., LDQ >= N.
127 Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
129 On entry, if WANTZ = .TRUE., Z is an N-by-N matrix. On exit, Z
130 has been postmultiplied by the left orthogonal transformation
131 matrix which reorder (A, B); The leading M columns of Z form
132 orthonormal bases for the specified pair of left eigenspaces
133 (deflating subspaces). If WANTZ = .FALSE., Z is not refer‐
134 enced.
135 LDZ (input) INTEGER
137 The leading dimension of the array Z. LDZ >= 1; If WANTZ =
138 .TRUE., LDZ >= N.
139 M (output) INTEGER
141 The dimension of the specified pair of left and right eigen-
142 spaces (deflating subspaces). 0 <= M <= N.
143 PL (output) DOUBLE PRECISION
145 PR (output) DOUBLE PRECISION If IJOB = 1, 4 or 5, PL, PR
146 are lower bounds on the reciprocal of the norm of "projections"
147 onto left and right eigenspaces with respect to the selected
148 cluster. 0 < PL, PR <= 1. If M = 0 or M = N, PL = PR = 1.
149 If IJOB = 0, 2 or 3, PL and PR are not referenced.
150 DIF (output) DOUBLE PRECISION array, dimension (2).
152 If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
153 If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
154 Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
155 estimates of Difu and Difl. If M = 0 or N, DIF(1:2) = F-
156 norm([A, B]). If IJOB = 0 or 1, DIF is not referenced.
157 WORK (workspace/output) DOUBLE PRECISION array,
159 dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns
160 the optimal LWORK.
161 LWORK (input) INTEGER
163 The dimension of the array WORK. LWORK >= 4*N+16. If IJOB =
164 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)). If IJOB = 3 or 5,
165 LWORK >= MAX(4*N+16, 4*M*(N-M)).
166 If LWORK = -1, then a workspace query is assumed; the routine
168 only calculates the optimal size of the WORK array, returns
169 this value as the first entry of the WORK array, and no error
170 message related to LWORK is issued by XERBLA.
171 IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
173 IF IJOB = 0, IWORK is not referenced. Otherwise, on exit, if
174 INFO = 0, IWORK(1) returns the optimal LIWORK.
175 LIWORK (input) INTEGER
177 The dimension of the array IWORK. LIWORK >= 1. If IJOB = 1, 2
178 or 4, LIWORK >= N+6. If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-
179 M), N+6).
180 If LIWORK = -1, then a workspace query is assumed; the routine
182 only calculates the optimal size of the IWORK array, returns
183 this value as the first entry of the IWORK array, and no error
184 message related to LIWORK is issued by XERBLA.
185 INFO (output) INTEGER
187 =0: Successful exit.
188 <0: If INFO = -i, the i-th argument had an illegal value.
189 =1: Reordering of (A, B) failed because the transformed matrix
190 pair (A, B) would be too far from generalized Schur form; the
191 problem is very ill-conditioned. (A, B) may have been par‐
192 tially reordered. If requested, 0 is returned in DIF(*), PL
193 and PR.
194
196 DTGSEN first collects the selected eigenvalues by computing orthogonal
197 U and W that move them to the top left corner of (A, B). In other
198 words, the selected eigenvalues are the eigenvalues of (A11, B11) in:
199 U'*(A, B)*W = (A11 A12) (B11 B12) n1
201 ( 0 A22),( 0 B22) n2
202 n1 n2 n1 n2
203 where N = n1+n2 and U' means the transpose of U. The first n1 columns
205 of U and W span the specified pair of left and right eigenspaces
206 (deflating subspaces) of (A, B).
207 If (A, B) has been obtained from the generalized real Schur decomposi‐
209 tion of a matrix pair (C, D) = Q*(A, B)*Z', then the reordered general‐
210 ized real Schur form of (C, D) is given by
211 (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
213 and the first n1 columns of Q*U and Z*W span the corresponding deflat‐
215 ing subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
216 Note that if the selected eigenvalue is sufficiently ill-conditioned,
218 then its value may differ significantly from its value before reorder‐
219 ing.
220 The reciprocal condition numbers of the left and right eigenspaces
222 spanned by the first n1 columns of U and W (or Q*U and Z*W) may be
223 returned in DIF(1:2), corresponding to Difu and Difl, resp.
224 The Difu and Difl are defined as:
226 Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
228 and
229 Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
230 where sigma-min(Zu) is the smallest singular value of the
232 (2*n1*n2)-by-(2*n1*n2) matrix
233 Zu = [ kron(In2, A11) -kron(A22', In1) ]
235 [ kron(In2, B11) -kron(B22', In1) ].
236 Here, Inx is the identity matrix of size nx and A22' is the transpose
238 of A22. kron(X, Y) is the Kronecker product between the matrices X and
239 Y.
240 When DIF(2) is small, small changes in (A, B) can cause large changes
242 in the deflating subspace. An approximate (asymptotic) bound on the
243 maximum angular error in the computed deflating subspaces is
244 EPS * norm((A, B)) / DIF(2),
246 where EPS is the machine precision.
248 The reciprocal norm of the projectors on the left and right eigenspaces
250 associated with (A11, B11) may be returned in PL and PR. They are com‐
251 puted as follows. First we compute L and R so that P*(A, B)*Q is block
252 diagonal, where
253 P = ( I -L ) n1 Q = ( I R ) n1
255 ( 0 I ) n2 and ( 0 I ) n2
256 n1 n2 n1 n2
257 and (L, R) is the solution to the generalized Sylvester equation
259 A11*R - L*A22 = -A12
261 B11*R - L*B22 = -B12
262 Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
264 An approximate (asymptotic) bound on the average absolute error of the
265 selected eigenvalues is
266 EPS * norm((A, B)) / PL.
268 There are also global error bounds which valid for perturbations up to
270 a certain restriction: A lower bound (x) on the smallest F-norm(E,F)
271 for which an eigenvalue of (A11, B11) may move and coalesce with an ei‐
272 genvalue of (A22, B22) under perturbation (E,F), (i.e. (A + E, B + F),
273 is
274 x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
276 An approximate bound on x can be computed from DIF(1:2), PL and PR.
278 If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed (L',
280 R') and unperturbed (L, R) left and right deflating subspaces associ‐
281 ated with the selected cluster in the (1,1)-blocks can be bounded as
282 max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
284 max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
285 See LAPACK User's Guide section 4.11 or the following references for
287 more information.
288 Note that if the default method for computing the Frobenius-norm- based
290 estimate DIF is not wanted (see DLATDF), then the parameter IDIFJB (see
291 below) should be changed from 3 to 4 (routine DLATDF (IJOB = 2 will be
292 used)). See DTGSYL for more details.
293 Based on contributions by
295 Bo Kagstrom and Peter Poromaa, Department of Computing Science,
296 Umea University, S-901 87 Umea, Sweden.
297 References
299 ==========
300 [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
302 Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
303 M.S. Moonen et al (eds), Linear Algebra for Large Scale and
304 Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
305 [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
307 Eigenvalues of a Regular Matrix Pair (A, B) and Condition
308 Estimation: Theory, Algorithms and Software,
309 Report UMINF - 94.04, Department of Computing Science, Umea
310 University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
311 Note 87. To appear in Numerical Algorithms, 1996.
312 [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
314 for Solving the Generalized Sylvester Equation and Estimating the
315 Separation between Regular Matrix Pairs, Report UMINF - 93.23,
316 Department of Computing Science, Umea University, S-901 87 Umea,
317 Sweden, December 1993, Revised April 1994, Also as LAPACK Working
318 Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
319 1996.
320 LAPACK routine (version 3.1.1) February 2007 DTGSEN(1)