1: SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
2: $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
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 DLACN2 in place of DLACON, 5 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: INTEGER IWORK( * )
19: DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
20: $ VR( LDVR, * ), WORK( LDWORK, * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * DTRSNA estimates reciprocal condition numbers for specified
27: * eigenvalues and/or right eigenvectors of a real upper
28: * quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
29: * orthogonal).
30: *
31: * T must be in Schur canonical form (as returned by DHSEQR), that is,
32: * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
33: * 2-by-2 diagonal block has its diagonal elements equal and its
34: * off-diagonal elements of opposite sign.
35: *
36: * Arguments
37: * =========
38: *
39: * JOB (input) CHARACTER*1
40: * Specifies whether condition numbers are required for
41: * eigenvalues (S) or eigenvectors (SEP):
42: * = 'E': for eigenvalues only (S);
43: * = 'V': for eigenvectors only (SEP);
44: * = 'B': for both eigenvalues and eigenvectors (S and SEP).
45: *
46: * HOWMNY (input) CHARACTER*1
47: * = 'A': compute condition numbers for all eigenpairs;
48: * = 'S': compute condition numbers for selected eigenpairs
49: * specified by the array SELECT.
50: *
51: * SELECT (input) LOGICAL array, dimension (N)
52: * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
53: * condition numbers are required. To select condition numbers
54: * for the eigenpair corresponding to a real eigenvalue w(j),
55: * SELECT(j) must be set to .TRUE.. To select condition numbers
56: * corresponding to a complex conjugate pair of eigenvalues w(j)
57: * and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
58: * set to .TRUE..
59: * If HOWMNY = 'A', SELECT is not referenced.
60: *
61: * N (input) INTEGER
62: * The order of the matrix T. N >= 0.
63: *
64: * T (input) DOUBLE PRECISION array, dimension (LDT,N)
65: * The upper quasi-triangular matrix T, in Schur canonical form.
66: *
67: * LDT (input) INTEGER
68: * The leading dimension of the array T. LDT >= max(1,N).
69: *
70: * VL (input) DOUBLE PRECISION array, dimension (LDVL,M)
71: * If JOB = 'E' or 'B', VL must contain left eigenvectors of T
72: * (or of any Q*T*Q**T with Q orthogonal), corresponding to the
73: * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
74: * must be stored in consecutive columns of VL, as returned by
75: * DHSEIN or DTREVC.
76: * If JOB = 'V', VL is not referenced.
77: *
78: * LDVL (input) INTEGER
79: * The leading dimension of the array VL.
80: * LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
81: *
82: * VR (input) DOUBLE PRECISION array, dimension (LDVR,M)
83: * If JOB = 'E' or 'B', VR must contain right eigenvectors of T
84: * (or of any Q*T*Q**T with Q orthogonal), corresponding to the
85: * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
86: * must be stored in consecutive columns of VR, as returned by
87: * DHSEIN or DTREVC.
88: * If JOB = 'V', VR is not referenced.
89: *
90: * LDVR (input) INTEGER
91: * The leading dimension of the array VR.
92: * LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
93: *
94: * S (output) DOUBLE PRECISION array, dimension (MM)
95: * If JOB = 'E' or 'B', the reciprocal condition numbers of the
96: * selected eigenvalues, stored in consecutive elements of the
97: * array. For a complex conjugate pair of eigenvalues two
98: * consecutive elements of S are set to the same value. Thus
99: * S(j), SEP(j), and the j-th columns of VL and VR all
100: * correspond to the same eigenpair (but not in general the
101: * j-th eigenpair, unless all eigenpairs are selected).
102: * If JOB = 'V', S is not referenced.
103: *
104: * SEP (output) DOUBLE PRECISION array, dimension (MM)
105: * If JOB = 'V' or 'B', the estimated reciprocal condition
106: * numbers of the selected eigenvectors, stored in consecutive
107: * elements of the array. For a complex eigenvector two
108: * consecutive elements of SEP are set to the same value. If
109: * the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
110: * is set to 0; this can only occur when the true value would be
111: * very small anyway.
112: * If JOB = 'E', SEP is not referenced.
113: *
114: * MM (input) INTEGER
115: * The number of elements in the arrays S (if JOB = 'E' or 'B')
116: * and/or SEP (if JOB = 'V' or 'B'). MM >= M.
117: *
118: * M (output) INTEGER
119: * The number of elements of the arrays S and/or SEP actually
120: * used to store the estimated condition numbers.
121: * If HOWMNY = 'A', M is set to N.
122: *
123: * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,N+6)
124: * If JOB = 'E', WORK is not referenced.
125: *
126: * LDWORK (input) INTEGER
127: * The leading dimension of the array WORK.
128: * LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
129: *
130: * IWORK (workspace) INTEGER array, dimension (2*(N-1))
131: * If JOB = 'E', IWORK is not referenced.
132: *
133: * INFO (output) INTEGER
134: * = 0: successful exit
135: * < 0: if INFO = -i, the i-th argument had an illegal value
136: *
137: * Further Details
138: * ===============
139: *
140: * The reciprocal of the condition number of an eigenvalue lambda is
141: * defined as
142: *
143: * S(lambda) = |v'*u| / (norm(u)*norm(v))
144: *
145: * where u and v are the right and left eigenvectors of T corresponding
146: * to lambda; v' denotes the conjugate-transpose of v, and norm(u)
147: * denotes the Euclidean norm. These reciprocal condition numbers always
148: * lie between zero (very badly conditioned) and one (very well
149: * conditioned). If n = 1, S(lambda) is defined to be 1.
150: *
151: * An approximate error bound for a computed eigenvalue W(i) is given by
152: *
153: * EPS * norm(T) / S(i)
154: *
155: * where EPS is the machine precision.
156: *
157: * The reciprocal of the condition number of the right eigenvector u
158: * corresponding to lambda is defined as follows. Suppose
159: *
160: * T = ( lambda c )
161: * ( 0 T22 )
162: *
163: * Then the reciprocal condition number is
164: *
165: * SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
166: *
167: * where sigma-min denotes the smallest singular value. We approximate
168: * the smallest singular value by the reciprocal of an estimate of the
169: * one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
170: * defined to be abs(T(1,1)).
171: *
172: * An approximate error bound for a computed right eigenvector VR(i)
173: * is given by
174: *
175: * EPS * norm(T) / SEP(i)
176: *
177: * =====================================================================
178: *
179: * .. Parameters ..
180: DOUBLE PRECISION ZERO, ONE, TWO
181: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
182: * ..
183: * .. Local Scalars ..
184: LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP
185: INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
186: DOUBLE PRECISION BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
187: $ MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
188: * ..
189: * .. Local Arrays ..
190: INTEGER ISAVE( 3 )
191: DOUBLE PRECISION DUMMY( 1 )
192: * ..
193: * .. External Functions ..
194: LOGICAL LSAME
195: DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
196: EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
197: * ..
198: * .. External Subroutines ..
199: EXTERNAL DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
200: * ..
201: * .. Intrinsic Functions ..
202: INTRINSIC ABS, MAX, SQRT
203: * ..
204: * .. Executable Statements ..
205: *
206: * Decode and test the input parameters
207: *
208: WANTBH = LSAME( JOB, 'B' )
209: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
210: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
211: *
212: SOMCON = LSAME( HOWMNY, 'S' )
213: *
214: INFO = 0
215: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
216: INFO = -1
217: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
218: INFO = -2
219: ELSE IF( N.LT.0 ) THEN
220: INFO = -4
221: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
222: INFO = -6
223: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
224: INFO = -8
225: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
226: INFO = -10
227: ELSE
228: *
229: * Set M to the number of eigenpairs for which condition numbers
230: * are required, and test MM.
231: *
232: IF( SOMCON ) THEN
233: M = 0
234: PAIR = .FALSE.
235: DO 10 K = 1, N
236: IF( PAIR ) THEN
237: PAIR = .FALSE.
238: ELSE
239: IF( K.LT.N ) THEN
240: IF( T( K+1, K ).EQ.ZERO ) THEN
241: IF( SELECT( K ) )
242: $ M = M + 1
243: ELSE
244: PAIR = .TRUE.
245: IF( SELECT( K ) .OR. SELECT( K+1 ) )
246: $ M = M + 2
247: END IF
248: ELSE
249: IF( SELECT( N ) )
250: $ M = M + 1
251: END IF
252: END IF
253: 10 CONTINUE
254: ELSE
255: M = N
256: END IF
257: *
258: IF( MM.LT.M ) THEN
259: INFO = -13
260: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
261: INFO = -16
262: END IF
263: END IF
264: IF( INFO.NE.0 ) THEN
265: CALL XERBLA( 'DTRSNA', -INFO )
266: RETURN
267: END IF
268: *
269: * Quick return if possible
270: *
271: IF( N.EQ.0 )
272: $ RETURN
273: *
274: IF( N.EQ.1 ) THEN
275: IF( SOMCON ) THEN
276: IF( .NOT.SELECT( 1 ) )
277: $ RETURN
278: END IF
279: IF( WANTS )
280: $ S( 1 ) = ONE
281: IF( WANTSP )
282: $ SEP( 1 ) = ABS( T( 1, 1 ) )
283: RETURN
284: END IF
285: *
286: * Get machine constants
287: *
288: EPS = DLAMCH( 'P' )
289: SMLNUM = DLAMCH( 'S' ) / EPS
290: BIGNUM = ONE / SMLNUM
291: CALL DLABAD( SMLNUM, BIGNUM )
292: *
293: KS = 0
294: PAIR = .FALSE.
295: DO 60 K = 1, N
296: *
297: * Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
298: *
299: IF( PAIR ) THEN
300: PAIR = .FALSE.
301: GO TO 60
302: ELSE
303: IF( K.LT.N )
304: $ PAIR = T( K+1, K ).NE.ZERO
305: END IF
306: *
307: * Determine whether condition numbers are required for the k-th
308: * eigenpair.
309: *
310: IF( SOMCON ) THEN
311: IF( PAIR ) THEN
312: IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
313: $ GO TO 60
314: ELSE
315: IF( .NOT.SELECT( K ) )
316: $ GO TO 60
317: END IF
318: END IF
319: *
320: KS = KS + 1
321: *
322: IF( WANTS ) THEN
323: *
324: * Compute the reciprocal condition number of the k-th
325: * eigenvalue.
326: *
327: IF( .NOT.PAIR ) THEN
328: *
329: * Real eigenvalue.
330: *
331: PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
332: RNRM = DNRM2( N, VR( 1, KS ), 1 )
333: LNRM = DNRM2( N, VL( 1, KS ), 1 )
334: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
335: ELSE
336: *
337: * Complex eigenvalue.
338: *
339: PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
340: PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
341: $ 1 )
342: PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
343: PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
344: $ 1 )
345: RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
346: $ DNRM2( N, VR( 1, KS+1 ), 1 ) )
347: LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
348: $ DNRM2( N, VL( 1, KS+1 ), 1 ) )
349: COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
350: S( KS ) = COND
351: S( KS+1 ) = COND
352: END IF
353: END IF
354: *
355: IF( WANTSP ) THEN
356: *
357: * Estimate the reciprocal condition number of the k-th
358: * eigenvector.
359: *
360: * Copy the matrix T to the array WORK and swap the diagonal
361: * block beginning at T(k,k) to the (1,1) position.
362: *
363: CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
364: IFST = K
365: ILST = 1
366: CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
367: $ WORK( 1, N+1 ), IERR )
368: *
369: IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
370: *
371: * Could not swap because blocks not well separated
372: *
373: SCALE = ONE
374: EST = BIGNUM
375: ELSE
376: *
377: * Reordering successful
378: *
379: IF( WORK( 2, 1 ).EQ.ZERO ) THEN
380: *
381: * Form C = T22 - lambda*I in WORK(2:N,2:N).
382: *
383: DO 20 I = 2, N
384: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
385: 20 CONTINUE
386: N2 = 1
387: NN = N - 1
388: ELSE
389: *
390: * Triangularize the 2 by 2 block by unitary
391: * transformation U = [ cs i*ss ]
392: * [ i*ss cs ].
393: * such that the (1,1) position of WORK is complex
394: * eigenvalue lambda with positive imaginary part. (2,2)
395: * position of WORK is the complex eigenvalue lambda
396: * with negative imaginary part.
397: *
398: MU = SQRT( ABS( WORK( 1, 2 ) ) )*
399: $ SQRT( ABS( WORK( 2, 1 ) ) )
400: DELTA = DLAPY2( MU, WORK( 2, 1 ) )
401: CS = MU / DELTA
402: SN = -WORK( 2, 1 ) / DELTA
403: *
404: * Form
405: *
406: * C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
407: * [ mu ]
408: * [ .. ]
409: * [ .. ]
410: * [ mu ]
411: * where C' is conjugate transpose of complex matrix C,
412: * and RWORK is stored starting in the N+1-st column of
413: * WORK.
414: *
415: DO 30 J = 3, N
416: WORK( 2, J ) = CS*WORK( 2, J )
417: WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
418: 30 CONTINUE
419: WORK( 2, 2 ) = ZERO
420: *
421: WORK( 1, N+1 ) = TWO*MU
422: DO 40 I = 2, N - 1
423: WORK( I, N+1 ) = SN*WORK( 1, I+1 )
424: 40 CONTINUE
425: N2 = 2
426: NN = 2*( N-1 )
427: END IF
428: *
429: * Estimate norm(inv(C'))
430: *
431: EST = ZERO
432: KASE = 0
433: 50 CONTINUE
434: CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
435: $ EST, KASE, ISAVE )
436: IF( KASE.NE.0 ) THEN
437: IF( KASE.EQ.1 ) THEN
438: IF( N2.EQ.1 ) THEN
439: *
440: * Real eigenvalue: solve C'*x = scale*c.
441: *
442: CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
443: $ LDWORK, DUMMY, DUMM, SCALE,
444: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
445: $ IERR )
446: ELSE
447: *
448: * Complex eigenvalue: solve
449: * C'*(p+iq) = scale*(c+id) in real arithmetic.
450: *
451: CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
452: $ LDWORK, WORK( 1, N+1 ), MU, SCALE,
453: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
454: $ IERR )
455: END IF
456: ELSE
457: IF( N2.EQ.1 ) THEN
458: *
459: * Real eigenvalue: solve C*x = scale*c.
460: *
461: CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
462: $ LDWORK, DUMMY, DUMM, SCALE,
463: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
464: $ IERR )
465: ELSE
466: *
467: * Complex eigenvalue: solve
468: * C*(p+iq) = scale*(c+id) in real arithmetic.
469: *
470: CALL DLAQTR( .FALSE., .FALSE., N-1,
471: $ WORK( 2, 2 ), LDWORK,
472: $ WORK( 1, N+1 ), MU, SCALE,
473: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
474: $ IERR )
475: *
476: END IF
477: END IF
478: *
479: GO TO 50
480: END IF
481: END IF
482: *
483: SEP( KS ) = SCALE / MAX( EST, SMLNUM )
484: IF( PAIR )
485: $ SEP( KS+1 ) = SEP( KS )
486: END IF
487: *
488: IF( PAIR )
489: $ KS = KS + 1
490: *
491: 60 CONTINUE
492: RETURN
493: *
494: * End of DTRSNA
495: *
496: 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