1STRSEN(1) LAPACK routine (version 3.2) STRSEN(1)
2
6 STRSEN - reorders the real Schur factorization of a real matrix A =
7 Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
8 leading diagonal 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
34 cluster of eigenvalues and/or the invariant subspace. T must be in
35 Schur canonical form (as returned by SHSEQR), that is, block upper tri‐
36 angular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal
37 block has its diagonal elemnts equal and its off-diagonal elements of
38 opposite sign.
39
41 JOB (input) CHARACTER*1
42 Specifies whether condition numbers are required for the clus‐
43 ter of eigenvalues (S) or the invariant subspace (SEP):
44 = 'N': none;
45 = 'E': for eigenvalues only (S);
46 = 'V': for invariant subspace only (SEP);
47 = 'B': for both eigenvalues and invariant subspace (S and SEP).
48 COMPQ (input) CHARACTER*1
50 = 'V': update the matrix Q of Schur vectors;
51 = 'N': do not update Q.
52 SELECT (input) LOGICAL array, dimension (N)
54 SELECT specifies the eigenvalues in the selected cluster. To
55 select a real eigenvalue w(j), SELECT(j) must be set to .TRUE..
56 To select a complex conjugate pair of eigenvalues w(j) and
57 w(j+1), corresponding to a 2-by-2 diagonal block, either
58 SELECT(j) or SELECT(j+1) or both must be set to .TRUE.; a com‐
59 plex conjugate pair of eigenvalues must be either both included
60 in the cluster or both excluded.
61 N (input) INTEGER
63 The order of the matrix T. N >= 0.
64 T (input/output) REAL array, dimension (LDT,N)
66 On entry, the upper quasi-triangular matrix T, in Schur canoni‐
67 cal form. On exit, T is overwritten by the reordered matrix T,
68 again in Schur canonical form, with the selected eigenvalues in
69 the leading diagonal blocks.
70 LDT (input) INTEGER
72 The leading dimension of the array T. LDT >= max(1,N).
73 Q (input/output) REAL array, dimension (LDQ,N)
75 On entry, if COMPQ = 'V', the matrix Q of Schur vectors. On
76 exit, if COMPQ = 'V', Q has been postmultiplied by the orthogo‐
77 nal transformation matrix which reorders T; the leading M col‐
78 umns of Q form an orthonormal basis for the specified invariant
79 subspace. If COMPQ = 'N', Q is not referenced.
80 LDQ (input) INTEGER
82 The leading dimension of the array Q. LDQ >= 1; and if COMPQ =
83 'V', LDQ >= N.
84 WR (output) REAL array, dimension (N)
86 WI (output) REAL array, dimension (N) The real and imagi‐
87 nary parts, respectively, of the reordered eigenvalues of T.
88 The eigenvalues are stored in the same order as on the diagonal
89 of T, with WR(i) = T(i,i) and, if T(i:i+1,i:i+1) is a 2-by-2
90 diagonal block, WI(i) > 0 and WI(i+1) = -WI(i). Note that if a
91 complex eigenvalue is sufficiently ill-conditioned, then its
92 value may differ significantly from its value before reorder‐
93 ing.
94 M (output) INTEGER
96 The dimension of the specified invariant subspace. 0 < = M <=
97 N.
98 S (output) REAL
100 If JOB = 'E' or 'B', S is a lower bound on the reciprocal con‐
101 dition number for the selected cluster of eigenvalues. S can‐
102 not underestimate the true reciprocal condition number by more
103 than a factor of sqrt(N). If M = 0 or N, S = 1. If JOB = 'N'
104 or 'V', S is not referenced.
105 SEP (output) REAL
107 If JOB = 'V' or 'B', SEP is the estimated reciprocal condition
108 number of the specified invariant subspace. If M = 0 or N, SEP
109 = norm(T). If JOB = 'N' or 'E', SEP is not referenced.
110 WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
112 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
113 LWORK (input) INTEGER
115 The dimension of the array WORK. If JOB = 'N', LWORK >=
116 max(1,N); if JOB = 'E', LWORK >= max(1,M*(N-M)); if JOB = 'V'
117 or 'B', LWORK >= max(1,2*M*(N-M)). If LWORK = -1, then a
118 workspace query is assumed; the routine only calculates the
119 optimal size of the WORK array, returns this value as the first
120 entry of the WORK array, and no error message related to LWORK
121 is issued by XERBLA.
122 IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
124 On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
125 LIWORK (input) INTEGER
127 The dimension of the array IWORK. If JOB = 'N' or 'E', LIWORK
128 >= 1; if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)). If LIWORK
129 = -1, then a workspace query is assumed; the routine only cal‐
130 culates the optimal size of the IWORK array, returns this value
131 as the first entry of the IWORK array, and no error message
132 related to LIWORK is issued by XERBLA.
133 INFO (output) INTEGER
135 = 0: successful exit
136 < 0: if INFO = -i, the i-th argument had an illegal value
137 = 1: reordering of T failed because some eigenvalues are too
138 close to separate (the problem is very ill-conditioned); T may
139 have been partially reordered, and WR and WI contain the eigen‐
140 values in the same order as in T; S and SEP (if requested) are
141 set to zero.
142
144 STRSEN first collects the selected eigenvalues by computing an orthogo‐
145 nal transformation Z to move them to the top left corner of T. In
146 other words, the selected eigenvalues are the eigenvalues of T11 in:
147 Z'*T*Z = ( T11 T12 ) n1
148 ( 0 T22 ) n2
149 n1 n2
150 where N = n1+n2 and Z' means the transpose of Z. The first n1 columns
151 of Z span the specified invariant subspace of T.
152 If T has been obtained from the real Schur factorization of a matrix A
153 = Q*T*Q', then the reordered real Schur factorization of A is given by
154 A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
155 corresponding invariant subspace of A.
156 The reciprocal condition number of the average of the eigenvalues of
157 T11 may be returned in S. S lies between 0 (very badly conditioned) and
158 1 (very well conditioned). It is computed as follows. First we compute
159 R so that
160 P = ( I R ) n1
161 ( 0 0 ) n2
162 n1 n2
163 is the projector on the invariant subspace associated with T11. R is
164 the solution of the Sylvester equation:
165 T11*R - R*T22 = T12.
166 Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote the
167 two-norm of M. Then S is computed as the lower bound
168 (1 + F-norm(R)**2)**(-1/2)
169 on the reciprocal of 2-norm(P), the true reciprocal condition number.
170 S cannot underestimate 1 / 2-norm(P) by more than a factor of sqrt(N).
171 An approximate error bound for the computed average of the eigenvalues
172 of T11 is
173 EPS * norm(T) / S
174 where EPS is the machine precision.
175 The reciprocal condition number of the right invariant subspace spanned
176 by the first n1 columns of Z (or of Q*Z) is returned in SEP. SEP is
177 defined as the separation of T11 and T22:
178 sep( T11, T22 ) = sigma-min( C )
179 where sigma-min(C) is the smallest singular value of the
180 n1*n2-by-n1*n2 matrix
181 C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) I(m) is an
182 m by m identity matrix, and kprod denotes the Kronecker product. We
183 estimate sigma-min(C) by the reciprocal of an estimate of the 1-norm of
184 inverse(C). The true reciprocal 1-norm of inverse(C) cannot differ from
185 sigma-min(C) by more than a factor of sqrt(n1*n2). When SEP is small,
186 small changes in T can cause large changes in the invariant subspace.
187 An approximate bound on the maximum angular error in the computed right
188 invariant subspace is
189 EPS * norm(T) / SEP
190 LAPACK routine (version 3.2) November 2008 STRSEN(1)