Annotation of rpl/lapack/lapack/ztrsna.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZTRSNA
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZTRSNA + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsna.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsna.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsna.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
! 22: * LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
! 23: * INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER HOWMNY, JOB
! 27: * INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
! 28: * ..
! 29: * .. Array Arguments ..
! 30: * LOGICAL SELECT( * )
! 31: * DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
! 32: * COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
! 33: * $ WORK( LDWORK, * )
! 34: * ..
! 35: *
! 36: *
! 37: *> \par Purpose:
! 38: * =============
! 39: *>
! 40: *> \verbatim
! 41: *>
! 42: *> ZTRSNA estimates reciprocal condition numbers for specified
! 43: *> eigenvalues and/or right eigenvectors of a complex upper triangular
! 44: *> matrix T (or of any matrix Q*T*Q**H with Q unitary).
! 45: *> \endverbatim
! 46: *
! 47: * Arguments:
! 48: * ==========
! 49: *
! 50: *> \param[in] JOB
! 51: *> \verbatim
! 52: *> JOB is CHARACTER*1
! 53: *> Specifies whether condition numbers are required for
! 54: *> eigenvalues (S) or eigenvectors (SEP):
! 55: *> = 'E': for eigenvalues only (S);
! 56: *> = 'V': for eigenvectors only (SEP);
! 57: *> = 'B': for both eigenvalues and eigenvectors (S and SEP).
! 58: *> \endverbatim
! 59: *>
! 60: *> \param[in] HOWMNY
! 61: *> \verbatim
! 62: *> HOWMNY is CHARACTER*1
! 63: *> = 'A': compute condition numbers for all eigenpairs;
! 64: *> = 'S': compute condition numbers for selected eigenpairs
! 65: *> specified by the array SELECT.
! 66: *> \endverbatim
! 67: *>
! 68: *> \param[in] SELECT
! 69: *> \verbatim
! 70: *> SELECT is LOGICAL array, dimension (N)
! 71: *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
! 72: *> condition numbers are required. To select condition numbers
! 73: *> for the j-th eigenpair, SELECT(j) must be set to .TRUE..
! 74: *> If HOWMNY = 'A', SELECT is not referenced.
! 75: *> \endverbatim
! 76: *>
! 77: *> \param[in] N
! 78: *> \verbatim
! 79: *> N is INTEGER
! 80: *> The order of the matrix T. N >= 0.
! 81: *> \endverbatim
! 82: *>
! 83: *> \param[in] T
! 84: *> \verbatim
! 85: *> T is COMPLEX*16 array, dimension (LDT,N)
! 86: *> The upper triangular matrix T.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] LDT
! 90: *> \verbatim
! 91: *> LDT is INTEGER
! 92: *> The leading dimension of the array T. LDT >= max(1,N).
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[in] VL
! 96: *> \verbatim
! 97: *> VL is COMPLEX*16 array, dimension (LDVL,M)
! 98: *> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
! 99: *> (or of any Q*T*Q**H with Q unitary), corresponding to the
! 100: *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
! 101: *> must be stored in consecutive columns of VL, as returned by
! 102: *> ZHSEIN or ZTREVC.
! 103: *> If JOB = 'V', VL is not referenced.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[in] LDVL
! 107: *> \verbatim
! 108: *> LDVL is INTEGER
! 109: *> The leading dimension of the array VL.
! 110: *> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in] VR
! 114: *> \verbatim
! 115: *> VR is COMPLEX*16 array, dimension (LDVR,M)
! 116: *> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
! 117: *> (or of any Q*T*Q**H with Q unitary), corresponding to the
! 118: *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
! 119: *> must be stored in consecutive columns of VR, as returned by
! 120: *> ZHSEIN or ZTREVC.
! 121: *> If JOB = 'V', VR is not referenced.
! 122: *> \endverbatim
! 123: *>
! 124: *> \param[in] LDVR
! 125: *> \verbatim
! 126: *> LDVR is INTEGER
! 127: *> The leading dimension of the array VR.
! 128: *> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
! 129: *> \endverbatim
! 130: *>
! 131: *> \param[out] S
! 132: *> \verbatim
! 133: *> S is DOUBLE PRECISION array, dimension (MM)
! 134: *> If JOB = 'E' or 'B', the reciprocal condition numbers of the
! 135: *> selected eigenvalues, stored in consecutive elements of the
! 136: *> array. Thus S(j), SEP(j), and the j-th columns of VL and VR
! 137: *> all correspond to the same eigenpair (but not in general the
! 138: *> j-th eigenpair, unless all eigenpairs are selected).
! 139: *> If JOB = 'V', S is not referenced.
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[out] SEP
! 143: *> \verbatim
! 144: *> SEP is DOUBLE PRECISION array, dimension (MM)
! 145: *> If JOB = 'V' or 'B', the estimated reciprocal condition
! 146: *> numbers of the selected eigenvectors, stored in consecutive
! 147: *> elements of the array.
! 148: *> If JOB = 'E', SEP is not referenced.
! 149: *> \endverbatim
! 150: *>
! 151: *> \param[in] MM
! 152: *> \verbatim
! 153: *> MM is INTEGER
! 154: *> The number of elements in the arrays S (if JOB = 'E' or 'B')
! 155: *> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
! 156: *> \endverbatim
! 157: *>
! 158: *> \param[out] M
! 159: *> \verbatim
! 160: *> M is INTEGER
! 161: *> The number of elements of the arrays S and/or SEP actually
! 162: *> used to store the estimated condition numbers.
! 163: *> If HOWMNY = 'A', M is set to N.
! 164: *> \endverbatim
! 165: *>
! 166: *> \param[out] WORK
! 167: *> \verbatim
! 168: *> WORK is COMPLEX*16 array, dimension (LDWORK,N+6)
! 169: *> If JOB = 'E', WORK is not referenced.
! 170: *> \endverbatim
! 171: *>
! 172: *> \param[in] LDWORK
! 173: *> \verbatim
! 174: *> LDWORK is INTEGER
! 175: *> The leading dimension of the array WORK.
! 176: *> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
! 177: *> \endverbatim
! 178: *>
! 179: *> \param[out] RWORK
! 180: *> \verbatim
! 181: *> RWORK is DOUBLE PRECISION array, dimension (N)
! 182: *> If JOB = 'E', RWORK is not referenced.
! 183: *> \endverbatim
! 184: *>
! 185: *> \param[out] INFO
! 186: *> \verbatim
! 187: *> INFO is INTEGER
! 188: *> = 0: successful exit
! 189: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 190: *> \endverbatim
! 191: *
! 192: * Authors:
! 193: * ========
! 194: *
! 195: *> \author Univ. of Tennessee
! 196: *> \author Univ. of California Berkeley
! 197: *> \author Univ. of Colorado Denver
! 198: *> \author NAG Ltd.
! 199: *
! 200: *> \date November 2011
! 201: *
! 202: *> \ingroup complex16OTHERcomputational
! 203: *
! 204: *> \par Further Details:
! 205: * =====================
! 206: *>
! 207: *> \verbatim
! 208: *>
! 209: *> The reciprocal of the condition number of an eigenvalue lambda is
! 210: *> defined as
! 211: *>
! 212: *> S(lambda) = |v**H*u| / (norm(u)*norm(v))
! 213: *>
! 214: *> where u and v are the right and left eigenvectors of T corresponding
! 215: *> to lambda; v**H denotes the conjugate transpose of v, and norm(u)
! 216: *> denotes the Euclidean norm. These reciprocal condition numbers always
! 217: *> lie between zero (very badly conditioned) and one (very well
! 218: *> conditioned). If n = 1, S(lambda) is defined to be 1.
! 219: *>
! 220: *> An approximate error bound for a computed eigenvalue W(i) is given by
! 221: *>
! 222: *> EPS * norm(T) / S(i)
! 223: *>
! 224: *> where EPS is the machine precision.
! 225: *>
! 226: *> The reciprocal of the condition number of the right eigenvector u
! 227: *> corresponding to lambda is defined as follows. Suppose
! 228: *>
! 229: *> T = ( lambda c )
! 230: *> ( 0 T22 )
! 231: *>
! 232: *> Then the reciprocal condition number is
! 233: *>
! 234: *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
! 235: *>
! 236: *> where sigma-min denotes the smallest singular value. We approximate
! 237: *> the smallest singular value by the reciprocal of an estimate of the
! 238: *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
! 239: *> defined to be abs(T(1,1)).
! 240: *>
! 241: *> An approximate error bound for a computed right eigenvector VR(i)
! 242: *> is given by
! 243: *>
! 244: *> EPS * norm(T) / SEP(i)
! 245: *> \endverbatim
! 246: *>
! 247: * =====================================================================
1.1 bertrand 248: SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
249: $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
250: $ INFO )
251: *
1.9 ! bertrand 252: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 253: * -- LAPACK is a software package provided by Univ. of Tennessee, --
254: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 255: * November 2011
1.1 bertrand 256: *
257: * .. Scalar Arguments ..
258: CHARACTER HOWMNY, JOB
259: INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
260: * ..
261: * .. Array Arguments ..
262: LOGICAL SELECT( * )
263: DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
264: COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
265: $ WORK( LDWORK, * )
266: * ..
267: *
268: * =====================================================================
269: *
270: * .. Parameters ..
271: DOUBLE PRECISION ZERO, ONE
272: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
273: * ..
274: * .. Local Scalars ..
275: LOGICAL SOMCON, WANTBH, WANTS, WANTSP
276: CHARACTER NORMIN
277: INTEGER I, IERR, IX, J, K, KASE, KS
278: DOUBLE PRECISION BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
279: $ XNORM
280: COMPLEX*16 CDUM, PROD
281: * ..
282: * .. Local Arrays ..
283: INTEGER ISAVE( 3 )
284: COMPLEX*16 DUMMY( 1 )
285: * ..
286: * .. External Functions ..
287: LOGICAL LSAME
288: INTEGER IZAMAX
289: DOUBLE PRECISION DLAMCH, DZNRM2
290: COMPLEX*16 ZDOTC
291: EXTERNAL LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
292: * ..
293: * .. External Subroutines ..
294: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC
295: * ..
296: * .. Intrinsic Functions ..
297: INTRINSIC ABS, DBLE, DIMAG, MAX
298: * ..
299: * .. Statement Functions ..
300: DOUBLE PRECISION CABS1
301: * ..
302: * .. Statement Function definitions ..
303: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
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: *
313: SOMCON = LSAME( HOWMNY, 'S' )
314: *
315: * Set M to the number of eigenpairs for which condition numbers are
316: * to be computed.
317: *
318: IF( SOMCON ) THEN
319: M = 0
320: DO 10 J = 1, N
321: IF( SELECT( J ) )
322: $ M = M + 1
323: 10 CONTINUE
324: ELSE
325: M = N
326: END IF
327: *
328: INFO = 0
329: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
330: INFO = -1
331: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
332: INFO = -2
333: ELSE IF( N.LT.0 ) THEN
334: INFO = -4
335: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
336: INFO = -6
337: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
338: INFO = -8
339: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
340: INFO = -10
341: ELSE IF( MM.LT.M ) THEN
342: INFO = -13
343: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
344: INFO = -16
345: END IF
346: IF( INFO.NE.0 ) THEN
347: CALL XERBLA( 'ZTRSNA', -INFO )
348: RETURN
349: END IF
350: *
351: * Quick return if possible
352: *
353: IF( N.EQ.0 )
354: $ RETURN
355: *
356: IF( N.EQ.1 ) THEN
357: IF( SOMCON ) THEN
358: IF( .NOT.SELECT( 1 ) )
359: $ RETURN
360: END IF
361: IF( WANTS )
362: $ S( 1 ) = ONE
363: IF( WANTSP )
364: $ SEP( 1 ) = ABS( T( 1, 1 ) )
365: RETURN
366: END IF
367: *
368: * Get machine constants
369: *
370: EPS = DLAMCH( 'P' )
371: SMLNUM = DLAMCH( 'S' ) / EPS
372: BIGNUM = ONE / SMLNUM
373: CALL DLABAD( SMLNUM, BIGNUM )
374: *
375: KS = 1
376: DO 50 K = 1, N
377: *
378: IF( SOMCON ) THEN
379: IF( .NOT.SELECT( K ) )
380: $ GO TO 50
381: END IF
382: *
383: IF( WANTS ) THEN
384: *
385: * Compute the reciprocal condition number of the k-th
386: * eigenvalue.
387: *
388: PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
389: RNRM = DZNRM2( N, VR( 1, KS ), 1 )
390: LNRM = DZNRM2( N, VL( 1, KS ), 1 )
391: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
392: *
393: END IF
394: *
395: IF( WANTSP ) THEN
396: *
397: * Estimate the reciprocal condition number of the k-th
398: * eigenvector.
399: *
400: * Copy the matrix T to the array WORK and swap the k-th
401: * diagonal element to the (1,1) position.
402: *
403: CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
404: CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
405: *
406: * Form C = T22 - lambda*I in WORK(2:N,2:N).
407: *
408: DO 20 I = 2, N
409: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
410: 20 CONTINUE
411: *
1.8 bertrand 412: * Estimate a lower bound for the 1-norm of inv(C**H). The 1st
1.1 bertrand 413: * and (N+1)th columns of WORK are used to store work vectors.
414: *
415: SEP( KS ) = ZERO
416: EST = ZERO
417: KASE = 0
418: NORMIN = 'N'
419: 30 CONTINUE
420: CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
421: *
422: IF( KASE.NE.0 ) THEN
423: IF( KASE.EQ.1 ) THEN
424: *
1.8 bertrand 425: * Solve C**H*x = scale*b
1.1 bertrand 426: *
427: CALL ZLATRS( 'Upper', 'Conjugate transpose',
428: $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
429: $ LDWORK, WORK, SCALE, RWORK, IERR )
430: ELSE
431: *
432: * Solve C*x = scale*b
433: *
434: CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
435: $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
436: $ SCALE, RWORK, IERR )
437: END IF
438: NORMIN = 'Y'
439: IF( SCALE.NE.ONE ) THEN
440: *
441: * Multiply by 1/SCALE if doing so will not cause
442: * overflow.
443: *
444: IX = IZAMAX( N-1, WORK, 1 )
445: XNORM = CABS1( WORK( IX, 1 ) )
446: IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
447: $ GO TO 40
448: CALL ZDRSCL( N, SCALE, WORK, 1 )
449: END IF
450: GO TO 30
451: END IF
452: *
453: SEP( KS ) = ONE / MAX( EST, SMLNUM )
454: END IF
455: *
456: 40 CONTINUE
457: KS = KS + 1
458: 50 CONTINUE
459: RETURN
460: *
461: * End of ZTRSNA
462: *
463: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztrsna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack