Annotation of rpl/lapack/lapack/ztrsen.f, revision 1.14
1.9 bertrand 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: *> \date November 2011
186: *
187: *> \ingroup complex16OTHERcomputational
188: *
189: *> \par Further Details:
190: * =====================
191: *>
192: *> \verbatim
193: *>
194: *> ZTRSEN first collects the selected eigenvalues by computing a unitary
195: *> transformation Z to move them to the top left corner of T. In other
196: *> words, the selected eigenvalues are the eigenvalues of T11 in:
197: *>
198: *> Z**H * T * Z = ( T11 T12 ) n1
199: *> ( 0 T22 ) n2
200: *> n1 n2
201: *>
202: *> where N = n1+n2. The first
203: *> n1 columns of Z span the specified invariant subspace of T.
204: *>
205: *> If T has been obtained from the Schur factorization of a matrix
206: *> A = Q*T*Q**H, then the reordered Schur factorization of A is given by
207: *> A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
208: *> corresponding invariant subspace of A.
209: *>
210: *> The reciprocal condition number of the average of the eigenvalues of
211: *> T11 may be returned in S. S lies between 0 (very badly conditioned)
212: *> and 1 (very well conditioned). It is computed as follows. First we
213: *> compute R so that
214: *>
215: *> P = ( I R ) n1
216: *> ( 0 0 ) n2
217: *> n1 n2
218: *>
219: *> is the projector on the invariant subspace associated with T11.
220: *> R is the solution of the Sylvester equation:
221: *>
222: *> T11*R - R*T22 = T12.
223: *>
224: *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
225: *> the two-norm of M. Then S is computed as the lower bound
226: *>
227: *> (1 + F-norm(R)**2)**(-1/2)
228: *>
229: *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
230: *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
231: *> sqrt(N).
232: *>
233: *> An approximate error bound for the computed average of the
234: *> eigenvalues of T11 is
235: *>
236: *> EPS * norm(T) / S
237: *>
238: *> where EPS is the machine precision.
239: *>
240: *> The reciprocal condition number of the right invariant subspace
241: *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
242: *> SEP is defined as the separation of T11 and T22:
243: *>
244: *> sep( T11, T22 ) = sigma-min( C )
245: *>
246: *> where sigma-min(C) is the smallest singular value of the
247: *> n1*n2-by-n1*n2 matrix
248: *>
249: *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
250: *>
251: *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
252: *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
253: *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
254: *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
255: *>
256: *> When SEP is small, small changes in T can cause large changes in
257: *> the invariant subspace. An approximate bound on the maximum angular
258: *> error in the computed right invariant subspace is
259: *>
260: *> EPS * norm(T) / SEP
261: *> \endverbatim
262: *>
263: * =====================================================================
1.1 bertrand 264: SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
265: $ SEP, WORK, LWORK, INFO )
266: *
1.9 bertrand 267: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 268: * -- LAPACK is a software package provided by Univ. of Tennessee, --
269: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 bertrand 270: * November 2011
1.1 bertrand 271: *
272: * .. Scalar Arguments ..
273: CHARACTER COMPQ, JOB
274: INTEGER INFO, LDQ, LDT, LWORK, M, N
275: DOUBLE PRECISION S, SEP
276: * ..
277: * .. Array Arguments ..
278: LOGICAL SELECT( * )
279: COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
280: * ..
281: *
282: * =====================================================================
283: *
284: * .. Parameters ..
285: DOUBLE PRECISION ZERO, ONE
286: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
287: * ..
288: * .. Local Scalars ..
289: LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
290: INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
291: DOUBLE PRECISION EST, RNORM, SCALE
292: * ..
293: * .. Local Arrays ..
294: INTEGER ISAVE( 3 )
295: DOUBLE PRECISION RWORK( 1 )
296: * ..
297: * .. External Functions ..
298: LOGICAL LSAME
299: DOUBLE PRECISION ZLANGE
300: EXTERNAL LSAME, ZLANGE
301: * ..
302: * .. External Subroutines ..
303: EXTERNAL XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL
304: * ..
305: * .. Intrinsic Functions ..
306: INTRINSIC MAX, SQRT
307: * ..
308: * .. Executable Statements ..
309: *
310: * Decode and test the input parameters.
311: *
312: WANTBH = LSAME( JOB, 'B' )
313: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
314: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
315: WANTQ = LSAME( COMPQ, 'V' )
316: *
317: * Set M to the number of selected eigenvalues.
318: *
319: M = 0
320: DO 10 K = 1, N
321: IF( SELECT( K ) )
322: $ M = M + 1
323: 10 CONTINUE
324: *
325: N1 = M
326: N2 = N - M
327: NN = N1*N2
328: *
329: INFO = 0
330: LQUERY = ( LWORK.EQ.-1 )
331: *
332: IF( WANTSP ) THEN
333: LWMIN = MAX( 1, 2*NN )
334: ELSE IF( LSAME( JOB, 'N' ) ) THEN
335: LWMIN = 1
336: ELSE IF( LSAME( JOB, 'E' ) ) THEN
337: LWMIN = MAX( 1, NN )
338: END IF
339: *
340: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
341: $ THEN
342: INFO = -1
343: ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
344: INFO = -2
345: ELSE IF( N.LT.0 ) THEN
346: INFO = -4
347: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
348: INFO = -6
349: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
350: INFO = -8
351: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
352: INFO = -14
353: END IF
354: *
355: IF( INFO.EQ.0 ) THEN
356: WORK( 1 ) = LWMIN
357: END IF
358: *
359: IF( INFO.NE.0 ) THEN
360: CALL XERBLA( 'ZTRSEN', -INFO )
361: RETURN
362: ELSE IF( LQUERY ) THEN
363: RETURN
364: END IF
365: *
366: * Quick return if possible
367: *
368: IF( M.EQ.N .OR. M.EQ.0 ) THEN
369: IF( WANTS )
370: $ S = ONE
371: IF( WANTSP )
372: $ SEP = ZLANGE( '1', N, N, T, LDT, RWORK )
373: GO TO 40
374: END IF
375: *
376: * Collect the selected eigenvalues at the top left corner of T.
377: *
378: KS = 0
379: DO 20 K = 1, N
380: IF( SELECT( K ) ) THEN
381: KS = KS + 1
382: *
383: * Swap the K-th eigenvalue to position KS.
384: *
385: IF( K.NE.KS )
386: $ CALL ZTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
387: END IF
388: 20 CONTINUE
389: *
390: IF( WANTS ) THEN
391: *
392: * Solve the Sylvester equation for R:
393: *
394: * T11*R - R*T22 = scale*T12
395: *
396: CALL ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
397: CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
398: $ LDT, WORK, N1, SCALE, IERR )
399: *
400: * Estimate the reciprocal of the condition number of the cluster
401: * of eigenvalues.
402: *
403: RNORM = ZLANGE( 'F', N1, N2, WORK, N1, RWORK )
404: IF( RNORM.EQ.ZERO ) THEN
405: S = ONE
406: ELSE
407: S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
408: $ SQRT( RNORM ) )
409: END IF
410: END IF
411: *
412: IF( WANTSP ) THEN
413: *
414: * Estimate sep(T11,T22).
415: *
416: EST = ZERO
417: KASE = 0
418: 30 CONTINUE
419: CALL ZLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE )
420: IF( KASE.NE.0 ) THEN
421: IF( KASE.EQ.1 ) THEN
422: *
423: * Solve T11*R - R*T22 = scale*X.
424: *
425: CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
426: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
427: $ IERR )
428: ELSE
429: *
1.8 bertrand 430: * Solve T11**H*R - R*T22**H = scale*X.
1.1 bertrand 431: *
432: CALL ZTRSYL( 'C', 'C', -1, N1, N2, T, LDT,
433: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
434: $ IERR )
435: END IF
436: GO TO 30
437: END IF
438: *
439: SEP = SCALE / EST
440: END IF
441: *
442: 40 CONTINUE
443: *
444: * Copy reordered eigenvalues to W.
445: *
446: DO 50 K = 1, N
447: W( K ) = T( K, K )
448: 50 CONTINUE
449: *
450: WORK( 1 ) = LWMIN
451: *
452: RETURN
453: *
454: * End of ZTRSEN
455: *
456: 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