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