1: SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
2: $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
3: $ INFO )
4: *
5: * -- LAPACK routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
11: *
12: * .. Scalar Arguments ..
13: CHARACTER HOWMNY, JOB
14: INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
15: * ..
16: * .. Array Arguments ..
17: LOGICAL SELECT( * )
18: DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
19: COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
20: $ WORK( LDWORK, * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * ZTRSNA estimates reciprocal condition numbers for specified
27: * eigenvalues and/or right eigenvectors of a complex upper triangular
28: * matrix T (or of any matrix Q*T*Q**H with Q unitary).
29: *
30: * Arguments
31: * =========
32: *
33: * JOB (input) CHARACTER*1
34: * Specifies whether condition numbers are required for
35: * eigenvalues (S) or eigenvectors (SEP):
36: * = 'E': for eigenvalues only (S);
37: * = 'V': for eigenvectors only (SEP);
38: * = 'B': for both eigenvalues and eigenvectors (S and SEP).
39: *
40: * HOWMNY (input) CHARACTER*1
41: * = 'A': compute condition numbers for all eigenpairs;
42: * = 'S': compute condition numbers for selected eigenpairs
43: * specified by the array SELECT.
44: *
45: * SELECT (input) LOGICAL array, dimension (N)
46: * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
47: * condition numbers are required. To select condition numbers
48: * for the j-th eigenpair, SELECT(j) must be set to .TRUE..
49: * If HOWMNY = 'A', SELECT is not referenced.
50: *
51: * N (input) INTEGER
52: * The order of the matrix T. N >= 0.
53: *
54: * T (input) COMPLEX*16 array, dimension (LDT,N)
55: * The upper triangular matrix T.
56: *
57: * LDT (input) INTEGER
58: * The leading dimension of the array T. LDT >= max(1,N).
59: *
60: * VL (input) COMPLEX*16 array, dimension (LDVL,M)
61: * If JOB = 'E' or 'B', VL must contain left eigenvectors of T
62: * (or of any Q*T*Q**H with Q unitary), corresponding to the
63: * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
64: * must be stored in consecutive columns of VL, as returned by
65: * ZHSEIN or ZTREVC.
66: * If JOB = 'V', VL is not referenced.
67: *
68: * LDVL (input) INTEGER
69: * The leading dimension of the array VL.
70: * LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
71: *
72: * VR (input) COMPLEX*16 array, dimension (LDVR,M)
73: * If JOB = 'E' or 'B', VR must contain right eigenvectors of T
74: * (or of any Q*T*Q**H with Q unitary), corresponding to the
75: * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
76: * must be stored in consecutive columns of VR, as returned by
77: * ZHSEIN or ZTREVC.
78: * If JOB = 'V', VR is not referenced.
79: *
80: * LDVR (input) INTEGER
81: * The leading dimension of the array VR.
82: * LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
83: *
84: * S (output) DOUBLE PRECISION array, dimension (MM)
85: * If JOB = 'E' or 'B', the reciprocal condition numbers of the
86: * selected eigenvalues, stored in consecutive elements of the
87: * array. Thus S(j), SEP(j), and the j-th columns of VL and VR
88: * all correspond to the same eigenpair (but not in general the
89: * j-th eigenpair, unless all eigenpairs are selected).
90: * If JOB = 'V', S is not referenced.
91: *
92: * SEP (output) DOUBLE PRECISION array, dimension (MM)
93: * If JOB = 'V' or 'B', the estimated reciprocal condition
94: * numbers of the selected eigenvectors, stored in consecutive
95: * elements of the array.
96: * If JOB = 'E', SEP is not referenced.
97: *
98: * MM (input) INTEGER
99: * The number of elements in the arrays S (if JOB = 'E' or 'B')
100: * and/or SEP (if JOB = 'V' or 'B'). MM >= M.
101: *
102: * M (output) INTEGER
103: * The number of elements of the arrays S and/or SEP actually
104: * used to store the estimated condition numbers.
105: * If HOWMNY = 'A', M is set to N.
106: *
107: * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N+6)
108: * If JOB = 'E', WORK is not referenced.
109: *
110: * LDWORK (input) INTEGER
111: * The leading dimension of the array WORK.
112: * LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
113: *
114: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
115: * If JOB = 'E', RWORK is not referenced.
116: *
117: * INFO (output) INTEGER
118: * = 0: successful exit
119: * < 0: if INFO = -i, the i-th argument had an illegal value
120: *
121: * Further Details
122: * ===============
123: *
124: * The reciprocal of the condition number of an eigenvalue lambda is
125: * defined as
126: *
127: * S(lambda) = |v'*u| / (norm(u)*norm(v))
128: *
129: * where u and v are the right and left eigenvectors of T corresponding
130: * to lambda; v' denotes the conjugate transpose of v, and norm(u)
131: * denotes the Euclidean norm. These reciprocal condition numbers always
132: * lie between zero (very badly conditioned) and one (very well
133: * conditioned). If n = 1, S(lambda) is defined to be 1.
134: *
135: * An approximate error bound for a computed eigenvalue W(i) is given by
136: *
137: * EPS * norm(T) / S(i)
138: *
139: * where EPS is the machine precision.
140: *
141: * The reciprocal of the condition number of the right eigenvector u
142: * corresponding to lambda is defined as follows. Suppose
143: *
144: * T = ( lambda c )
145: * ( 0 T22 )
146: *
147: * Then the reciprocal condition number is
148: *
149: * SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
150: *
151: * where sigma-min denotes the smallest singular value. We approximate
152: * the smallest singular value by the reciprocal of an estimate of the
153: * one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
154: * defined to be abs(T(1,1)).
155: *
156: * An approximate error bound for a computed right eigenvector VR(i)
157: * is given by
158: *
159: * EPS * norm(T) / SEP(i)
160: *
161: * =====================================================================
162: *
163: * .. Parameters ..
164: DOUBLE PRECISION ZERO, ONE
165: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
166: * ..
167: * .. Local Scalars ..
168: LOGICAL SOMCON, WANTBH, WANTS, WANTSP
169: CHARACTER NORMIN
170: INTEGER I, IERR, IX, J, K, KASE, KS
171: DOUBLE PRECISION BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
172: $ XNORM
173: COMPLEX*16 CDUM, PROD
174: * ..
175: * .. Local Arrays ..
176: INTEGER ISAVE( 3 )
177: COMPLEX*16 DUMMY( 1 )
178: * ..
179: * .. External Functions ..
180: LOGICAL LSAME
181: INTEGER IZAMAX
182: DOUBLE PRECISION DLAMCH, DZNRM2
183: COMPLEX*16 ZDOTC
184: EXTERNAL LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
185: * ..
186: * .. External Subroutines ..
187: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC
188: * ..
189: * .. Intrinsic Functions ..
190: INTRINSIC ABS, DBLE, DIMAG, MAX
191: * ..
192: * .. Statement Functions ..
193: DOUBLE PRECISION CABS1
194: * ..
195: * .. Statement Function definitions ..
196: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
197: * ..
198: * .. Executable Statements ..
199: *
200: * Decode and test the input parameters
201: *
202: WANTBH = LSAME( JOB, 'B' )
203: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
204: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
205: *
206: SOMCON = LSAME( HOWMNY, 'S' )
207: *
208: * Set M to the number of eigenpairs for which condition numbers are
209: * to be computed.
210: *
211: IF( SOMCON ) THEN
212: M = 0
213: DO 10 J = 1, N
214: IF( SELECT( J ) )
215: $ M = M + 1
216: 10 CONTINUE
217: ELSE
218: M = N
219: END IF
220: *
221: INFO = 0
222: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
223: INFO = -1
224: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
225: INFO = -2
226: ELSE IF( N.LT.0 ) THEN
227: INFO = -4
228: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
229: INFO = -6
230: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
231: INFO = -8
232: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
233: INFO = -10
234: ELSE IF( MM.LT.M ) THEN
235: INFO = -13
236: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
237: INFO = -16
238: END IF
239: IF( INFO.NE.0 ) THEN
240: CALL XERBLA( 'ZTRSNA', -INFO )
241: RETURN
242: END IF
243: *
244: * Quick return if possible
245: *
246: IF( N.EQ.0 )
247: $ RETURN
248: *
249: IF( N.EQ.1 ) THEN
250: IF( SOMCON ) THEN
251: IF( .NOT.SELECT( 1 ) )
252: $ RETURN
253: END IF
254: IF( WANTS )
255: $ S( 1 ) = ONE
256: IF( WANTSP )
257: $ SEP( 1 ) = ABS( T( 1, 1 ) )
258: RETURN
259: END IF
260: *
261: * Get machine constants
262: *
263: EPS = DLAMCH( 'P' )
264: SMLNUM = DLAMCH( 'S' ) / EPS
265: BIGNUM = ONE / SMLNUM
266: CALL DLABAD( SMLNUM, BIGNUM )
267: *
268: KS = 1
269: DO 50 K = 1, N
270: *
271: IF( SOMCON ) THEN
272: IF( .NOT.SELECT( K ) )
273: $ GO TO 50
274: END IF
275: *
276: IF( WANTS ) THEN
277: *
278: * Compute the reciprocal condition number of the k-th
279: * eigenvalue.
280: *
281: PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
282: RNRM = DZNRM2( N, VR( 1, KS ), 1 )
283: LNRM = DZNRM2( N, VL( 1, KS ), 1 )
284: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
285: *
286: END IF
287: *
288: IF( WANTSP ) THEN
289: *
290: * Estimate the reciprocal condition number of the k-th
291: * eigenvector.
292: *
293: * Copy the matrix T to the array WORK and swap the k-th
294: * diagonal element to the (1,1) position.
295: *
296: CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
297: CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
298: *
299: * Form C = T22 - lambda*I in WORK(2:N,2:N).
300: *
301: DO 20 I = 2, N
302: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
303: 20 CONTINUE
304: *
305: * Estimate a lower bound for the 1-norm of inv(C'). The 1st
306: * and (N+1)th columns of WORK are used to store work vectors.
307: *
308: SEP( KS ) = ZERO
309: EST = ZERO
310: KASE = 0
311: NORMIN = 'N'
312: 30 CONTINUE
313: CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
314: *
315: IF( KASE.NE.0 ) THEN
316: IF( KASE.EQ.1 ) THEN
317: *
318: * Solve C'*x = scale*b
319: *
320: CALL ZLATRS( 'Upper', 'Conjugate transpose',
321: $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
322: $ LDWORK, WORK, SCALE, RWORK, IERR )
323: ELSE
324: *
325: * Solve C*x = scale*b
326: *
327: CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
328: $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
329: $ SCALE, RWORK, IERR )
330: END IF
331: NORMIN = 'Y'
332: IF( SCALE.NE.ONE ) THEN
333: *
334: * Multiply by 1/SCALE if doing so will not cause
335: * overflow.
336: *
337: IX = IZAMAX( N-1, WORK, 1 )
338: XNORM = CABS1( WORK( IX, 1 ) )
339: IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
340: $ GO TO 40
341: CALL ZDRSCL( N, SCALE, WORK, 1 )
342: END IF
343: GO TO 30
344: END IF
345: *
346: SEP( KS ) = ONE / MAX( EST, SMLNUM )
347: END IF
348: *
349: 40 CONTINUE
350: KS = KS + 1
351: 50 CONTINUE
352: RETURN
353: *
354: * End of ZTRSNA
355: *
356: 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