1: *> \brief \b DTRSNA
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTRSNA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22: * LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
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: * INTEGER IWORK( * )
32: * DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
33: * $ VR( LDVR, * ), WORK( LDWORK, * )
34: * ..
35: *
36: *
37: *> \par Purpose:
38: * =============
39: *>
40: *> \verbatim
41: *>
42: *> DTRSNA estimates reciprocal condition numbers for specified
43: *> eigenvalues and/or right eigenvectors of a real upper
44: *> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
45: *> orthogonal).
46: *>
47: *> T must be in Schur canonical form (as returned by DHSEQR), that is,
48: *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
49: *> 2-by-2 diagonal block has its diagonal elements equal and its
50: *> off-diagonal elements of opposite sign.
51: *> \endverbatim
52: *
53: * Arguments:
54: * ==========
55: *
56: *> \param[in] JOB
57: *> \verbatim
58: *> JOB is CHARACTER*1
59: *> Specifies whether condition numbers are required for
60: *> eigenvalues (S) or eigenvectors (SEP):
61: *> = 'E': for eigenvalues only (S);
62: *> = 'V': for eigenvectors only (SEP);
63: *> = 'B': for both eigenvalues and eigenvectors (S and SEP).
64: *> \endverbatim
65: *>
66: *> \param[in] HOWMNY
67: *> \verbatim
68: *> HOWMNY is CHARACTER*1
69: *> = 'A': compute condition numbers for all eigenpairs;
70: *> = 'S': compute condition numbers for selected eigenpairs
71: *> specified by the array SELECT.
72: *> \endverbatim
73: *>
74: *> \param[in] SELECT
75: *> \verbatim
76: *> SELECT is LOGICAL array, dimension (N)
77: *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
78: *> condition numbers are required. To select condition numbers
79: *> for the eigenpair corresponding to a real eigenvalue w(j),
80: *> SELECT(j) must be set to .TRUE.. To select condition numbers
81: *> corresponding to a complex conjugate pair of eigenvalues w(j)
82: *> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
83: *> set to .TRUE..
84: *> If HOWMNY = 'A', SELECT is not referenced.
85: *> \endverbatim
86: *>
87: *> \param[in] N
88: *> \verbatim
89: *> N is INTEGER
90: *> The order of the matrix T. N >= 0.
91: *> \endverbatim
92: *>
93: *> \param[in] T
94: *> \verbatim
95: *> T is DOUBLE PRECISION array, dimension (LDT,N)
96: *> The upper quasi-triangular matrix T, in Schur canonical form.
97: *> \endverbatim
98: *>
99: *> \param[in] LDT
100: *> \verbatim
101: *> LDT is INTEGER
102: *> The leading dimension of the array T. LDT >= max(1,N).
103: *> \endverbatim
104: *>
105: *> \param[in] VL
106: *> \verbatim
107: *> VL is DOUBLE PRECISION array, dimension (LDVL,M)
108: *> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
109: *> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
110: *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
111: *> must be stored in consecutive columns of VL, as returned by
112: *> DHSEIN or DTREVC.
113: *> If JOB = 'V', VL is not referenced.
114: *> \endverbatim
115: *>
116: *> \param[in] LDVL
117: *> \verbatim
118: *> LDVL is INTEGER
119: *> The leading dimension of the array VL.
120: *> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
121: *> \endverbatim
122: *>
123: *> \param[in] VR
124: *> \verbatim
125: *> VR is DOUBLE PRECISION array, dimension (LDVR,M)
126: *> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
127: *> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
128: *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
129: *> must be stored in consecutive columns of VR, as returned by
130: *> DHSEIN or DTREVC.
131: *> If JOB = 'V', VR is not referenced.
132: *> \endverbatim
133: *>
134: *> \param[in] LDVR
135: *> \verbatim
136: *> LDVR is INTEGER
137: *> The leading dimension of the array VR.
138: *> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
139: *> \endverbatim
140: *>
141: *> \param[out] S
142: *> \verbatim
143: *> S is DOUBLE PRECISION array, dimension (MM)
144: *> If JOB = 'E' or 'B', the reciprocal condition numbers of the
145: *> selected eigenvalues, stored in consecutive elements of the
146: *> array. For a complex conjugate pair of eigenvalues two
147: *> consecutive elements of S are set to the same value. Thus
148: *> S(j), SEP(j), and the j-th columns of VL and VR all
149: *> correspond to the same eigenpair (but not in general the
150: *> j-th eigenpair, unless all eigenpairs are selected).
151: *> If JOB = 'V', S is not referenced.
152: *> \endverbatim
153: *>
154: *> \param[out] SEP
155: *> \verbatim
156: *> SEP is DOUBLE PRECISION array, dimension (MM)
157: *> If JOB = 'V' or 'B', the estimated reciprocal condition
158: *> numbers of the selected eigenvectors, stored in consecutive
159: *> elements of the array. For a complex eigenvector two
160: *> consecutive elements of SEP are set to the same value. If
161: *> the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
162: *> is set to 0; this can only occur when the true value would be
163: *> very small anyway.
164: *> If JOB = 'E', SEP is not referenced.
165: *> \endverbatim
166: *>
167: *> \param[in] MM
168: *> \verbatim
169: *> MM is INTEGER
170: *> The number of elements in the arrays S (if JOB = 'E' or 'B')
171: *> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
172: *> \endverbatim
173: *>
174: *> \param[out] M
175: *> \verbatim
176: *> M is INTEGER
177: *> The number of elements of the arrays S and/or SEP actually
178: *> used to store the estimated condition numbers.
179: *> If HOWMNY = 'A', M is set to N.
180: *> \endverbatim
181: *>
182: *> \param[out] WORK
183: *> \verbatim
184: *> WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6)
185: *> If JOB = 'E', WORK is not referenced.
186: *> \endverbatim
187: *>
188: *> \param[in] LDWORK
189: *> \verbatim
190: *> LDWORK is INTEGER
191: *> The leading dimension of the array WORK.
192: *> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
193: *> \endverbatim
194: *>
195: *> \param[out] IWORK
196: *> \verbatim
197: *> IWORK is INTEGER array, dimension (2*(N-1))
198: *> If JOB = 'E', IWORK is not referenced.
199: *> \endverbatim
200: *>
201: *> \param[out] INFO
202: *> \verbatim
203: *> INFO is INTEGER
204: *> = 0: successful exit
205: *> < 0: if INFO = -i, the i-th argument had an illegal value
206: *> \endverbatim
207: *
208: * Authors:
209: * ========
210: *
211: *> \author Univ. of Tennessee
212: *> \author Univ. of California Berkeley
213: *> \author Univ. of Colorado Denver
214: *> \author NAG Ltd.
215: *
216: *> \ingroup doubleOTHERcomputational
217: *
218: *> \par Further Details:
219: * =====================
220: *>
221: *> \verbatim
222: *>
223: *> The reciprocal of the condition number of an eigenvalue lambda is
224: *> defined as
225: *>
226: *> S(lambda) = |v**T*u| / (norm(u)*norm(v))
227: *>
228: *> where u and v are the right and left eigenvectors of T corresponding
229: *> to lambda; v**T denotes the transpose of v, and norm(u)
230: *> denotes the Euclidean norm. These reciprocal condition numbers always
231: *> lie between zero (very badly conditioned) and one (very well
232: *> conditioned). If n = 1, S(lambda) is defined to be 1.
233: *>
234: *> An approximate error bound for a computed eigenvalue W(i) is given by
235: *>
236: *> EPS * norm(T) / S(i)
237: *>
238: *> where EPS is the machine precision.
239: *>
240: *> The reciprocal of the condition number of the right eigenvector u
241: *> corresponding to lambda is defined as follows. Suppose
242: *>
243: *> T = ( lambda c )
244: *> ( 0 T22 )
245: *>
246: *> Then the reciprocal condition number is
247: *>
248: *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
249: *>
250: *> where sigma-min denotes the smallest singular value. We approximate
251: *> the smallest singular value by the reciprocal of an estimate of the
252: *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
253: *> defined to be abs(T(1,1)).
254: *>
255: *> An approximate error bound for a computed right eigenvector VR(i)
256: *> is given by
257: *>
258: *> EPS * norm(T) / SEP(i)
259: *> \endverbatim
260: *>
261: * =====================================================================
262: SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
263: $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
264: $ INFO )
265: *
266: * -- LAPACK computational routine --
267: * -- LAPACK is a software package provided by Univ. of Tennessee, --
268: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269: *
270: * .. Scalar Arguments ..
271: CHARACTER HOWMNY, JOB
272: INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
273: * ..
274: * .. Array Arguments ..
275: LOGICAL SELECT( * )
276: INTEGER IWORK( * )
277: DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
278: $ VR( LDVR, * ), WORK( LDWORK, * )
279: * ..
280: *
281: * =====================================================================
282: *
283: * .. Parameters ..
284: DOUBLE PRECISION ZERO, ONE, TWO
285: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
286: * ..
287: * .. Local Scalars ..
288: LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP
289: INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
290: DOUBLE PRECISION BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
291: $ MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
292: * ..
293: * .. Local Arrays ..
294: INTEGER ISAVE( 3 )
295: DOUBLE PRECISION DUMMY( 1 )
296: * ..
297: * .. External Functions ..
298: LOGICAL LSAME
299: DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
300: EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
301: * ..
302: * .. External Subroutines ..
303: EXTERNAL DLABAD, DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
304: * ..
305: * .. Intrinsic Functions ..
306: INTRINSIC ABS, MAX, SQRT
307: * ..
308: * .. Executable Statements ..
309: *
310: * Decode and test the input parameters
311: *
312: WANTBH = LSAME( JOB, 'B' )
313: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
314: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
315: *
316: SOMCON = LSAME( HOWMNY, 'S' )
317: *
318: INFO = 0
319: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
320: INFO = -1
321: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
322: INFO = -2
323: ELSE IF( N.LT.0 ) THEN
324: INFO = -4
325: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
326: INFO = -6
327: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
328: INFO = -8
329: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
330: INFO = -10
331: ELSE
332: *
333: * Set M to the number of eigenpairs for which condition numbers
334: * are required, and test MM.
335: *
336: IF( SOMCON ) THEN
337: M = 0
338: PAIR = .FALSE.
339: DO 10 K = 1, N
340: IF( PAIR ) THEN
341: PAIR = .FALSE.
342: ELSE
343: IF( K.LT.N ) THEN
344: IF( T( K+1, K ).EQ.ZERO ) THEN
345: IF( SELECT( K ) )
346: $ M = M + 1
347: ELSE
348: PAIR = .TRUE.
349: IF( SELECT( K ) .OR. SELECT( K+1 ) )
350: $ M = M + 2
351: END IF
352: ELSE
353: IF( SELECT( N ) )
354: $ M = M + 1
355: END IF
356: END IF
357: 10 CONTINUE
358: ELSE
359: M = N
360: END IF
361: *
362: IF( MM.LT.M ) THEN
363: INFO = -13
364: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
365: INFO = -16
366: END IF
367: END IF
368: IF( INFO.NE.0 ) THEN
369: CALL XERBLA( 'DTRSNA', -INFO )
370: RETURN
371: END IF
372: *
373: * Quick return if possible
374: *
375: IF( N.EQ.0 )
376: $ RETURN
377: *
378: IF( N.EQ.1 ) THEN
379: IF( SOMCON ) THEN
380: IF( .NOT.SELECT( 1 ) )
381: $ RETURN
382: END IF
383: IF( WANTS )
384: $ S( 1 ) = ONE
385: IF( WANTSP )
386: $ SEP( 1 ) = ABS( T( 1, 1 ) )
387: RETURN
388: END IF
389: *
390: * Get machine constants
391: *
392: EPS = DLAMCH( 'P' )
393: SMLNUM = DLAMCH( 'S' ) / EPS
394: BIGNUM = ONE / SMLNUM
395: CALL DLABAD( SMLNUM, BIGNUM )
396: *
397: KS = 0
398: PAIR = .FALSE.
399: DO 60 K = 1, N
400: *
401: * Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
402: *
403: IF( PAIR ) THEN
404: PAIR = .FALSE.
405: GO TO 60
406: ELSE
407: IF( K.LT.N )
408: $ PAIR = T( K+1, K ).NE.ZERO
409: END IF
410: *
411: * Determine whether condition numbers are required for the k-th
412: * eigenpair.
413: *
414: IF( SOMCON ) THEN
415: IF( PAIR ) THEN
416: IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
417: $ GO TO 60
418: ELSE
419: IF( .NOT.SELECT( K ) )
420: $ GO TO 60
421: END IF
422: END IF
423: *
424: KS = KS + 1
425: *
426: IF( WANTS ) THEN
427: *
428: * Compute the reciprocal condition number of the k-th
429: * eigenvalue.
430: *
431: IF( .NOT.PAIR ) THEN
432: *
433: * Real eigenvalue.
434: *
435: PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
436: RNRM = DNRM2( N, VR( 1, KS ), 1 )
437: LNRM = DNRM2( N, VL( 1, KS ), 1 )
438: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
439: ELSE
440: *
441: * Complex eigenvalue.
442: *
443: PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
444: PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
445: $ 1 )
446: PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
447: PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
448: $ 1 )
449: RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
450: $ DNRM2( N, VR( 1, KS+1 ), 1 ) )
451: LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
452: $ DNRM2( N, VL( 1, KS+1 ), 1 ) )
453: COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
454: S( KS ) = COND
455: S( KS+1 ) = COND
456: END IF
457: END IF
458: *
459: IF( WANTSP ) THEN
460: *
461: * Estimate the reciprocal condition number of the k-th
462: * eigenvector.
463: *
464: * Copy the matrix T to the array WORK and swap the diagonal
465: * block beginning at T(k,k) to the (1,1) position.
466: *
467: CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
468: IFST = K
469: ILST = 1
470: CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
471: $ WORK( 1, N+1 ), IERR )
472: *
473: IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
474: *
475: * Could not swap because blocks not well separated
476: *
477: SCALE = ONE
478: EST = BIGNUM
479: ELSE
480: *
481: * Reordering successful
482: *
483: IF( WORK( 2, 1 ).EQ.ZERO ) THEN
484: *
485: * Form C = T22 - lambda*I in WORK(2:N,2:N).
486: *
487: DO 20 I = 2, N
488: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
489: 20 CONTINUE
490: N2 = 1
491: NN = N - 1
492: ELSE
493: *
494: * Triangularize the 2 by 2 block by unitary
495: * transformation U = [ cs i*ss ]
496: * [ i*ss cs ].
497: * such that the (1,1) position of WORK is complex
498: * eigenvalue lambda with positive imaginary part. (2,2)
499: * position of WORK is the complex eigenvalue lambda
500: * with negative imaginary part.
501: *
502: MU = SQRT( ABS( WORK( 1, 2 ) ) )*
503: $ SQRT( ABS( WORK( 2, 1 ) ) )
504: DELTA = DLAPY2( MU, WORK( 2, 1 ) )
505: CS = MU / DELTA
506: SN = -WORK( 2, 1 ) / DELTA
507: *
508: * Form
509: *
510: * C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
511: * [ mu ]
512: * [ .. ]
513: * [ .. ]
514: * [ mu ]
515: * where C**T is transpose of matrix C,
516: * and RWORK is stored starting in the N+1-st column of
517: * WORK.
518: *
519: DO 30 J = 3, N
520: WORK( 2, J ) = CS*WORK( 2, J )
521: WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
522: 30 CONTINUE
523: WORK( 2, 2 ) = ZERO
524: *
525: WORK( 1, N+1 ) = TWO*MU
526: DO 40 I = 2, N - 1
527: WORK( I, N+1 ) = SN*WORK( 1, I+1 )
528: 40 CONTINUE
529: N2 = 2
530: NN = 2*( N-1 )
531: END IF
532: *
533: * Estimate norm(inv(C**T))
534: *
535: EST = ZERO
536: KASE = 0
537: 50 CONTINUE
538: CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
539: $ EST, KASE, ISAVE )
540: IF( KASE.NE.0 ) THEN
541: IF( KASE.EQ.1 ) THEN
542: IF( N2.EQ.1 ) THEN
543: *
544: * Real eigenvalue: solve C**T*x = scale*c.
545: *
546: CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
547: $ LDWORK, DUMMY, DUMM, SCALE,
548: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
549: $ IERR )
550: ELSE
551: *
552: * Complex eigenvalue: solve
553: * C**T*(p+iq) = scale*(c+id) in real arithmetic.
554: *
555: CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
556: $ LDWORK, WORK( 1, N+1 ), MU, SCALE,
557: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
558: $ IERR )
559: END IF
560: ELSE
561: IF( N2.EQ.1 ) THEN
562: *
563: * Real eigenvalue: solve C*x = scale*c.
564: *
565: CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
566: $ LDWORK, DUMMY, DUMM, SCALE,
567: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
568: $ IERR )
569: ELSE
570: *
571: * Complex eigenvalue: solve
572: * C*(p+iq) = scale*(c+id) in real arithmetic.
573: *
574: CALL DLAQTR( .FALSE., .FALSE., N-1,
575: $ WORK( 2, 2 ), LDWORK,
576: $ WORK( 1, N+1 ), MU, SCALE,
577: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
578: $ IERR )
579: *
580: END IF
581: END IF
582: *
583: GO TO 50
584: END IF
585: END IF
586: *
587: SEP( KS ) = SCALE / MAX( EST, SMLNUM )
588: IF( PAIR )
589: $ SEP( KS+1 ) = SEP( KS )
590: END IF
591: *
592: IF( PAIR )
593: $ KS = KS + 1
594: *
595: 60 CONTINUE
596: RETURN
597: *
598: * End of DTRSNA
599: *
600: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrsna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack