1CTGSEN(1) LAPACK routine (version 3.1.1) CTGSEN(1)
2
6 CTGSEN - the generalized Schur decomposition of a complex matrix pair
7 (A, B) (in terms of an unitary equivalence trans- formation Q' * (A, B)
8 * Z), so that a selected cluster of eigenvalues appears in the leading
9 diagonal blocks of the pair (A,B)
10
12 SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
13 ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK,
14 LWORK, IWORK, LIWORK, INFO )
15 LOGICAL WANTQ, WANTZ
17 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, M, N
19 REAL PL, PR
21 LOGICAL SELECT( * )
23 INTEGER IWORK( * )
25 REAL DIF( * )
27 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), BETA( * ), Q(
29 LDQ, * ), WORK( * ), Z( LDZ, * )
30
32 CTGSEN reorders the generalized Schur decomposition of a complex matrix
33 pair (A, B) (in terms of an unitary equivalence trans- formation Q' *
34 (A, B) * Z), so that a selected cluster of eigenvalues appears in the
35 leading diagonal blocks of the pair (A,B). The leading columns of Q and
36 Z form unitary bases of the corresponding left and right eigenspaces
37 (deflating subspaces). (A, B) must be in generalized Schur canonical
38 form, that is, A and B are both upper triangular.
39 CTGSEN also computes the generalized eigenvalues
41 w(j)= ALPHA(j) / BETA(j)
43 of the reordered matrix pair (A, B).
45 Optionally, the routine computes estimates of reciprocal condition num‐
47 bers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
48 (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
49 between the matrix pairs (A11, B11) and (A22,B22) that correspond to
50 the selected cluster and the eigenvalues outside the cluster, resp.,
51 and norms of "projections" onto left and right eigenspaces w.r.t. the
52 selected cluster in the (1,1)-block.
53
57 IJOB (input) integer
58 Specifies whether condition numbers are required for the clus‐
59 ter of eigenvalues (PL and PR) or the deflating subspaces (Difu
60 and Difl):
61 =0: Only reorder w.r.t. SELECT. No extras.
62 =1: Reciprocal of norms of "projections" onto left and right
63 eigenspaces w.r.t. the selected cluster (PL and PR). =2: Upper
64 bounds on Difu and Difl. F-norm-based estimate
65 (DIF(1:2)).
66 =3: Estimate of Difu and Difl. 1-norm-based estimate
67 (DIF(1:2)). About 5 times as expensive as IJOB = 2. =4: Com‐
68 pute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic version
69 to get it all. =5: Compute PL, PR and DIF (i.e. 0, 1 and 3
70 above)
71 WANTQ (input) LOGICAL
73 WANTZ (input) LOGICAL
75 SELECT (input) LOGICAL array, dimension (N)
77 SELECT specifies the eigenvalues in the selected cluster. To
78 select an eigenvalue w(j), SELECT(j) must be set to
79 N (input) INTEGER
81 The order of the matrices A and B. N >= 0.
82 A (input/output) COMPLEX array, dimension(LDA,N)
84 On entry, the upper triangular matrix A, in generalized Schur
85 canonical form. On exit, A is overwritten by the reordered
86 matrix A.
87 LDA (input) INTEGER
89 The leading dimension of the array A. LDA >= max(1,N).
90 B (input/output) COMPLEX array, dimension(LDB,N)
92 On entry, the upper triangular matrix B, in generalized Schur
93 canonical form. On exit, B is overwritten by the reordered
94 matrix B.
95 LDB (input) INTEGER
97 The leading dimension of the array B. LDB >= max(1,N).
98 ALPHA (output) COMPLEX array, dimension (N)
100 BETA (output) COMPLEX array, dimension (N) The diagonal ele‐
101 ments of A and B, respectively, when the pair (A,B) has been
102 reduced to generalized Schur form. ALPHA(i)/BETA(i) i=1,...,N
103 are the generalized eigenvalues.
104 Q (input/output) COMPLEX array, dimension (LDQ,N)
106 On entry, if WANTQ = .TRUE., Q is an N-by-N matrix. On exit, Q
107 has been postmultiplied by the left unitary transformation
108 matrix which reorder (A, B); The leading M columns of Q form
109 orthonormal bases for the specified pair of left eigenspaces
110 (deflating subspaces). If WANTQ = .FALSE., Q is not refer‐
111 enced.
112 LDQ (input) INTEGER
114 The leading dimension of the array Q. LDQ >= 1. If WANTQ =
115 .TRUE., LDQ >= N.
116 Z (input/output) COMPLEX array, dimension (LDZ,N)
118 On entry, if WANTZ = .TRUE., Z is an N-by-N matrix. On exit, Z
119 has been postmultiplied by the left unitary transformation
120 matrix which reorder (A, B); The leading M columns of Z form
121 orthonormal bases for the specified pair of left eigenspaces
122 (deflating subspaces). If WANTZ = .FALSE., Z is not refer‐
123 enced.
124 LDZ (input) INTEGER
126 The leading dimension of the array Z. LDZ >= 1. If WANTZ =
127 .TRUE., LDZ >= N.
128 M (output) INTEGER
130 The dimension of the specified pair of left and right
131 eigenspaces, (deflating subspaces) 0 <= M <= N.
132 PL (output) REAL
134 PR (output) REAL If IJOB = 1, 4 or 5, PL, PR are lower
135 bounds on the reciprocal of the norm of "projections" onto left
136 and right eigenspace with respect to the selected cluster. 0 <
137 PL, PR <= 1. If M = 0 or M = N, PL = PR = 1. If IJOB = 0, 2 or
138 3 PL, PR are not referenced.
139 DIF (output) REAL array, dimension (2).
141 If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
142 If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
143 Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
144 estimates of Difu and Difl, computed using reversed communica‐
145 tion with CLACN2. If M = 0 or N, DIF(1:2) = F-norm([A, B]).
146 If IJOB = 0 or 1, DIF is not referenced.
147 WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
149 IF IJOB = 0, WORK is not referenced. Otherwise, on exit, if
150 INFO = 0, WORK(1) returns the optimal LWORK.
151 LWORK (input) INTEGER
153 The dimension of the array WORK. LWORK >= 1 If IJOB = 1, 2 or
154 4, LWORK >= 2*M*(N-M) If IJOB = 3 or 5, LWORK >= 4*M*(N-M)
155 If LWORK = -1, then a workspace query is assumed; the routine
157 only calculates the optimal size of the WORK array, returns
158 this value as the first entry of the WORK array, and no error
159 message related to LWORK is issued by XERBLA.
160 IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
162 IF IJOB = 0, IWORK is not referenced. Otherwise, on exit, if
163 INFO = 0, IWORK(1) returns the optimal LIWORK.
164 LIWORK (input) INTEGER
166 The dimension of the array IWORK. LIWORK >= 1. If IJOB = 1, 2
167 or 4, LIWORK >= N+2; If IJOB = 3 or 5, LIWORK >= MAX(N+2,
168 2*M*(N-M));
169 If LIWORK = -1, then a workspace query is assumed; the routine
171 only calculates the optimal size of the IWORK array, returns
172 this value as the first entry of the IWORK array, and no error
173 message related to LIWORK is issued by XERBLA.
174 INFO (output) INTEGER
176 =0: Successful exit.
177 <0: If INFO = -i, the i-th argument had an illegal value.
178 =1: Reordering of (A, B) failed because the transformed matrix
179 pair (A, B) would be too far from generalized Schur form; the
180 problem is very ill-conditioned. (A, B) may have been par‐
181 tially reordered. If requested, 0 is returned in DIF(*), PL
182 and PR.
183
185 CTGSEN first collects the selected eigenvalues by computing unitary U
186 and W that move them to the top left corner of (A, B). In other words,
187 the selected eigenvalues are the eigenvalues of (A11, B11) in
188 U'*(A, B)*W = (A11 A12) (B11 B12) n1
190 ( 0 A22),( 0 B22) n2
191 n1 n2 n1 n2
192 where N = n1+n2 and U' means the conjugate transpose of U. The first n1
194 columns of U and W span the specified pair of left and right
195 eigenspaces (deflating subspaces) of (A, B).
196 If (A, B) has been obtained from the generalized real Schur decomposi‐
198 tion of a matrix pair (C, D) = Q*(A, B)*Z', then the reordered general‐
199 ized Schur form of (C, D) is given by
200 (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
202 and the first n1 columns of Q*U and Z*W span the corresponding deflat‐
204 ing subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
205 Note that if the selected eigenvalue is sufficiently ill-conditioned,
207 then its value may differ significantly from its value before reorder‐
208 ing.
209 The reciprocal condition numbers of the left and right eigenspaces
211 spanned by the first n1 columns of U and W (or Q*U and Z*W) may be
212 returned in DIF(1:2), corresponding to Difu and Difl, resp.
213 The Difu and Difl are defined as:
215 Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
217 and
218 Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
219 where sigma-min(Zu) is the smallest singular value of the
221 (2*n1*n2)-by-(2*n1*n2) matrix
222 Zu = [ kron(In2, A11) -kron(A22', In1) ]
224 [ kron(In2, B11) -kron(B22', In1) ].
225 Here, Inx is the identity matrix of size nx and A22' is the transpose
227 of A22. kron(X, Y) is the Kronecker product between the matrices X and
228 Y.
229 When DIF(2) is small, small changes in (A, B) can cause large changes
231 in the deflating subspace. An approximate (asymptotic) bound on the
232 maximum angular error in the computed deflating subspaces is
233 EPS * norm((A, B)) / DIF(2),
235 where EPS is the machine precision.
237 The reciprocal norm of the projectors on the left and right eigenspaces
239 associated with (A11, B11) may be returned in PL and PR. They are com‐
240 puted as follows. First we compute L and R so that P*(A, B)*Q is block
241 diagonal, where
242 P = ( I -L ) n1 Q = ( I R ) n1
244 ( 0 I ) n2 and ( 0 I ) n2
245 n1 n2 n1 n2
246 and (L, R) is the solution to the generalized Sylvester equation
248 A11*R - L*A22 = -A12
250 B11*R - L*B22 = -B12
251 Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
253 An approximate (asymptotic) bound on the average absolute error of the
254 selected eigenvalues is
255 EPS * norm((A, B)) / PL.
257 There are also global error bounds which valid for perturbations up to
259 a certain restriction: A lower bound (x) on the smallest F-norm(E,F)
260 for which an eigenvalue of (A11, B11) may move and coalesce with an ei‐
261 genvalue of (A22, B22) under perturbation (E,F), (i.e. (A + E, B + F),
262 is
263 x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
265 An approximate bound on x can be computed from DIF(1:2), PL and PR.
267 If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed (L',
269 R') and unperturbed (L, R) left and right deflating subspaces associ‐
270 ated with the selected cluster in the (1,1)-blocks can be bounded as
271 max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
273 max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
274 See LAPACK User's Guide section 4.11 or the following references for
276 more information.
277 Note that if the default method for computing the Frobenius-norm- based
279 estimate DIF is not wanted (see CLATDF), then the parameter IDIFJB (see
280 below) should be changed from 3 to 4 (routine CLATDF (IJOB = 2 will be
281 used)). See CTGSYL for more details.
282 Based on contributions by
284 Bo Kagstrom and Peter Poromaa, Department of Computing Science,
285 Umea University, S-901 87 Umea, Sweden.
286 References
288 ==========
289 [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
291 Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
292 M.S. Moonen et al (eds), Linear Algebra for Large Scale and
293 Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
294 [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
296 Eigenvalues of a Regular Matrix Pair (A, B) and Condition
297 Estimation: Theory, Algorithms and Software, Report
298 UMINF - 94.04, Department of Computing Science, Umea University,
299 S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
300 To appear in Numerical Algorithms, 1996.
301 [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
303 for Solving the Generalized Sylvester Equation and Estimating the
304 Separation between Regular Matrix Pairs, Report UMINF - 93.23,
305 Department of Computing Science, Umea University, S-901 87 Umea,
306 Sweden, December 1993, Revised April 1994, Also as LAPACK working
307 Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
308 1996.
309 LAPACK routine (version 3.1.1) February 2007 CTGSEN(1)