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.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
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**T*u| / (norm(u)*norm(v))
144 *
145 * where u and v are the right and left eigenvectors of T corresponding
146 * to lambda; v**T denotes the 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**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
407 * [ mu ]
408 * [ .. ]
409 * [ .. ]
410 * [ mu ]
411 * where C**T is transpose of 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**T))
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**T*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**T*(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
2 $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
3 $ INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
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**T*u| / (norm(u)*norm(v))
144 *
145 * where u and v are the right and left eigenvectors of T corresponding
146 * to lambda; v**T denotes the 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**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
407 * [ mu ]
408 * [ .. ]
409 * [ .. ]
410 * [ mu ]
411 * where C**T is transpose of 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**T))
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**T*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**T*(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