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: * =====================================================================
248: SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
249: $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
250: $ INFO )
251: *
252: * -- LAPACK computational routine (version 3.4.0) --
253: * -- LAPACK is a software package provided by Univ. of Tennessee, --
254: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
255: * November 2011
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: *
412: * Estimate a lower bound for the 1-norm of inv(C**H). The 1st
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: *
425: * Solve C**H*x = scale*b
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