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: *> \ingroup complex16OTHERcomputational
201: *
202: *> \par Further Details:
203: * =====================
204: *>
205: *> \verbatim
206: *>
207: *> The reciprocal of the condition number of an eigenvalue lambda is
208: *> defined as
209: *>
210: *> S(lambda) = |v**H*u| / (norm(u)*norm(v))
211: *>
212: *> where u and v are the right and left eigenvectors of T corresponding
213: *> to lambda; v**H denotes the conjugate transpose of v, and norm(u)
214: *> denotes the Euclidean norm. These reciprocal condition numbers always
215: *> lie between zero (very badly conditioned) and one (very well
216: *> conditioned). If n = 1, S(lambda) is defined to be 1.
217: *>
218: *> An approximate error bound for a computed eigenvalue W(i) is given by
219: *>
220: *> EPS * norm(T) / S(i)
221: *>
222: *> where EPS is the machine precision.
223: *>
224: *> The reciprocal of the condition number of the right eigenvector u
225: *> corresponding to lambda is defined as follows. Suppose
226: *>
227: *> T = ( lambda c )
228: *> ( 0 T22 )
229: *>
230: *> Then the reciprocal condition number is
231: *>
232: *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
233: *>
234: *> where sigma-min denotes the smallest singular value. We approximate
235: *> the smallest singular value by the reciprocal of an estimate of the
236: *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
237: *> defined to be abs(T(1,1)).
238: *>
239: *> An approximate error bound for a computed right eigenvector VR(i)
240: *> is given by
241: *>
242: *> EPS * norm(T) / SEP(i)
243: *> \endverbatim
244: *>
245: * =====================================================================
246: SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
247: $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
248: $ INFO )
249: *
250: * -- LAPACK computational routine --
251: * -- LAPACK is a software package provided by Univ. of Tennessee, --
252: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253: *
254: * .. Scalar Arguments ..
255: CHARACTER HOWMNY, JOB
256: INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
257: * ..
258: * .. Array Arguments ..
259: LOGICAL SELECT( * )
260: DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
261: COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
262: $ WORK( LDWORK, * )
263: * ..
264: *
265: * =====================================================================
266: *
267: * .. Parameters ..
268: DOUBLE PRECISION ZERO, ONE
269: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
270: * ..
271: * .. Local Scalars ..
272: LOGICAL SOMCON, WANTBH, WANTS, WANTSP
273: CHARACTER NORMIN
274: INTEGER I, IERR, IX, J, K, KASE, KS
275: DOUBLE PRECISION BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
276: $ XNORM
277: COMPLEX*16 CDUM, PROD
278: * ..
279: * .. Local Arrays ..
280: INTEGER ISAVE( 3 )
281: COMPLEX*16 DUMMY( 1 )
282: * ..
283: * .. External Functions ..
284: LOGICAL LSAME
285: INTEGER IZAMAX
286: DOUBLE PRECISION DLAMCH, DZNRM2
287: COMPLEX*16 ZDOTC
288: EXTERNAL LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
289: * ..
290: * .. External Subroutines ..
291: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC,
292: $ DLABAD
293: * ..
294: * .. Intrinsic Functions ..
295: INTRINSIC ABS, DBLE, DIMAG, MAX
296: * ..
297: * .. Statement Functions ..
298: DOUBLE PRECISION CABS1
299: * ..
300: * .. Statement Function definitions ..
301: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
302: * ..
303: * .. Executable Statements ..
304: *
305: * Decode and test the input parameters
306: *
307: WANTBH = LSAME( JOB, 'B' )
308: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
309: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
310: *
311: SOMCON = LSAME( HOWMNY, 'S' )
312: *
313: * Set M to the number of eigenpairs for which condition numbers are
314: * to be computed.
315: *
316: IF( SOMCON ) THEN
317: M = 0
318: DO 10 J = 1, N
319: IF( SELECT( J ) )
320: $ M = M + 1
321: 10 CONTINUE
322: ELSE
323: M = N
324: END IF
325: *
326: INFO = 0
327: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
328: INFO = -1
329: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
330: INFO = -2
331: ELSE IF( N.LT.0 ) THEN
332: INFO = -4
333: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
334: INFO = -6
335: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
336: INFO = -8
337: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
338: INFO = -10
339: ELSE IF( MM.LT.M ) THEN
340: INFO = -13
341: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
342: INFO = -16
343: END IF
344: IF( INFO.NE.0 ) THEN
345: CALL XERBLA( 'ZTRSNA', -INFO )
346: RETURN
347: END IF
348: *
349: * Quick return if possible
350: *
351: IF( N.EQ.0 )
352: $ RETURN
353: *
354: IF( N.EQ.1 ) THEN
355: IF( SOMCON ) THEN
356: IF( .NOT.SELECT( 1 ) )
357: $ RETURN
358: END IF
359: IF( WANTS )
360: $ S( 1 ) = ONE
361: IF( WANTSP )
362: $ SEP( 1 ) = ABS( T( 1, 1 ) )
363: RETURN
364: END IF
365: *
366: * Get machine constants
367: *
368: EPS = DLAMCH( 'P' )
369: SMLNUM = DLAMCH( 'S' ) / EPS
370: BIGNUM = ONE / SMLNUM
371: CALL DLABAD( SMLNUM, BIGNUM )
372: *
373: KS = 1
374: DO 50 K = 1, N
375: *
376: IF( SOMCON ) THEN
377: IF( .NOT.SELECT( K ) )
378: $ GO TO 50
379: END IF
380: *
381: IF( WANTS ) THEN
382: *
383: * Compute the reciprocal condition number of the k-th
384: * eigenvalue.
385: *
386: PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
387: RNRM = DZNRM2( N, VR( 1, KS ), 1 )
388: LNRM = DZNRM2( N, VL( 1, KS ), 1 )
389: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
390: *
391: END IF
392: *
393: IF( WANTSP ) THEN
394: *
395: * Estimate the reciprocal condition number of the k-th
396: * eigenvector.
397: *
398: * Copy the matrix T to the array WORK and swap the k-th
399: * diagonal element to the (1,1) position.
400: *
401: CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
402: CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
403: *
404: * Form C = T22 - lambda*I in WORK(2:N,2:N).
405: *
406: DO 20 I = 2, N
407: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
408: 20 CONTINUE
409: *
410: * Estimate a lower bound for the 1-norm of inv(C**H). The 1st
411: * and (N+1)th columns of WORK are used to store work vectors.
412: *
413: SEP( KS ) = ZERO
414: EST = ZERO
415: KASE = 0
416: NORMIN = 'N'
417: 30 CONTINUE
418: CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
419: *
420: IF( KASE.NE.0 ) THEN
421: IF( KASE.EQ.1 ) THEN
422: *
423: * Solve C**H*x = scale*b
424: *
425: CALL ZLATRS( 'Upper', 'Conjugate transpose',
426: $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
427: $ LDWORK, WORK, SCALE, RWORK, IERR )
428: ELSE
429: *
430: * Solve C*x = scale*b
431: *
432: CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
433: $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
434: $ SCALE, RWORK, IERR )
435: END IF
436: NORMIN = 'Y'
437: IF( SCALE.NE.ONE ) THEN
438: *
439: * Multiply by 1/SCALE if doing so will not cause
440: * overflow.
441: *
442: IX = IZAMAX( N-1, WORK, 1 )
443: XNORM = CABS1( WORK( IX, 1 ) )
444: IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
445: $ GO TO 40
446: CALL ZDRSCL( N, SCALE, WORK, 1 )
447: END IF
448: GO TO 30
449: END IF
450: *
451: SEP( KS ) = ONE / MAX( EST, SMLNUM )
452: END IF
453: *
454: 40 CONTINUE
455: KS = KS + 1
456: 50 CONTINUE
457: RETURN
458: *
459: * End of ZTRSNA
460: *
461: 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