1DTRSEN(1) LAPACK routine (version 3.1) DTRSEN(1)
2
6 DTRSEN - 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 DTRSEN( 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 DOUBLE PRECISION S, SEP
19 LOGICAL SELECT( * )
21 INTEGER IWORK( * )
23 DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( *
25 ), WR( * )
26
28 DTRSEN 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 DHSEQR), 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N)
87 WI (output) DOUBLE PRECISION array, dimension (N) The real
88 and imaginary parts, respectively, of the reordered eigenvalues
89 of T. The eigenvalues are stored in the same order as on the
90 diagonal of T, with WR(i) = T(i,i) and, if T(i:i+1,i:i+1) is a
91 2-by-2 diagonal block, WI(i) > 0 and WI(i+1) = -WI(i). Note
92 that if a complex eigenvalue is sufficiently ill-conditioned,
93 then its value may differ significantly from its value before
94 reordering.
95 M (output) INTEGER
97 The dimension of the specified invariant subspace. 0 < = M <=
98 N.
99 S (output) DOUBLE PRECISION
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) DOUBLE PRECISION
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) DOUBLE PRECISION array, dimension
113 (MAX(1,LWORK))
114 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
115 LWORK (input) INTEGER
117 The dimension of the array WORK. If JOB = 'N', LWORK >=
118 max(1,N); if JOB = 'E', LWORK >= max(1,M*(N-M)); if JOB = 'V'
119 or 'B', LWORK >= max(1,2*M*(N-M)).
120 If LWORK = -1, then a workspace query is assumed; the routine
122 only calculates the optimal size of the WORK array, returns
123 this value as the first entry of the WORK array, and no error
124 message related to LWORK is issued by XERBLA.
125 IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
127 On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
128 LIWORK (input) INTEGER
130 The dimension of the array IWORK. If JOB = 'N' or 'E', LIWORK
131 >= 1; if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
132 If LIWORK = -1, then a workspace query is assumed; the routine
134 only calculates the optimal size of the IWORK array, returns
135 this value as the first entry of the IWORK array, and no error
136 message related to LIWORK is issued by XERBLA.
137 INFO (output) INTEGER
139 = 0: successful exit
140 < 0: if INFO = -i, the i-th argument had an illegal value
141 = 1: reordering of T failed because some eigenvalues are too
142 close to separate (the problem is very ill-conditioned); T may
143 have been partially reordered, and WR and WI contain the eigen‐
144 values in the same order as in T; S and SEP (if requested) are
145 set to zero.
146
148 DTRSEN first collects the selected eigenvalues by computing an orthogo‐
149 nal transformation Z to move them to the top left corner of T. In
150 other words, the selected eigenvalues are the eigenvalues of T11 in:
151 Z'*T*Z = ( T11 T12 ) n1
153 ( 0 T22 ) n2
154 n1 n2
155 where N = n1+n2 and Z' means the transpose of Z. The first n1 columns
157 of Z span the specified invariant subspace of T.
158 If T has been obtained from the real Schur factorization of a matrix A
160 = Q*T*Q', then the reordered real Schur factorization of A is given by
161 A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
162 corresponding invariant subspace of A.
163 The reciprocal condition number of the average of the eigenvalues of
165 T11 may be returned in S. S lies between 0 (very badly conditioned) and
166 1 (very well conditioned). It is computed as follows. First we compute
167 R so that
168 P = ( I R ) n1
170 ( 0 0 ) n2
171 n1 n2
172 is the projector on the invariant subspace associated with T11. R is
174 the solution of the Sylvester equation:
175 T11*R - R*T22 = T12.
177 Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote the
179 two-norm of M. Then S is computed as the lower bound
180 (1 + F-norm(R)**2)**(-1/2)
182 on the reciprocal of 2-norm(P), the true reciprocal condition number.
184 S cannot underestimate 1 / 2-norm(P) by more than a factor of sqrt(N).
185 An approximate error bound for the computed average of the eigenvalues
187 of T11 is
188 EPS * norm(T) / S
190 where EPS is the machine precision.
192 The reciprocal condition number of the right invariant subspace spanned
194 by the first n1 columns of Z (or of Q*Z) is returned in SEP. SEP is
195 defined as the separation of T11 and T22:
196 sep( T11, T22 ) = sigma-min( C )
198 where sigma-min(C) is the smallest singular value of the
200 n1*n2-by-n1*n2 matrix
201 C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
203 I(m) is an m by m identity matrix, and kprod denotes the Kronecker
205 product. We estimate sigma-min(C) by the reciprocal of an estimate of
206 the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) can‐
207 not differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
208 When SEP is small, small changes in T can cause large changes in the
210 invariant subspace. An approximate bound on the maximum angular error
211 in the computed right invariant subspace is
212 EPS * norm(T) / SEP
214 LAPACK routine (version 3.1) November 2006 DTRSEN(1)