1: *> \brief \b ZTRSEN
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTRSEN + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsen.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsen.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsen.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
22: * SEP, WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER COMPQ, JOB
26: * INTEGER INFO, LDQ, LDT, LWORK, M, N
27: * DOUBLE PRECISION S, SEP
28: * ..
29: * .. Array Arguments ..
30: * LOGICAL SELECT( * )
31: * COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> ZTRSEN reorders the Schur factorization of a complex matrix
41: *> A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
42: *> the leading positions on the diagonal of the upper triangular matrix
43: *> T, and the leading columns of Q form an orthonormal basis of the
44: *> corresponding right invariant subspace.
45: *>
46: *> Optionally the routine computes the reciprocal condition numbers of
47: *> the cluster of eigenvalues and/or the invariant subspace.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] JOB
54: *> \verbatim
55: *> JOB is CHARACTER*1
56: *> Specifies whether condition numbers are required for the
57: *> cluster of eigenvalues (S) or the invariant subspace (SEP):
58: *> = 'N': none;
59: *> = 'E': for eigenvalues only (S);
60: *> = 'V': for invariant subspace only (SEP);
61: *> = 'B': for both eigenvalues and invariant subspace (S and
62: *> SEP).
63: *> \endverbatim
64: *>
65: *> \param[in] COMPQ
66: *> \verbatim
67: *> COMPQ is CHARACTER*1
68: *> = 'V': update the matrix Q of Schur vectors;
69: *> = 'N': do not update Q.
70: *> \endverbatim
71: *>
72: *> \param[in] SELECT
73: *> \verbatim
74: *> SELECT is LOGICAL array, dimension (N)
75: *> SELECT specifies the eigenvalues in the selected cluster. To
76: *> select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
77: *> \endverbatim
78: *>
79: *> \param[in] N
80: *> \verbatim
81: *> N is INTEGER
82: *> The order of the matrix T. N >= 0.
83: *> \endverbatim
84: *>
85: *> \param[in,out] T
86: *> \verbatim
87: *> T is COMPLEX*16 array, dimension (LDT,N)
88: *> On entry, the upper triangular matrix T.
89: *> On exit, T is overwritten by the reordered matrix T, with the
90: *> selected eigenvalues as the leading diagonal elements.
91: *> \endverbatim
92: *>
93: *> \param[in] LDT
94: *> \verbatim
95: *> LDT is INTEGER
96: *> The leading dimension of the array T. LDT >= max(1,N).
97: *> \endverbatim
98: *>
99: *> \param[in,out] Q
100: *> \verbatim
101: *> Q is COMPLEX*16 array, dimension (LDQ,N)
102: *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
103: *> On exit, if COMPQ = 'V', Q has been postmultiplied by the
104: *> unitary transformation matrix which reorders T; the leading M
105: *> columns of Q form an orthonormal basis for the specified
106: *> invariant subspace.
107: *> If COMPQ = 'N', Q is not referenced.
108: *> \endverbatim
109: *>
110: *> \param[in] LDQ
111: *> \verbatim
112: *> LDQ is INTEGER
113: *> The leading dimension of the array Q.
114: *> LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
115: *> \endverbatim
116: *>
117: *> \param[out] W
118: *> \verbatim
119: *> W is COMPLEX*16 array, dimension (N)
120: *> The reordered eigenvalues of T, in the same order as they
121: *> appear on the diagonal of T.
122: *> \endverbatim
123: *>
124: *> \param[out] M
125: *> \verbatim
126: *> M is INTEGER
127: *> The dimension of the specified invariant subspace.
128: *> 0 <= M <= N.
129: *> \endverbatim
130: *>
131: *> \param[out] S
132: *> \verbatim
133: *> S is DOUBLE PRECISION
134: *> If JOB = 'E' or 'B', S is a lower bound on the reciprocal
135: *> condition number for the selected cluster of eigenvalues.
136: *> S cannot underestimate the true reciprocal condition number
137: *> by more than a factor of sqrt(N). If M = 0 or N, S = 1.
138: *> If JOB = 'N' or 'V', S is not referenced.
139: *> \endverbatim
140: *>
141: *> \param[out] SEP
142: *> \verbatim
143: *> SEP is DOUBLE PRECISION
144: *> If JOB = 'V' or 'B', SEP is the estimated reciprocal
145: *> condition number of the specified invariant subspace. If
146: *> M = 0 or N, SEP = norm(T).
147: *> If JOB = 'N' or 'E', SEP is not referenced.
148: *> \endverbatim
149: *>
150: *> \param[out] WORK
151: *> \verbatim
152: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
153: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
154: *> \endverbatim
155: *>
156: *> \param[in] LWORK
157: *> \verbatim
158: *> LWORK is INTEGER
159: *> The dimension of the array WORK.
160: *> If JOB = 'N', LWORK >= 1;
161: *> if JOB = 'E', LWORK = max(1,M*(N-M));
162: *> if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
163: *>
164: *> If LWORK = -1, then a workspace query is assumed; the routine
165: *> only calculates the optimal size of the WORK array, returns
166: *> this value as the first entry of the WORK array, and no error
167: *> message related to LWORK is issued by XERBLA.
168: *> \endverbatim
169: *>
170: *> \param[out] INFO
171: *> \verbatim
172: *> INFO is INTEGER
173: *> = 0: successful exit
174: *> < 0: if INFO = -i, the i-th argument had an illegal value
175: *> \endverbatim
176: *
177: * Authors:
178: * ========
179: *
180: *> \author Univ. of Tennessee
181: *> \author Univ. of California Berkeley
182: *> \author Univ. of Colorado Denver
183: *> \author NAG Ltd.
184: *
185: *> \ingroup complex16OTHERcomputational
186: *
187: *> \par Further Details:
188: * =====================
189: *>
190: *> \verbatim
191: *>
192: *> ZTRSEN first collects the selected eigenvalues by computing a unitary
193: *> transformation Z to move them to the top left corner of T. In other
194: *> words, the selected eigenvalues are the eigenvalues of T11 in:
195: *>
196: *> Z**H * T * Z = ( T11 T12 ) n1
197: *> ( 0 T22 ) n2
198: *> n1 n2
199: *>
200: *> where N = n1+n2. The first
201: *> n1 columns of Z span the specified invariant subspace of T.
202: *>
203: *> If T has been obtained from the Schur factorization of a matrix
204: *> A = Q*T*Q**H, then the reordered Schur factorization of A is given by
205: *> A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
206: *> corresponding invariant subspace of A.
207: *>
208: *> The reciprocal condition number of the average of the eigenvalues of
209: *> T11 may be returned in S. S lies between 0 (very badly conditioned)
210: *> and 1 (very well conditioned). It is computed as follows. First we
211: *> compute R so that
212: *>
213: *> P = ( I R ) n1
214: *> ( 0 0 ) n2
215: *> n1 n2
216: *>
217: *> is the projector on the invariant subspace associated with T11.
218: *> R is the solution of the Sylvester equation:
219: *>
220: *> T11*R - R*T22 = T12.
221: *>
222: *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
223: *> the two-norm of M. Then S is computed as the lower bound
224: *>
225: *> (1 + F-norm(R)**2)**(-1/2)
226: *>
227: *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
228: *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
229: *> sqrt(N).
230: *>
231: *> An approximate error bound for the computed average of the
232: *> eigenvalues of T11 is
233: *>
234: *> EPS * norm(T) / S
235: *>
236: *> where EPS is the machine precision.
237: *>
238: *> The reciprocal condition number of the right invariant subspace
239: *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
240: *> SEP is defined as the separation of T11 and T22:
241: *>
242: *> sep( T11, T22 ) = sigma-min( C )
243: *>
244: *> where sigma-min(C) is the smallest singular value of the
245: *> n1*n2-by-n1*n2 matrix
246: *>
247: *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
248: *>
249: *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
250: *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
251: *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
252: *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
253: *>
254: *> When SEP is small, small changes in T can cause large changes in
255: *> the invariant subspace. An approximate bound on the maximum angular
256: *> error in the computed right invariant subspace is
257: *>
258: *> EPS * norm(T) / SEP
259: *> \endverbatim
260: *>
261: * =====================================================================
262: SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
263: $ SEP, WORK, LWORK, INFO )
264: *
265: * -- LAPACK computational routine --
266: * -- LAPACK is a software package provided by Univ. of Tennessee, --
267: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268: *
269: * .. Scalar Arguments ..
270: CHARACTER COMPQ, JOB
271: INTEGER INFO, LDQ, LDT, LWORK, M, N
272: DOUBLE PRECISION S, SEP
273: * ..
274: * .. Array Arguments ..
275: LOGICAL SELECT( * )
276: COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
277: * ..
278: *
279: * =====================================================================
280: *
281: * .. Parameters ..
282: DOUBLE PRECISION ZERO, ONE
283: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
284: * ..
285: * .. Local Scalars ..
286: LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
287: INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
288: DOUBLE PRECISION EST, RNORM, SCALE
289: * ..
290: * .. Local Arrays ..
291: INTEGER ISAVE( 3 )
292: DOUBLE PRECISION RWORK( 1 )
293: * ..
294: * .. External Functions ..
295: LOGICAL LSAME
296: DOUBLE PRECISION ZLANGE
297: EXTERNAL LSAME, ZLANGE
298: * ..
299: * .. External Subroutines ..
300: EXTERNAL XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL
301: * ..
302: * .. Intrinsic Functions ..
303: INTRINSIC MAX, SQRT
304: * ..
305: * .. Executable Statements ..
306: *
307: * Decode and test the input parameters.
308: *
309: WANTBH = LSAME( JOB, 'B' )
310: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
311: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
312: WANTQ = LSAME( COMPQ, 'V' )
313: *
314: * Set M to the number of selected eigenvalues.
315: *
316: M = 0
317: DO 10 K = 1, N
318: IF( SELECT( K ) )
319: $ M = M + 1
320: 10 CONTINUE
321: *
322: N1 = M
323: N2 = N - M
324: NN = N1*N2
325: *
326: INFO = 0
327: LQUERY = ( LWORK.EQ.-1 )
328: *
329: IF( WANTSP ) THEN
330: LWMIN = MAX( 1, 2*NN )
331: ELSE IF( LSAME( JOB, 'N' ) ) THEN
332: LWMIN = 1
333: ELSE IF( LSAME( JOB, 'E' ) ) THEN
334: LWMIN = MAX( 1, NN )
335: END IF
336: *
337: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
338: $ THEN
339: INFO = -1
340: ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
341: INFO = -2
342: ELSE IF( N.LT.0 ) THEN
343: INFO = -4
344: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
345: INFO = -6
346: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
347: INFO = -8
348: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
349: INFO = -14
350: END IF
351: *
352: IF( INFO.EQ.0 ) THEN
353: WORK( 1 ) = LWMIN
354: END IF
355: *
356: IF( INFO.NE.0 ) THEN
357: CALL XERBLA( 'ZTRSEN', -INFO )
358: RETURN
359: ELSE IF( LQUERY ) THEN
360: RETURN
361: END IF
362: *
363: * Quick return if possible
364: *
365: IF( M.EQ.N .OR. M.EQ.0 ) THEN
366: IF( WANTS )
367: $ S = ONE
368: IF( WANTSP )
369: $ SEP = ZLANGE( '1', N, N, T, LDT, RWORK )
370: GO TO 40
371: END IF
372: *
373: * Collect the selected eigenvalues at the top left corner of T.
374: *
375: KS = 0
376: DO 20 K = 1, N
377: IF( SELECT( K ) ) THEN
378: KS = KS + 1
379: *
380: * Swap the K-th eigenvalue to position KS.
381: *
382: IF( K.NE.KS )
383: $ CALL ZTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
384: END IF
385: 20 CONTINUE
386: *
387: IF( WANTS ) THEN
388: *
389: * Solve the Sylvester equation for R:
390: *
391: * T11*R - R*T22 = scale*T12
392: *
393: CALL ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
394: CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
395: $ LDT, WORK, N1, SCALE, IERR )
396: *
397: * Estimate the reciprocal of the condition number of the cluster
398: * of eigenvalues.
399: *
400: RNORM = ZLANGE( 'F', N1, N2, WORK, N1, RWORK )
401: IF( RNORM.EQ.ZERO ) THEN
402: S = ONE
403: ELSE
404: S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
405: $ SQRT( RNORM ) )
406: END IF
407: END IF
408: *
409: IF( WANTSP ) THEN
410: *
411: * Estimate sep(T11,T22).
412: *
413: EST = ZERO
414: KASE = 0
415: 30 CONTINUE
416: CALL ZLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE )
417: IF( KASE.NE.0 ) THEN
418: IF( KASE.EQ.1 ) THEN
419: *
420: * Solve T11*R - R*T22 = scale*X.
421: *
422: CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
423: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
424: $ IERR )
425: ELSE
426: *
427: * Solve T11**H*R - R*T22**H = scale*X.
428: *
429: CALL ZTRSYL( 'C', 'C', -1, N1, N2, T, LDT,
430: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
431: $ IERR )
432: END IF
433: GO TO 30
434: END IF
435: *
436: SEP = SCALE / EST
437: END IF
438: *
439: 40 CONTINUE
440: *
441: * Copy reordered eigenvalues to W.
442: *
443: DO 50 K = 1, N
444: W( K ) = T( K, K )
445: 50 CONTINUE
446: *
447: WORK( 1 ) = LWMIN
448: *
449: RETURN
450: *
451: * End of ZTRSEN
452: *
453: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztrsen.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack