Annotation of rpl/lapack/lapack/ztrsen.f, revision 1.9
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