Mercurial > octave-nkf

comparison libcruft/lapack/ctrsen.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE CTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
2 $ SEP, WORK, LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH.
9 *
10 * .. Scalar Arguments ..
11 CHARACTER COMPQ, JOB
12 INTEGER INFO, LDQ, LDT, LWORK, M, N
13 REAL S, SEP
14 * ..
15 * .. Array Arguments ..
16 LOGICAL SELECT( * )
17 COMPLEX Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * CTRSEN reorders the Schur factorization of a complex matrix
24 * A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
25 * the leading positions on the diagonal of the upper triangular matrix
26 * T, and the leading columns of Q form an orthonormal basis of the
27 * corresponding right invariant subspace.
28 *
29 * Optionally the routine computes the reciprocal condition numbers of
30 * the cluster of eigenvalues and/or the invariant subspace.
31 *
32 * Arguments
33 * =========
34 *
35 * JOB (input) CHARACTER*1
36 * Specifies whether condition numbers are required for the
37 * cluster of eigenvalues (S) or the invariant subspace (SEP):
38 * = 'N': none;
39 * = 'E': for eigenvalues only (S);
40 * = 'V': for invariant subspace only (SEP);
41 * = 'B': for both eigenvalues and invariant subspace (S and
42 * SEP).
43 *
44 * COMPQ (input) CHARACTER*1
45 * = 'V': update the matrix Q of Schur vectors;
46 * = 'N': do not update Q.
47 *
48 * SELECT (input) LOGICAL array, dimension (N)
49 * SELECT specifies the eigenvalues in the selected cluster. To
50 * select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
51 *
52 * N (input) INTEGER
53 * The order of the matrix T. N >= 0.
54 *
55 * T (input/output) COMPLEX array, dimension (LDT,N)
56 * On entry, the upper triangular matrix T.
57 * On exit, T is overwritten by the reordered matrix T, with the
58 * selected eigenvalues as the leading diagonal elements.
59 *
60 * LDT (input) INTEGER
61 * The leading dimension of the array T. LDT >= max(1,N).
62 *
63 * Q (input/output) COMPLEX array, dimension (LDQ,N)
64 * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
65 * On exit, if COMPQ = 'V', Q has been postmultiplied by the
66 * unitary transformation matrix which reorders T; the leading M
67 * columns of Q form an orthonormal basis for the specified
68 * invariant subspace.
69 * If COMPQ = 'N', Q is not referenced.
70 *
71 * LDQ (input) INTEGER
72 * The leading dimension of the array Q.
73 * LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
74 *
75 * W (output) COMPLEX array, dimension (N)
76 * The reordered eigenvalues of T, in the same order as they
77 * appear on the diagonal of T.
78 *
79 * M (output) INTEGER
80 * The dimension of the specified invariant subspace.
81 * 0 <= M <= N.
82 *
83 * S (output) REAL
84 * If JOB = 'E' or 'B', S is a lower bound on the reciprocal
85 * condition number for the selected cluster of eigenvalues.
86 * S cannot underestimate the true reciprocal condition number
87 * by more than a factor of sqrt(N). If M = 0 or N, S = 1.
88 * If JOB = 'N' or 'V', S is not referenced.
89 *
90 * SEP (output) REAL
91 * If JOB = 'V' or 'B', SEP is the estimated reciprocal
92 * condition number of the specified invariant subspace. If
93 * M = 0 or N, SEP = norm(T).
94 * If JOB = 'N' or 'E', SEP is not referenced.
95 *
96 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
97 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
98 *
99 * LWORK (input) INTEGER
100 * The dimension of the array WORK.
101 * If JOB = 'N', LWORK >= 1;
102 * if JOB = 'E', LWORK = max(1,M*(N-M));
103 * if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
104 *
105 * If LWORK = -1, then a workspace query is assumed; the routine
106 * only calculates the optimal size of the WORK array, returns
107 * this value as the first entry of the WORK array, and no error
108 * message related to LWORK is issued by XERBLA.
109 *
110 * INFO (output) INTEGER
111 * = 0: successful exit
112 * < 0: if INFO = -i, the i-th argument had an illegal value
113 *
114 * Further Details
115 * ===============
116 *
117 * CTRSEN first collects the selected eigenvalues by computing a unitary
118 * transformation Z to move them to the top left corner of T. In other
119 * words, the selected eigenvalues are the eigenvalues of T11 in:
120 *
121 * Z'*T*Z = ( T11 T12 ) n1
122 * ( 0 T22 ) n2
123 * n1 n2
124 *
125 * where N = n1+n2 and Z' means the conjugate transpose of Z. The first
126 * n1 columns of Z span the specified invariant subspace of T.
127 *
128 * If T has been obtained from the Schur factorization of a matrix
129 * A = Q*T*Q', then the reordered Schur factorization of A is given by
130 * A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
131 * corresponding invariant subspace of A.
132 *
133 * The reciprocal condition number of the average of the eigenvalues of
134 * T11 may be returned in S. S lies between 0 (very badly conditioned)
135 * and 1 (very well conditioned). It is computed as follows. First we
136 * compute R so that
137 *
138 * P = ( I R ) n1
139 * ( 0 0 ) n2
140 * n1 n2
141 *
142 * is the projector on the invariant subspace associated with T11.
143 * R is the solution of the Sylvester equation:
144 *
145 * T11*R - R*T22 = T12.
146 *
147 * Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
148 * the two-norm of M. Then S is computed as the lower bound
149 *
150 * (1 + F-norm(R)**2)**(-1/2)
151 *
152 * on the reciprocal of 2-norm(P), the true reciprocal condition number.
153 * S cannot underestimate 1 / 2-norm(P) by more than a factor of
154 * sqrt(N).
155 *
156 * An approximate error bound for the computed average of the
157 * eigenvalues of T11 is
158 *
159 * EPS * norm(T) / S
160 *
161 * where EPS is the machine precision.
162 *
163 * The reciprocal condition number of the right invariant subspace
164 * spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
165 * SEP is defined as the separation of T11 and T22:
166 *
167 * sep( T11, T22 ) = sigma-min( C )
168 *
169 * where sigma-min(C) is the smallest singular value of the
170 * n1*n2-by-n1*n2 matrix
171 *
172 * C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
173 *
174 * I(m) is an m by m identity matrix, and kprod denotes the Kronecker
175 * product. We estimate sigma-min(C) by the reciprocal of an estimate of
176 * the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
177 * cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
178 *
179 * When SEP is small, small changes in T can cause large changes in
180 * the invariant subspace. An approximate bound on the maximum angular
181 * error in the computed right invariant subspace is
182 *
183 * EPS * norm(T) / SEP
184 *
185 * =====================================================================
186 *
187 * .. Parameters ..
188 REAL ZERO, ONE
189 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
190 * ..
191 * .. Local Scalars ..
192 LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
193 INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
194 REAL EST, RNORM, SCALE
195 * ..
196 * .. Local Arrays ..
197 INTEGER ISAVE( 3 )
198 REAL RWORK( 1 )
199 * ..
200 * .. External Functions ..
201 LOGICAL LSAME
202 REAL CLANGE
203 EXTERNAL LSAME, CLANGE
204 * ..
205 * .. External Subroutines ..
206 EXTERNAL CLACN2, CLACPY, CTREXC, CTRSYL, XERBLA
207 * ..
208 * .. Intrinsic Functions ..
209 INTRINSIC MAX, SQRT
210 * ..
211 * .. Executable Statements ..
212 *
213 * Decode and test the input parameters.
214 *
215 WANTBH = LSAME( JOB, 'B' )
216 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
217 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
218 WANTQ = LSAME( COMPQ, 'V' )
219 *
220 * Set M to the number of selected eigenvalues.
221 *
222 M = 0
223 DO 10 K = 1, N
224 IF( SELECT( K ) )
225 $ M = M + 1
226 10 CONTINUE
227 *
228 N1 = M
229 N2 = N - M
230 NN = N1*N2
231 *
232 INFO = 0
233 LQUERY = ( LWORK.EQ.-1 )
234 *
235 IF( WANTSP ) THEN
236 LWMIN = MAX( 1, 2*NN )
237 ELSE IF( LSAME( JOB, 'N' ) ) THEN
238 LWMIN = 1
239 ELSE IF( LSAME( JOB, 'E' ) ) THEN
240 LWMIN = MAX( 1, NN )
241 END IF
242 *
243 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
244 $ THEN
245 INFO = -1
246 ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
247 INFO = -2
248 ELSE IF( N.LT.0 ) THEN
249 INFO = -4
250 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
251 INFO = -6
252 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
253 INFO = -8
254 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
255 INFO = -14
256 END IF
257 *
258 IF( INFO.EQ.0 ) THEN
259 WORK( 1 ) = LWMIN
260 END IF
261 *
262 IF( INFO.NE.0 ) THEN
263 CALL XERBLA( 'CTRSEN', -INFO )
264 RETURN
265 ELSE IF( LQUERY ) THEN
266 RETURN
267 END IF
268 *
269 * Quick return if possible
270 *
271 IF( M.EQ.N .OR. M.EQ.0 ) THEN
272 IF( WANTS )
273 $ S = ONE
274 IF( WANTSP )
275 $ SEP = CLANGE( '1', N, N, T, LDT, RWORK )
276 GO TO 40
277 END IF
278 *
279 * Collect the selected eigenvalues at the top left corner of T.
280 *
281 KS = 0
282 DO 20 K = 1, N
283 IF( SELECT( K ) ) THEN
284 KS = KS + 1
285 *
286 * Swap the K-th eigenvalue to position KS.
287 *
288 IF( K.NE.KS )
289 $ CALL CTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
290 END IF
291 20 CONTINUE
292 *
293 IF( WANTS ) THEN
294 *
295 * Solve the Sylvester equation for R:
296 *
297 * T11*R - R*T22 = scale*T12
298 *
299 CALL CLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
300 CALL CTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
301 $ LDT, WORK, N1, SCALE, IERR )
302 *
303 * Estimate the reciprocal of the condition number of the cluster
304 * of eigenvalues.
305 *
306 RNORM = CLANGE( 'F', N1, N2, WORK, N1, RWORK )
307 IF( RNORM.EQ.ZERO ) THEN
308 S = ONE
309 ELSE
310 S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
311 $ SQRT( RNORM ) )
312 END IF
313 END IF
314 *
315 IF( WANTSP ) THEN
316 *
317 * Estimate sep(T11,T22).
318 *
319 EST = ZERO
320 KASE = 0
321 30 CONTINUE
322 CALL CLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE )
323 IF( KASE.NE.0 ) THEN
324 IF( KASE.EQ.1 ) THEN
325 *
326 * Solve T11*R - R*T22 = scale*X.
327 *
328 CALL CTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
329 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
330 $ IERR )
331 ELSE
332 *
333 * Solve T11'*R - R*T22' = scale*X.
334 *
335 CALL CTRSYL( 'C', 'C', -1, N1, N2, T, LDT,
336 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
337 $ IERR )
338 END IF
339 GO TO 30
340 END IF
341 *
342 SEP = SCALE / EST
343 END IF
344 *
345 40 CONTINUE
346 *
347 * Copy reordered eigenvalues to W.
348 *
349 DO 50 K = 1, N
350 W( K ) = T( K, K )
351 50 CONTINUE
352 *
353 WORK( 1 ) = LWMIN
354 *
355 RETURN
356 *
357 * End of CTRSEN
358 *
359 END