1STRSEN(1) LAPACK routine (version 3.1) STRSEN(1)
2
6 STRSEN - the real Schur factorization of a real matrix A = Q*T*Q**T, so
7 that a selected cluster of eigenvalues appears in the leading diagonal
8 blocks of the upper quasi-triangular matrix T,
9
11 SUBROUTINE STRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S,
12 SEP, WORK, LWORK, IWORK, LIWORK, INFO )
13 CHARACTER COMPQ, JOB
15 INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
17 REAL S, SEP
19 LOGICAL SELECT( * )
21 INTEGER IWORK( * )
23 REAL Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ), WR( *
25 )
26
28 STRSEN reorders the real Schur factorization of a real matrix A =
29 Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
30 leading diagonal blocks of the upper quasi-triangular matrix T, and the
31 leading columns of Q form an orthonormal basis of the corresponding
32 right invariant subspace.
33 Optionally the routine computes the reciprocal condition numbers of the
35 cluster of eigenvalues and/or the invariant subspace.
36 T must be in Schur canonical form (as returned by SHSEQR), that is,
38 block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
39 2-by-2 diagonal block has its diagonal elemnts equal and its off-diago‐
40 nal elements of opposite sign.
41
44 JOB (input) CHARACTER*1
45 Specifies whether condition numbers are required for the clus‐
46 ter of eigenvalues (S) or the invariant subspace (SEP):
47 = 'N': none;
48 = 'E': for eigenvalues only (S);
49 = 'V': for invariant subspace only (SEP);
50 = 'B': for both eigenvalues and invariant subspace (S and SEP).
51 COMPQ (input) CHARACTER*1
53 = 'V': update the matrix Q of Schur vectors;
54 = 'N': do not update Q.
55 SELECT (input) LOGICAL array, dimension (N)
57 SELECT specifies the eigenvalues in the selected cluster. To
58 select a real eigenvalue w(j), SELECT(j) must be set to w(j)
59 and w(j+1), corresponding to a 2-by-2 diagonal block, either
60 SELECT(j) or SELECT(j+1) or both must be set to either both
61 included in the cluster or both excluded.
62 N (input) INTEGER
64 The order of the matrix T. N >= 0.
65 T (input/output) REAL array, dimension (LDT,N)
67 On entry, the upper quasi-triangular matrix T, in Schur canoni‐
68 cal form. On exit, T is overwritten by the reordered matrix T,
69 again in Schur canonical form, with the selected eigenvalues in
70 the leading diagonal blocks.
71 LDT (input) INTEGER
73 The leading dimension of the array T. LDT >= max(1,N).
74 Q (input/output) REAL array, dimension (LDQ,N)
76 On entry, if COMPQ = 'V', the matrix Q of Schur vectors. On
77 exit, if COMPQ = 'V', Q has been postmultiplied by the orthogo‐
78 nal transformation matrix which reorders T; the leading M col‐
79 umns of Q form an orthonormal basis for the specified invariant
80 subspace. If COMPQ = 'N', Q is not referenced.
81 LDQ (input) INTEGER
83 The leading dimension of the array Q. LDQ >= 1; and if COMPQ =
84 'V', LDQ >= N.
85 WR (output) REAL array, dimension (N)
87 WI (output) REAL array, dimension (N) The real and imagi‐
88 nary parts, respectively, of the reordered eigenvalues of T.
89 The eigenvalues are stored in the same order as on the diagonal
90 of T, with WR(i) = T(i,i) and, if T(i:i+1,i:i+1) is a 2-by-2
91 diagonal block, WI(i) > 0 and WI(i+1) = -WI(i). Note that if a
92 complex eigenvalue is sufficiently ill-conditioned, then its
93 value may differ significantly from its value before reorder‐
94 ing.
95 M (output) INTEGER
97 The dimension of the specified invariant subspace. 0 < = M <=
98 N.
99 S (output) REAL
101 If JOB = 'E' or 'B', S is a lower bound on the reciprocal con‐
102 dition number for the selected cluster of eigenvalues. S can‐
103 not underestimate the true reciprocal condition number by more
104 than a factor of sqrt(N). If M = 0 or N, S = 1. If JOB = 'N'
105 or 'V', S is not referenced.
106 SEP (output) REAL
108 If JOB = 'V' or 'B', SEP is the estimated reciprocal condition
109 number of the specified invariant subspace. If M = 0 or N, SEP
110 = norm(T). If JOB = 'N' or 'E', SEP is not referenced.
111 WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
113 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
114 LWORK (input) INTEGER
116 The dimension of the array WORK. If JOB = 'N', LWORK >=
117 max(1,N); if JOB = 'E', LWORK >= max(1,M*(N-M)); if JOB = 'V'
118 or 'B', LWORK >= max(1,2*M*(N-M)).
119 If LWORK = -1, then a workspace query is assumed; the routine
121 only calculates the optimal size of the WORK array, returns
122 this value as the first entry of the WORK array, and no error
123 message related to LWORK is issued by XERBLA.
124 IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
126 On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
127 LIWORK (input) INTEGER
129 The dimension of the array IWORK. If JOB = 'N' or 'E', LIWORK
130 >= 1; if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
131 If LIWORK = -1, then a workspace query is assumed; the routine
133 only calculates the optimal size of the IWORK array, returns
134 this value as the first entry of the IWORK array, and no error
135 message related to LIWORK is issued by XERBLA.
136 INFO (output) INTEGER
138 = 0: successful exit
139 < 0: if INFO = -i, the i-th argument had an illegal value
140 = 1: reordering of T failed because some eigenvalues are too
141 close to separate (the problem is very ill-conditioned); T may
142 have been partially reordered, and WR and WI contain the eigen‐
143 values in the same order as in T; S and SEP (if requested) are
144 set to zero.
145
147 STRSEN first collects the selected eigenvalues by computing an orthogo‐
148 nal transformation Z to move them to the top left corner of T. In
149 other words, the selected eigenvalues are the eigenvalues of T11 in:
150 Z'*T*Z = ( T11 T12 ) n1
152 ( 0 T22 ) n2
153 n1 n2
154 where N = n1+n2 and Z' means the transpose of Z. The first n1 columns
156 of Z span the specified invariant subspace of T.
157 If T has been obtained from the real Schur factorization of a matrix A
159 = Q*T*Q', then the reordered real Schur factorization of A is given by
160 A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
161 corresponding invariant subspace of A.
162 The reciprocal condition number of the average of the eigenvalues of
164 T11 may be returned in S. S lies between 0 (very badly conditioned) and
165 1 (very well conditioned). It is computed as follows. First we compute
166 R so that
167 P = ( I R ) n1
169 ( 0 0 ) n2
170 n1 n2
171 is the projector on the invariant subspace associated with T11. R is
173 the solution of the Sylvester equation:
174 T11*R - R*T22 = T12.
176 Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote the
178 two-norm of M. Then S is computed as the lower bound
179 (1 + F-norm(R)**2)**(-1/2)
181 on the reciprocal of 2-norm(P), the true reciprocal condition number.
183 S cannot underestimate 1 / 2-norm(P) by more than a factor of sqrt(N).
184 An approximate error bound for the computed average of the eigenvalues
186 of T11 is
187 EPS * norm(T) / S
189 where EPS is the machine precision.
191 The reciprocal condition number of the right invariant subspace spanned
193 by the first n1 columns of Z (or of Q*Z) is returned in SEP. SEP is
194 defined as the separation of T11 and T22:
195 sep( T11, T22 ) = sigma-min( C )
197 where sigma-min(C) is the smallest singular value of the
199 n1*n2-by-n1*n2 matrix
200 C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
202 I(m) is an m by m identity matrix, and kprod denotes the Kronecker
204 product. We estimate sigma-min(C) by the reciprocal of an estimate of
205 the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) can‐
206 not differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
207 When SEP is small, small changes in T can cause large changes in the
209 invariant subspace. An approximate bound on the maximum angular error
210 in the computed right invariant subspace is
211 EPS * norm(T) / SEP
213 LAPACK routine (version 3.1) November 2006 STRSEN(1)