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