1 SUBROUTINE DGET38( RMAX, LMAX, NINFO, KNT, NIN )
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER KNT, NIN
9 * ..
10 * .. Array Arguments ..
11 INTEGER LMAX( 3 ), NINFO( 3 )
12 DOUBLE PRECISION RMAX( 3 )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DGET38 tests DTRSEN, a routine for estimating condition numbers of a
19 * cluster of eigenvalues and/or its associated right invariant subspace
20 *
21 * The test matrices are read from a file with logical unit number NIN.
22 *
23 * Arguments
24 * ==========
25 *
26 * RMAX (output) DOUBLE PRECISION array, dimension (3)
27 * Values of the largest test ratios.
28 * RMAX(1) = largest residuals from DHST01 or comparing
29 * different calls to DTRSEN
30 * RMAX(2) = largest error in reciprocal condition
31 * numbers taking their conditioning into account
32 * RMAX(3) = largest error in reciprocal condition
33 * numbers not taking their conditioning into
34 * account (may be larger than RMAX(2))
35 *
36 * LMAX (output) INTEGER array, dimension (3)
37 * LMAX(i) is example number where largest test ratio
38 * RMAX(i) is achieved. Also:
39 * If DGEHRD returns INFO nonzero on example i, LMAX(1)=i
40 * If DHSEQR returns INFO nonzero on example i, LMAX(2)=i
41 * If DTRSEN returns INFO nonzero on example i, LMAX(3)=i
42 *
43 * NINFO (output) INTEGER array, dimension (3)
44 * NINFO(1) = No. of times DGEHRD returned INFO nonzero
45 * NINFO(2) = No. of times DHSEQR returned INFO nonzero
46 * NINFO(3) = No. of times DTRSEN returned INFO nonzero
47 *
48 * KNT (output) INTEGER
49 * Total number of examples tested.
50 *
51 * NIN (input) INTEGER
52 * Input logical unit number.
53 *
54 * =====================================================================
55 *
56 * .. Parameters ..
57 DOUBLE PRECISION ZERO, ONE, TWO
58 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
59 DOUBLE PRECISION EPSIN
60 PARAMETER ( EPSIN = 5.9605D-8 )
61 INTEGER LDT, LWORK
62 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
63 INTEGER LIWORK
64 PARAMETER ( LIWORK = LDT*LDT )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, INFO, ISCL, ITMP, J, KMIN, M, N, NDIM
68 DOUBLE PRECISION BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
69 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VIMIN, VMAX,
70 $ VMUL, VRMIN
71 * ..
72 * .. Local Arrays ..
73 LOGICAL SELECT( LDT )
74 INTEGER IPNT( LDT ), ISELEC( LDT ), IWORK( LIWORK )
75 DOUBLE PRECISION Q( LDT, LDT ), QSAV( LDT, LDT ),
76 $ QTMP( LDT, LDT ), RESULT( 2 ), T( LDT, LDT ),
77 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
78 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), VAL( 3 ),
79 $ WI( LDT ), WITMP( LDT ), WORK( LWORK ),
80 $ WR( LDT ), WRTMP( LDT )
81 * ..
82 * .. External Functions ..
83 DOUBLE PRECISION DLAMCH, DLANGE
84 EXTERNAL DLAMCH, DLANGE
85 * ..
86 * .. External Subroutines ..
87 EXTERNAL DCOPY, DGEHRD, DHSEQR, DHST01, DLABAD, DLACPY,
88 $ DORGHR, DSCAL, DTRSEN
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC DBLE, MAX, SQRT
92 * ..
93 * .. Executable Statements ..
94 *
95 EPS = DLAMCH( 'P' )
96 SMLNUM = DLAMCH( 'S' ) / EPS
97 BIGNUM = ONE / SMLNUM
98 CALL DLABAD( SMLNUM, BIGNUM )
99 *
100 * EPSIN = 2**(-24) = precision to which input data computed
101 *
102 EPS = MAX( EPS, EPSIN )
103 RMAX( 1 ) = ZERO
104 RMAX( 2 ) = ZERO
105 RMAX( 3 ) = ZERO
106 LMAX( 1 ) = 0
107 LMAX( 2 ) = 0
108 LMAX( 3 ) = 0
109 KNT = 0
110 NINFO( 1 ) = 0
111 NINFO( 2 ) = 0
112 NINFO( 3 ) = 0
113 *
114 VAL( 1 ) = SQRT( SMLNUM )
115 VAL( 2 ) = ONE
116 VAL( 3 ) = SQRT( SQRT( BIGNUM ) )
117 *
118 * Read input data until N=0. Assume input eigenvalues are sorted
119 * lexicographically (increasing by real part, then decreasing by
120 * imaginary part)
121 *
122 10 CONTINUE
123 READ( NIN, FMT = * )N, NDIM
124 IF( N.EQ.0 )
125 $ RETURN
126 READ( NIN, FMT = * )( ISELEC( I ), I = 1, NDIM )
127 DO 20 I = 1, N
128 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
129 20 CONTINUE
130 READ( NIN, FMT = * )SIN, SEPIN
131 *
132 TNRM = DLANGE( 'M', N, N, TMP, LDT, WORK )
133 DO 160 ISCL = 1, 3
134 *
135 * Scale input matrix
136 *
137 KNT = KNT + 1
138 CALL DLACPY( 'F', N, N, TMP, LDT, T, LDT )
139 VMUL = VAL( ISCL )
140 DO 30 I = 1, N
141 CALL DSCAL( N, VMUL, T( 1, I ), 1 )
142 30 CONTINUE
143 IF( TNRM.EQ.ZERO )
144 $ VMUL = ONE
145 CALL DLACPY( 'F', N, N, T, LDT, TSAV, LDT )
146 *
147 * Compute Schur form
148 *
149 CALL DGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
150 $ INFO )
151 IF( INFO.NE.0 ) THEN
152 LMAX( 1 ) = KNT
153 NINFO( 1 ) = NINFO( 1 ) + 1
154 GO TO 160
155 END IF
156 *
157 * Generate orthogonal matrix
158 *
159 CALL DLACPY( 'L', N, N, T, LDT, Q, LDT )
160 CALL DORGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
161 $ INFO )
162 *
163 * Compute Schur form
164 *
165 CALL DHSEQR( 'S', 'V', N, 1, N, T, LDT, WR, WI, Q, LDT, WORK,
166 $ LWORK, INFO )
167 IF( INFO.NE.0 ) THEN
168 LMAX( 2 ) = KNT
169 NINFO( 2 ) = NINFO( 2 ) + 1
170 GO TO 160
171 END IF
172 *
173 * Sort, select eigenvalues
174 *
175 DO 40 I = 1, N
176 IPNT( I ) = I
177 SELECT( I ) = .FALSE.
178 40 CONTINUE
179 CALL DCOPY( N, WR, 1, WRTMP, 1 )
180 CALL DCOPY( N, WI, 1, WITMP, 1 )
181 DO 60 I = 1, N - 1
182 KMIN = I
183 VRMIN = WRTMP( I )
184 VIMIN = WITMP( I )
185 DO 50 J = I + 1, N
186 IF( WRTMP( J ).LT.VRMIN ) THEN
187 KMIN = J
188 VRMIN = WRTMP( J )
189 VIMIN = WITMP( J )
190 END IF
191 50 CONTINUE
192 WRTMP( KMIN ) = WRTMP( I )
193 WITMP( KMIN ) = WITMP( I )
194 WRTMP( I ) = VRMIN
195 WITMP( I ) = VIMIN
196 ITMP = IPNT( I )
197 IPNT( I ) = IPNT( KMIN )
198 IPNT( KMIN ) = ITMP
199 60 CONTINUE
200 DO 70 I = 1, NDIM
201 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE.
202 70 CONTINUE
203 *
204 * Compute condition numbers
205 *
206 CALL DLACPY( 'F', N, N, Q, LDT, QSAV, LDT )
207 CALL DLACPY( 'F', N, N, T, LDT, TSAV1, LDT )
208 CALL DTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WRTMP, WITMP,
209 $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
210 IF( INFO.NE.0 ) THEN
211 LMAX( 3 ) = KNT
212 NINFO( 3 ) = NINFO( 3 ) + 1
213 GO TO 160
214 END IF
215 SEPTMP = SEP / VMUL
216 STMP = S
217 *
218 * Compute residuals
219 *
220 CALL DHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK,
221 $ RESULT )
222 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) )
223 IF( VMAX.GT.RMAX( 1 ) ) THEN
224 RMAX( 1 ) = VMAX
225 IF( NINFO( 1 ).EQ.0 )
226 $ LMAX( 1 ) = KNT
227 END IF
228 *
229 * Compare condition number for eigenvalue cluster
230 * taking its condition number into account
231 *
232 V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM )
233 IF( TNRM.EQ.ZERO )
234 $ V = ONE
235 IF( V.GT.SEPTMP ) THEN
236 TOL = ONE
237 ELSE
238 TOL = V / SEPTMP
239 END IF
240 IF( V.GT.SEPIN ) THEN
241 TOLIN = ONE
242 ELSE
243 TOLIN = V / SEPIN
244 END IF
245 TOL = MAX( TOL, SMLNUM / EPS )
246 TOLIN = MAX( TOLIN, SMLNUM / EPS )
247 IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN
248 VMAX = ONE / EPS
249 ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN
250 VMAX = ( SIN-TOLIN ) / ( STMP+TOL )
251 ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN
252 VMAX = ONE / EPS
253 ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN
254 VMAX = ( STMP-TOL ) / ( SIN+TOLIN )
255 ELSE
256 VMAX = ONE
257 END IF
258 IF( VMAX.GT.RMAX( 2 ) ) THEN
259 RMAX( 2 ) = VMAX
260 IF( NINFO( 2 ).EQ.0 )
261 $ LMAX( 2 ) = KNT
262 END IF
263 *
264 * Compare condition numbers for invariant subspace
265 * taking its condition number into account
266 *
267 IF( V.GT.SEPTMP*STMP ) THEN
268 TOL = SEPTMP
269 ELSE
270 TOL = V / STMP
271 END IF
272 IF( V.GT.SEPIN*SIN ) THEN
273 TOLIN = SEPIN
274 ELSE
275 TOLIN = V / SIN
276 END IF
277 TOL = MAX( TOL, SMLNUM / EPS )
278 TOLIN = MAX( TOLIN, SMLNUM / EPS )
279 IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN
280 VMAX = ONE / EPS
281 ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN
282 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL )
283 ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN
284 VMAX = ONE / EPS
285 ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN
286 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN )
287 ELSE
288 VMAX = ONE
289 END IF
290 IF( VMAX.GT.RMAX( 2 ) ) THEN
291 RMAX( 2 ) = VMAX
292 IF( NINFO( 2 ).EQ.0 )
293 $ LMAX( 2 ) = KNT
294 END IF
295 *
296 * Compare condition number for eigenvalue cluster
297 * without taking its condition number into account
298 *
299 IF( SIN.LE.DBLE( 2*N )*EPS .AND. STMP.LE.DBLE( 2*N )*EPS ) THEN
300 VMAX = ONE
301 ELSE IF( EPS*SIN.GT.STMP ) THEN
302 VMAX = ONE / EPS
303 ELSE IF( SIN.GT.STMP ) THEN
304 VMAX = SIN / STMP
305 ELSE IF( SIN.LT.EPS*STMP ) THEN
306 VMAX = ONE / EPS
307 ELSE IF( SIN.LT.STMP ) THEN
308 VMAX = STMP / SIN
309 ELSE
310 VMAX = ONE
311 END IF
312 IF( VMAX.GT.RMAX( 3 ) ) THEN
313 RMAX( 3 ) = VMAX
314 IF( NINFO( 3 ).EQ.0 )
315 $ LMAX( 3 ) = KNT
316 END IF
317 *
318 * Compare condition numbers for invariant subspace
319 * without taking its condition number into account
320 *
321 IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN
322 VMAX = ONE
323 ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN
324 VMAX = ONE / EPS
325 ELSE IF( SEPIN.GT.SEPTMP ) THEN
326 VMAX = SEPIN / SEPTMP
327 ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN
328 VMAX = ONE / EPS
329 ELSE IF( SEPIN.LT.SEPTMP ) THEN
330 VMAX = SEPTMP / SEPIN
331 ELSE
332 VMAX = ONE
333 END IF
334 IF( VMAX.GT.RMAX( 3 ) ) THEN
335 RMAX( 3 ) = VMAX
336 IF( NINFO( 3 ).EQ.0 )
337 $ LMAX( 3 ) = KNT
338 END IF
339 *
340 * Compute eigenvalue condition number only and compare
341 * Update Q
342 *
343 VMAX = ZERO
344 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
345 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
346 SEPTMP = -ONE
347 STMP = -ONE
348 CALL DTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
349 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
350 $ LIWORK, INFO )
351 IF( INFO.NE.0 ) THEN
352 LMAX( 3 ) = KNT
353 NINFO( 3 ) = NINFO( 3 ) + 1
354 GO TO 160
355 END IF
356 IF( S.NE.STMP )
357 $ VMAX = ONE / EPS
358 IF( -ONE.NE.SEPTMP )
359 $ VMAX = ONE / EPS
360 DO 90 I = 1, N
361 DO 80 J = 1, N
362 IF( TTMP( I, J ).NE.T( I, J ) )
363 $ VMAX = ONE / EPS
364 IF( QTMP( I, J ).NE.Q( I, J ) )
365 $ VMAX = ONE / EPS
366 80 CONTINUE
367 90 CONTINUE
368 *
369 * Compute invariant subspace condition number only and compare
370 * Update Q
371 *
372 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
373 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
374 SEPTMP = -ONE
375 STMP = -ONE
376 CALL DTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
377 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
378 $ LIWORK, INFO )
379 IF( INFO.NE.0 ) THEN
380 LMAX( 3 ) = KNT
381 NINFO( 3 ) = NINFO( 3 ) + 1
382 GO TO 160
383 END IF
384 IF( -ONE.NE.STMP )
385 $ VMAX = ONE / EPS
386 IF( SEP.NE.SEPTMP )
387 $ VMAX = ONE / EPS
388 DO 110 I = 1, N
389 DO 100 J = 1, N
390 IF( TTMP( I, J ).NE.T( I, J ) )
391 $ VMAX = ONE / EPS
392 IF( QTMP( I, J ).NE.Q( I, J ) )
393 $ VMAX = ONE / EPS
394 100 CONTINUE
395 110 CONTINUE
396 *
397 * Compute eigenvalue condition number only and compare
398 * Do not update Q
399 *
400 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
401 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
402 SEPTMP = -ONE
403 STMP = -ONE
404 CALL DTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
405 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
406 $ LIWORK, INFO )
407 IF( INFO.NE.0 ) THEN
408 LMAX( 3 ) = KNT
409 NINFO( 3 ) = NINFO( 3 ) + 1
410 GO TO 160
411 END IF
412 IF( S.NE.STMP )
413 $ VMAX = ONE / EPS
414 IF( -ONE.NE.SEPTMP )
415 $ VMAX = ONE / EPS
416 DO 130 I = 1, N
417 DO 120 J = 1, N
418 IF( TTMP( I, J ).NE.T( I, J ) )
419 $ VMAX = ONE / EPS
420 IF( QTMP( I, J ).NE.QSAV( I, J ) )
421 $ VMAX = ONE / EPS
422 120 CONTINUE
423 130 CONTINUE
424 *
425 * Compute invariant subspace condition number only and compare
426 * Do not update Q
427 *
428 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
429 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
430 SEPTMP = -ONE
431 STMP = -ONE
432 CALL DTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
433 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
434 $ LIWORK, INFO )
435 IF( INFO.NE.0 ) THEN
436 LMAX( 3 ) = KNT
437 NINFO( 3 ) = NINFO( 3 ) + 1
438 GO TO 160
439 END IF
440 IF( -ONE.NE.STMP )
441 $ VMAX = ONE / EPS
442 IF( SEP.NE.SEPTMP )
443 $ VMAX = ONE / EPS
444 DO 150 I = 1, N
445 DO 140 J = 1, N
446 IF( TTMP( I, J ).NE.T( I, J ) )
447 $ VMAX = ONE / EPS
448 IF( QTMP( I, J ).NE.QSAV( I, J ) )
449 $ VMAX = ONE / EPS
450 140 CONTINUE
451 150 CONTINUE
452 IF( VMAX.GT.RMAX( 1 ) ) THEN
453 RMAX( 1 ) = VMAX
454 IF( NINFO( 1 ).EQ.0 )
455 $ LMAX( 1 ) = KNT
456 END IF
457 160 CONTINUE
458 GO TO 10
459 *
460 * End of DGET38
461 *
462 END
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER KNT, NIN
9 * ..
10 * .. Array Arguments ..
11 INTEGER LMAX( 3 ), NINFO( 3 )
12 DOUBLE PRECISION RMAX( 3 )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DGET38 tests DTRSEN, a routine for estimating condition numbers of a
19 * cluster of eigenvalues and/or its associated right invariant subspace
20 *
21 * The test matrices are read from a file with logical unit number NIN.
22 *
23 * Arguments
24 * ==========
25 *
26 * RMAX (output) DOUBLE PRECISION array, dimension (3)
27 * Values of the largest test ratios.
28 * RMAX(1) = largest residuals from DHST01 or comparing
29 * different calls to DTRSEN
30 * RMAX(2) = largest error in reciprocal condition
31 * numbers taking their conditioning into account
32 * RMAX(3) = largest error in reciprocal condition
33 * numbers not taking their conditioning into
34 * account (may be larger than RMAX(2))
35 *
36 * LMAX (output) INTEGER array, dimension (3)
37 * LMAX(i) is example number where largest test ratio
38 * RMAX(i) is achieved. Also:
39 * If DGEHRD returns INFO nonzero on example i, LMAX(1)=i
40 * If DHSEQR returns INFO nonzero on example i, LMAX(2)=i
41 * If DTRSEN returns INFO nonzero on example i, LMAX(3)=i
42 *
43 * NINFO (output) INTEGER array, dimension (3)
44 * NINFO(1) = No. of times DGEHRD returned INFO nonzero
45 * NINFO(2) = No. of times DHSEQR returned INFO nonzero
46 * NINFO(3) = No. of times DTRSEN returned INFO nonzero
47 *
48 * KNT (output) INTEGER
49 * Total number of examples tested.
50 *
51 * NIN (input) INTEGER
52 * Input logical unit number.
53 *
54 * =====================================================================
55 *
56 * .. Parameters ..
57 DOUBLE PRECISION ZERO, ONE, TWO
58 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
59 DOUBLE PRECISION EPSIN
60 PARAMETER ( EPSIN = 5.9605D-8 )
61 INTEGER LDT, LWORK
62 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
63 INTEGER LIWORK
64 PARAMETER ( LIWORK = LDT*LDT )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, INFO, ISCL, ITMP, J, KMIN, M, N, NDIM
68 DOUBLE PRECISION BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
69 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VIMIN, VMAX,
70 $ VMUL, VRMIN
71 * ..
72 * .. Local Arrays ..
73 LOGICAL SELECT( LDT )
74 INTEGER IPNT( LDT ), ISELEC( LDT ), IWORK( LIWORK )
75 DOUBLE PRECISION Q( LDT, LDT ), QSAV( LDT, LDT ),
76 $ QTMP( LDT, LDT ), RESULT( 2 ), T( LDT, LDT ),
77 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
78 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), VAL( 3 ),
79 $ WI( LDT ), WITMP( LDT ), WORK( LWORK ),
80 $ WR( LDT ), WRTMP( LDT )
81 * ..
82 * .. External Functions ..
83 DOUBLE PRECISION DLAMCH, DLANGE
84 EXTERNAL DLAMCH, DLANGE
85 * ..
86 * .. External Subroutines ..
87 EXTERNAL DCOPY, DGEHRD, DHSEQR, DHST01, DLABAD, DLACPY,
88 $ DORGHR, DSCAL, DTRSEN
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC DBLE, MAX, SQRT
92 * ..
93 * .. Executable Statements ..
94 *
95 EPS = DLAMCH( 'P' )
96 SMLNUM = DLAMCH( 'S' ) / EPS
97 BIGNUM = ONE / SMLNUM
98 CALL DLABAD( SMLNUM, BIGNUM )
99 *
100 * EPSIN = 2**(-24) = precision to which input data computed
101 *
102 EPS = MAX( EPS, EPSIN )
103 RMAX( 1 ) = ZERO
104 RMAX( 2 ) = ZERO
105 RMAX( 3 ) = ZERO
106 LMAX( 1 ) = 0
107 LMAX( 2 ) = 0
108 LMAX( 3 ) = 0
109 KNT = 0
110 NINFO( 1 ) = 0
111 NINFO( 2 ) = 0
112 NINFO( 3 ) = 0
113 *
114 VAL( 1 ) = SQRT( SMLNUM )
115 VAL( 2 ) = ONE
116 VAL( 3 ) = SQRT( SQRT( BIGNUM ) )
117 *
118 * Read input data until N=0. Assume input eigenvalues are sorted
119 * lexicographically (increasing by real part, then decreasing by
120 * imaginary part)
121 *
122 10 CONTINUE
123 READ( NIN, FMT = * )N, NDIM
124 IF( N.EQ.0 )
125 $ RETURN
126 READ( NIN, FMT = * )( ISELEC( I ), I = 1, NDIM )
127 DO 20 I = 1, N
128 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
129 20 CONTINUE
130 READ( NIN, FMT = * )SIN, SEPIN
131 *
132 TNRM = DLANGE( 'M', N, N, TMP, LDT, WORK )
133 DO 160 ISCL = 1, 3
134 *
135 * Scale input matrix
136 *
137 KNT = KNT + 1
138 CALL DLACPY( 'F', N, N, TMP, LDT, T, LDT )
139 VMUL = VAL( ISCL )
140 DO 30 I = 1, N
141 CALL DSCAL( N, VMUL, T( 1, I ), 1 )
142 30 CONTINUE
143 IF( TNRM.EQ.ZERO )
144 $ VMUL = ONE
145 CALL DLACPY( 'F', N, N, T, LDT, TSAV, LDT )
146 *
147 * Compute Schur form
148 *
149 CALL DGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
150 $ INFO )
151 IF( INFO.NE.0 ) THEN
152 LMAX( 1 ) = KNT
153 NINFO( 1 ) = NINFO( 1 ) + 1
154 GO TO 160
155 END IF
156 *
157 * Generate orthogonal matrix
158 *
159 CALL DLACPY( 'L', N, N, T, LDT, Q, LDT )
160 CALL DORGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
161 $ INFO )
162 *
163 * Compute Schur form
164 *
165 CALL DHSEQR( 'S', 'V', N, 1, N, T, LDT, WR, WI, Q, LDT, WORK,
166 $ LWORK, INFO )
167 IF( INFO.NE.0 ) THEN
168 LMAX( 2 ) = KNT
169 NINFO( 2 ) = NINFO( 2 ) + 1
170 GO TO 160
171 END IF
172 *
173 * Sort, select eigenvalues
174 *
175 DO 40 I = 1, N
176 IPNT( I ) = I
177 SELECT( I ) = .FALSE.
178 40 CONTINUE
179 CALL DCOPY( N, WR, 1, WRTMP, 1 )
180 CALL DCOPY( N, WI, 1, WITMP, 1 )
181 DO 60 I = 1, N - 1
182 KMIN = I
183 VRMIN = WRTMP( I )
184 VIMIN = WITMP( I )
185 DO 50 J = I + 1, N
186 IF( WRTMP( J ).LT.VRMIN ) THEN
187 KMIN = J
188 VRMIN = WRTMP( J )
189 VIMIN = WITMP( J )
190 END IF
191 50 CONTINUE
192 WRTMP( KMIN ) = WRTMP( I )
193 WITMP( KMIN ) = WITMP( I )
194 WRTMP( I ) = VRMIN
195 WITMP( I ) = VIMIN
196 ITMP = IPNT( I )
197 IPNT( I ) = IPNT( KMIN )
198 IPNT( KMIN ) = ITMP
199 60 CONTINUE
200 DO 70 I = 1, NDIM
201 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE.
202 70 CONTINUE
203 *
204 * Compute condition numbers
205 *
206 CALL DLACPY( 'F', N, N, Q, LDT, QSAV, LDT )
207 CALL DLACPY( 'F', N, N, T, LDT, TSAV1, LDT )
208 CALL DTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WRTMP, WITMP,
209 $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
210 IF( INFO.NE.0 ) THEN
211 LMAX( 3 ) = KNT
212 NINFO( 3 ) = NINFO( 3 ) + 1
213 GO TO 160
214 END IF
215 SEPTMP = SEP / VMUL
216 STMP = S
217 *
218 * Compute residuals
219 *
220 CALL DHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK,
221 $ RESULT )
222 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) )
223 IF( VMAX.GT.RMAX( 1 ) ) THEN
224 RMAX( 1 ) = VMAX
225 IF( NINFO( 1 ).EQ.0 )
226 $ LMAX( 1 ) = KNT
227 END IF
228 *
229 * Compare condition number for eigenvalue cluster
230 * taking its condition number into account
231 *
232 V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM )
233 IF( TNRM.EQ.ZERO )
234 $ V = ONE
235 IF( V.GT.SEPTMP ) THEN
236 TOL = ONE
237 ELSE
238 TOL = V / SEPTMP
239 END IF
240 IF( V.GT.SEPIN ) THEN
241 TOLIN = ONE
242 ELSE
243 TOLIN = V / SEPIN
244 END IF
245 TOL = MAX( TOL, SMLNUM / EPS )
246 TOLIN = MAX( TOLIN, SMLNUM / EPS )
247 IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN
248 VMAX = ONE / EPS
249 ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN
250 VMAX = ( SIN-TOLIN ) / ( STMP+TOL )
251 ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN
252 VMAX = ONE / EPS
253 ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN
254 VMAX = ( STMP-TOL ) / ( SIN+TOLIN )
255 ELSE
256 VMAX = ONE
257 END IF
258 IF( VMAX.GT.RMAX( 2 ) ) THEN
259 RMAX( 2 ) = VMAX
260 IF( NINFO( 2 ).EQ.0 )
261 $ LMAX( 2 ) = KNT
262 END IF
263 *
264 * Compare condition numbers for invariant subspace
265 * taking its condition number into account
266 *
267 IF( V.GT.SEPTMP*STMP ) THEN
268 TOL = SEPTMP
269 ELSE
270 TOL = V / STMP
271 END IF
272 IF( V.GT.SEPIN*SIN ) THEN
273 TOLIN = SEPIN
274 ELSE
275 TOLIN = V / SIN
276 END IF
277 TOL = MAX( TOL, SMLNUM / EPS )
278 TOLIN = MAX( TOLIN, SMLNUM / EPS )
279 IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN
280 VMAX = ONE / EPS
281 ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN
282 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL )
283 ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN
284 VMAX = ONE / EPS
285 ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN
286 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN )
287 ELSE
288 VMAX = ONE
289 END IF
290 IF( VMAX.GT.RMAX( 2 ) ) THEN
291 RMAX( 2 ) = VMAX
292 IF( NINFO( 2 ).EQ.0 )
293 $ LMAX( 2 ) = KNT
294 END IF
295 *
296 * Compare condition number for eigenvalue cluster
297 * without taking its condition number into account
298 *
299 IF( SIN.LE.DBLE( 2*N )*EPS .AND. STMP.LE.DBLE( 2*N )*EPS ) THEN
300 VMAX = ONE
301 ELSE IF( EPS*SIN.GT.STMP ) THEN
302 VMAX = ONE / EPS
303 ELSE IF( SIN.GT.STMP ) THEN
304 VMAX = SIN / STMP
305 ELSE IF( SIN.LT.EPS*STMP ) THEN
306 VMAX = ONE / EPS
307 ELSE IF( SIN.LT.STMP ) THEN
308 VMAX = STMP / SIN
309 ELSE
310 VMAX = ONE
311 END IF
312 IF( VMAX.GT.RMAX( 3 ) ) THEN
313 RMAX( 3 ) = VMAX
314 IF( NINFO( 3 ).EQ.0 )
315 $ LMAX( 3 ) = KNT
316 END IF
317 *
318 * Compare condition numbers for invariant subspace
319 * without taking its condition number into account
320 *
321 IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN
322 VMAX = ONE
323 ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN
324 VMAX = ONE / EPS
325 ELSE IF( SEPIN.GT.SEPTMP ) THEN
326 VMAX = SEPIN / SEPTMP
327 ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN
328 VMAX = ONE / EPS
329 ELSE IF( SEPIN.LT.SEPTMP ) THEN
330 VMAX = SEPTMP / SEPIN
331 ELSE
332 VMAX = ONE
333 END IF
334 IF( VMAX.GT.RMAX( 3 ) ) THEN
335 RMAX( 3 ) = VMAX
336 IF( NINFO( 3 ).EQ.0 )
337 $ LMAX( 3 ) = KNT
338 END IF
339 *
340 * Compute eigenvalue condition number only and compare
341 * Update Q
342 *
343 VMAX = ZERO
344 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
345 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
346 SEPTMP = -ONE
347 STMP = -ONE
348 CALL DTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
349 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
350 $ LIWORK, INFO )
351 IF( INFO.NE.0 ) THEN
352 LMAX( 3 ) = KNT
353 NINFO( 3 ) = NINFO( 3 ) + 1
354 GO TO 160
355 END IF
356 IF( S.NE.STMP )
357 $ VMAX = ONE / EPS
358 IF( -ONE.NE.SEPTMP )
359 $ VMAX = ONE / EPS
360 DO 90 I = 1, N
361 DO 80 J = 1, N
362 IF( TTMP( I, J ).NE.T( I, J ) )
363 $ VMAX = ONE / EPS
364 IF( QTMP( I, J ).NE.Q( I, J ) )
365 $ VMAX = ONE / EPS
366 80 CONTINUE
367 90 CONTINUE
368 *
369 * Compute invariant subspace condition number only and compare
370 * Update Q
371 *
372 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
373 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
374 SEPTMP = -ONE
375 STMP = -ONE
376 CALL DTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
377 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
378 $ LIWORK, INFO )
379 IF( INFO.NE.0 ) THEN
380 LMAX( 3 ) = KNT
381 NINFO( 3 ) = NINFO( 3 ) + 1
382 GO TO 160
383 END IF
384 IF( -ONE.NE.STMP )
385 $ VMAX = ONE / EPS
386 IF( SEP.NE.SEPTMP )
387 $ VMAX = ONE / EPS
388 DO 110 I = 1, N
389 DO 100 J = 1, N
390 IF( TTMP( I, J ).NE.T( I, J ) )
391 $ VMAX = ONE / EPS
392 IF( QTMP( I, J ).NE.Q( I, J ) )
393 $ VMAX = ONE / EPS
394 100 CONTINUE
395 110 CONTINUE
396 *
397 * Compute eigenvalue condition number only and compare
398 * Do not update Q
399 *
400 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
401 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
402 SEPTMP = -ONE
403 STMP = -ONE
404 CALL DTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
405 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
406 $ LIWORK, INFO )
407 IF( INFO.NE.0 ) THEN
408 LMAX( 3 ) = KNT
409 NINFO( 3 ) = NINFO( 3 ) + 1
410 GO TO 160
411 END IF
412 IF( S.NE.STMP )
413 $ VMAX = ONE / EPS
414 IF( -ONE.NE.SEPTMP )
415 $ VMAX = ONE / EPS
416 DO 130 I = 1, N
417 DO 120 J = 1, N
418 IF( TTMP( I, J ).NE.T( I, J ) )
419 $ VMAX = ONE / EPS
420 IF( QTMP( I, J ).NE.QSAV( I, J ) )
421 $ VMAX = ONE / EPS
422 120 CONTINUE
423 130 CONTINUE
424 *
425 * Compute invariant subspace condition number only and compare
426 * Do not update Q
427 *
428 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
429 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
430 SEPTMP = -ONE
431 STMP = -ONE
432 CALL DTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
433 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
434 $ LIWORK, INFO )
435 IF( INFO.NE.0 ) THEN
436 LMAX( 3 ) = KNT
437 NINFO( 3 ) = NINFO( 3 ) + 1
438 GO TO 160
439 END IF
440 IF( -ONE.NE.STMP )
441 $ VMAX = ONE / EPS
442 IF( SEP.NE.SEPTMP )
443 $ VMAX = ONE / EPS
444 DO 150 I = 1, N
445 DO 140 J = 1, N
446 IF( TTMP( I, J ).NE.T( I, J ) )
447 $ VMAX = ONE / EPS
448 IF( QTMP( I, J ).NE.QSAV( I, J ) )
449 $ VMAX = ONE / EPS
450 140 CONTINUE
451 150 CONTINUE
452 IF( VMAX.GT.RMAX( 1 ) ) THEN
453 RMAX( 1 ) = VMAX
454 IF( NINFO( 1 ).EQ.0 )
455 $ LMAX( 1 ) = KNT
456 END IF
457 160 CONTINUE
458 GO TO 10
459 *
460 * End of DGET38
461 *
462 END