Annotation of rpl/lapack/lapack/ztrsna.f, revision 1.17
1.9 bertrand 1: *> \brief \b ZTRSNA
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.9 bertrand 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 )
1.15 bertrand 24: *
1.9 bertrand 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: * ..
1.15 bertrand 35: *
1.9 bertrand 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: *
1.15 bertrand 195: *> \author Univ. of Tennessee
196: *> \author Univ. of California Berkeley
197: *> \author Univ. of Colorado Denver
198: *> \author NAG Ltd.
1.9 bertrand 199: *
1.17 ! bertrand 200: *> \date November 2017
1.9 bertrand 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.17 ! bertrand 252: * -- LAPACK computational routine (version 3.8.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.17 ! bertrand 255: * November 2017
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 ..
1.17 ! bertrand 294: EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC,
! 295: $ DLABAD
1.1 bertrand 296: * ..
297: * .. Intrinsic Functions ..
298: INTRINSIC ABS, DBLE, DIMAG, MAX
299: * ..
300: * .. Statement Functions ..
301: DOUBLE PRECISION CABS1
302: * ..
303: * .. Statement Function definitions ..
304: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
305: * ..
306: * .. Executable Statements ..
307: *
308: * Decode and test the input parameters
309: *
310: WANTBH = LSAME( JOB, 'B' )
311: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
312: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
313: *
314: SOMCON = LSAME( HOWMNY, 'S' )
315: *
316: * Set M to the number of eigenpairs for which condition numbers are
317: * to be computed.
318: *
319: IF( SOMCON ) THEN
320: M = 0
321: DO 10 J = 1, N
322: IF( SELECT( J ) )
323: $ M = M + 1
324: 10 CONTINUE
325: ELSE
326: M = N
327: END IF
328: *
329: INFO = 0
330: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
331: INFO = -1
332: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
333: INFO = -2
334: ELSE IF( N.LT.0 ) THEN
335: INFO = -4
336: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
337: INFO = -6
338: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
339: INFO = -8
340: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
341: INFO = -10
342: ELSE IF( MM.LT.M ) THEN
343: INFO = -13
344: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
345: INFO = -16
346: END IF
347: IF( INFO.NE.0 ) THEN
348: CALL XERBLA( 'ZTRSNA', -INFO )
349: RETURN
350: END IF
351: *
352: * Quick return if possible
353: *
354: IF( N.EQ.0 )
355: $ RETURN
356: *
357: IF( N.EQ.1 ) THEN
358: IF( SOMCON ) THEN
359: IF( .NOT.SELECT( 1 ) )
360: $ RETURN
361: END IF
362: IF( WANTS )
363: $ S( 1 ) = ONE
364: IF( WANTSP )
365: $ SEP( 1 ) = ABS( T( 1, 1 ) )
366: RETURN
367: END IF
368: *
369: * Get machine constants
370: *
371: EPS = DLAMCH( 'P' )
372: SMLNUM = DLAMCH( 'S' ) / EPS
373: BIGNUM = ONE / SMLNUM
374: CALL DLABAD( SMLNUM, BIGNUM )
375: *
376: KS = 1
377: DO 50 K = 1, N
378: *
379: IF( SOMCON ) THEN
380: IF( .NOT.SELECT( K ) )
381: $ GO TO 50
382: END IF
383: *
384: IF( WANTS ) THEN
385: *
386: * Compute the reciprocal condition number of the k-th
387: * eigenvalue.
388: *
389: PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
390: RNRM = DZNRM2( N, VR( 1, KS ), 1 )
391: LNRM = DZNRM2( N, VL( 1, KS ), 1 )
392: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
393: *
394: END IF
395: *
396: IF( WANTSP ) THEN
397: *
398: * Estimate the reciprocal condition number of the k-th
399: * eigenvector.
400: *
401: * Copy the matrix T to the array WORK and swap the k-th
402: * diagonal element to the (1,1) position.
403: *
404: CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
405: CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
406: *
407: * Form C = T22 - lambda*I in WORK(2:N,2:N).
408: *
409: DO 20 I = 2, N
410: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
411: 20 CONTINUE
412: *
1.8 bertrand 413: * Estimate a lower bound for the 1-norm of inv(C**H). The 1st
1.1 bertrand 414: * and (N+1)th columns of WORK are used to store work vectors.
415: *
416: SEP( KS ) = ZERO
417: EST = ZERO
418: KASE = 0
419: NORMIN = 'N'
420: 30 CONTINUE
421: CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
422: *
423: IF( KASE.NE.0 ) THEN
424: IF( KASE.EQ.1 ) THEN
425: *
1.8 bertrand 426: * Solve C**H*x = scale*b
1.1 bertrand 427: *
428: CALL ZLATRS( 'Upper', 'Conjugate transpose',
429: $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
430: $ LDWORK, WORK, SCALE, RWORK, IERR )
431: ELSE
432: *
433: * Solve C*x = scale*b
434: *
435: CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
436: $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
437: $ SCALE, RWORK, IERR )
438: END IF
439: NORMIN = 'Y'
440: IF( SCALE.NE.ONE ) THEN
441: *
442: * Multiply by 1/SCALE if doing so will not cause
443: * overflow.
444: *
445: IX = IZAMAX( N-1, WORK, 1 )
446: XNORM = CABS1( WORK( IX, 1 ) )
447: IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
448: $ GO TO 40
449: CALL ZDRSCL( N, SCALE, WORK, 1 )
450: END IF
451: GO TO 30
452: END IF
453: *
454: SEP( KS ) = ONE / MAX( EST, SMLNUM )
455: END IF
456: *
457: 40 CONTINUE
458: KS = KS + 1
459: 50 CONTINUE
460: RETURN
461: *
462: * End of ZTRSNA
463: *
464: 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