1 SUBROUTINE CGET38( 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 REAL RMAX( 3 )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * CGET38 tests CTRSEN, 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) REAL array, dimension (3)
27 * Values of the largest test ratios.
28 * RMAX(1) = largest residuals from CHST01 or comparing
29 * different calls to CTRSEN
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 CGEHRD returns INFO nonzero on example i, LMAX(1)=i
40 * If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
41 * If CTRSEN returns INFO nonzero on example i, LMAX(3)=i
42 *
43 * NINFO (output) INTEGER array, dimension (3)
44 * NINFO(1) = No. of times CGEHRD returned INFO nonzero
45 * NINFO(2) = No. of times CHSEQR returned INFO nonzero
46 * NINFO(3) = No. of times CTRSEN 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 INTEGER LDT, LWORK
58 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
59 REAL ZERO, ONE, TWO
60 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
61 REAL EPSIN
62 PARAMETER ( EPSIN = 5.9605E-8 )
63 COMPLEX CZERO
64 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, INFO, ISCL, ISRT, ITMP, J, KMIN, M, N, NDIM
68 REAL BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
69 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VMAX, VMIN,
70 $ VMUL
71 * ..
72 * .. Local Arrays ..
73 LOGICAL SELECT( LDT )
74 INTEGER IPNT( LDT ), ISELEC( LDT )
75 REAL RESULT( 2 ), RWORK( LDT ), VAL( 3 ),
76 $ WSRT( LDT )
77 COMPLEX Q( LDT, LDT ), QSAV( LDT, LDT ),
78 $ QTMP( LDT, LDT ), T( LDT, LDT ),
79 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
80 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), W( LDT ),
81 $ WORK( LWORK ), WTMP( LDT )
82 * ..
83 * .. External Functions ..
84 REAL CLANGE, SLAMCH
85 EXTERNAL CLANGE, SLAMCH
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL CGEHRD, CHSEQR, CHST01, CLACPY, CSSCAL, CTRSEN,
89 $ CUNGHR, SLABAD
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC AIMAG, MAX, REAL, SQRT
93 * ..
94 * .. Executable Statements ..
95 *
96 EPS = SLAMCH( 'P' )
97 SMLNUM = SLAMCH( 'S' ) / EPS
98 BIGNUM = ONE / SMLNUM
99 CALL SLABAD( SMLNUM, BIGNUM )
100 *
101 * EPSIN = 2**(-24) = precision to which input data computed
102 *
103 EPS = MAX( EPS, EPSIN )
104 RMAX( 1 ) = ZERO
105 RMAX( 2 ) = ZERO
106 RMAX( 3 ) = ZERO
107 LMAX( 1 ) = 0
108 LMAX( 2 ) = 0
109 LMAX( 3 ) = 0
110 KNT = 0
111 NINFO( 1 ) = 0
112 NINFO( 2 ) = 0
113 NINFO( 3 ) = 0
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, ISRT
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 = CLANGE( 'M', N, N, TMP, LDT, RWORK )
133 DO 200 ISCL = 1, 3
134 *
135 * Scale input matrix
136 *
137 KNT = KNT + 1
138 CALL CLACPY( 'F', N, N, TMP, LDT, T, LDT )
139 VMUL = VAL( ISCL )
140 DO 30 I = 1, N
141 CALL CSSCAL( N, VMUL, T( 1, I ), 1 )
142 30 CONTINUE
143 IF( TNRM.EQ.ZERO )
144 $ VMUL = ONE
145 CALL CLACPY( 'F', N, N, T, LDT, TSAV, LDT )
146 *
147 * Compute Schur form
148 *
149 CALL CGEHRD( 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 200
155 END IF
156 *
157 * Generate unitary matrix
158 *
159 CALL CLACPY( 'L', N, N, T, LDT, Q, LDT )
160 CALL CUNGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
161 $ INFO )
162 *
163 * Compute Schur form
164 *
165 DO 50 J = 1, N - 2
166 DO 40 I = J + 2, N
167 T( I, J ) = CZERO
168 40 CONTINUE
169 50 CONTINUE
170 CALL CHSEQR( 'S', 'V', N, 1, N, T, LDT, W, Q, LDT, WORK, LWORK,
171 $ INFO )
172 IF( INFO.NE.0 ) THEN
173 LMAX( 2 ) = KNT
174 NINFO( 2 ) = NINFO( 2 ) + 1
175 GO TO 200
176 END IF
177 *
178 * Sort, select eigenvalues
179 *
180 DO 60 I = 1, N
181 IPNT( I ) = I
182 SELECT( I ) = .FALSE.
183 60 CONTINUE
184 IF( ISRT.EQ.0 ) THEN
185 DO 70 I = 1, N
186 WSRT( I ) = REAL( W( I ) )
187 70 CONTINUE
188 ELSE
189 DO 80 I = 1, N
190 WSRT( I ) = AIMAG( W( I ) )
191 80 CONTINUE
192 END IF
193 DO 100 I = 1, N - 1
194 KMIN = I
195 VMIN = WSRT( I )
196 DO 90 J = I + 1, N
197 IF( WSRT( J ).LT.VMIN ) THEN
198 KMIN = J
199 VMIN = WSRT( J )
200 END IF
201 90 CONTINUE
202 WSRT( KMIN ) = WSRT( I )
203 WSRT( I ) = VMIN
204 ITMP = IPNT( I )
205 IPNT( I ) = IPNT( KMIN )
206 IPNT( KMIN ) = ITMP
207 100 CONTINUE
208 DO 110 I = 1, NDIM
209 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE.
210 110 CONTINUE
211 *
212 * Compute condition numbers
213 *
214 CALL CLACPY( 'F', N, N, Q, LDT, QSAV, LDT )
215 CALL CLACPY( 'F', N, N, T, LDT, TSAV1, LDT )
216 CALL CTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WTMP, M, S,
217 $ SEP, WORK, LWORK, INFO )
218 IF( INFO.NE.0 ) THEN
219 LMAX( 3 ) = KNT
220 NINFO( 3 ) = NINFO( 3 ) + 1
221 GO TO 200
222 END IF
223 SEPTMP = SEP / VMUL
224 STMP = S
225 *
226 * Compute residuals
227 *
228 CALL CHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK,
229 $ RWORK, RESULT )
230 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) )
231 IF( VMAX.GT.RMAX( 1 ) ) THEN
232 RMAX( 1 ) = VMAX
233 IF( NINFO( 1 ).EQ.0 )
234 $ LMAX( 1 ) = KNT
235 END IF
236 *
237 * Compare condition number for eigenvalue cluster
238 * taking its condition number into account
239 *
240 V = MAX( TWO*REAL( N )*EPS*TNRM, SMLNUM )
241 IF( TNRM.EQ.ZERO )
242 $ V = ONE
243 IF( V.GT.SEPTMP ) THEN
244 TOL = ONE
245 ELSE
246 TOL = V / SEPTMP
247 END IF
248 IF( V.GT.SEPIN ) THEN
249 TOLIN = ONE
250 ELSE
251 TOLIN = V / SEPIN
252 END IF
253 TOL = MAX( TOL, SMLNUM / EPS )
254 TOLIN = MAX( TOLIN, SMLNUM / EPS )
255 IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN
256 VMAX = ONE / EPS
257 ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN
258 VMAX = ( SIN-TOLIN ) / ( STMP+TOL )
259 ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN
260 VMAX = ONE / EPS
261 ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN
262 VMAX = ( STMP-TOL ) / ( SIN+TOLIN )
263 ELSE
264 VMAX = ONE
265 END IF
266 IF( VMAX.GT.RMAX( 2 ) ) THEN
267 RMAX( 2 ) = VMAX
268 IF( NINFO( 2 ).EQ.0 )
269 $ LMAX( 2 ) = KNT
270 END IF
271 *
272 * Compare condition numbers for invariant subspace
273 * taking its condition number into account
274 *
275 IF( V.GT.SEPTMP*STMP ) THEN
276 TOL = SEPTMP
277 ELSE
278 TOL = V / STMP
279 END IF
280 IF( V.GT.SEPIN*SIN ) THEN
281 TOLIN = SEPIN
282 ELSE
283 TOLIN = V / SIN
284 END IF
285 TOL = MAX( TOL, SMLNUM / EPS )
286 TOLIN = MAX( TOLIN, SMLNUM / EPS )
287 IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN
288 VMAX = ONE / EPS
289 ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN
290 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL )
291 ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN
292 VMAX = ONE / EPS
293 ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN
294 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN )
295 ELSE
296 VMAX = ONE
297 END IF
298 IF( VMAX.GT.RMAX( 2 ) ) THEN
299 RMAX( 2 ) = VMAX
300 IF( NINFO( 2 ).EQ.0 )
301 $ LMAX( 2 ) = KNT
302 END IF
303 *
304 * Compare condition number for eigenvalue cluster
305 * without taking its condition number into account
306 *
307 IF( SIN.LE.REAL( 2*N )*EPS .AND. STMP.LE.REAL( 2*N )*EPS ) THEN
308 VMAX = ONE
309 ELSE IF( EPS*SIN.GT.STMP ) THEN
310 VMAX = ONE / EPS
311 ELSE IF( SIN.GT.STMP ) THEN
312 VMAX = SIN / STMP
313 ELSE IF( SIN.LT.EPS*STMP ) THEN
314 VMAX = ONE / EPS
315 ELSE IF( SIN.LT.STMP ) THEN
316 VMAX = STMP / SIN
317 ELSE
318 VMAX = ONE
319 END IF
320 IF( VMAX.GT.RMAX( 3 ) ) THEN
321 RMAX( 3 ) = VMAX
322 IF( NINFO( 3 ).EQ.0 )
323 $ LMAX( 3 ) = KNT
324 END IF
325 *
326 * Compare condition numbers for invariant subspace
327 * without taking its condition number into account
328 *
329 IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN
330 VMAX = ONE
331 ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN
332 VMAX = ONE / EPS
333 ELSE IF( SEPIN.GT.SEPTMP ) THEN
334 VMAX = SEPIN / SEPTMP
335 ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN
336 VMAX = ONE / EPS
337 ELSE IF( SEPIN.LT.SEPTMP ) THEN
338 VMAX = SEPTMP / SEPIN
339 ELSE
340 VMAX = ONE
341 END IF
342 IF( VMAX.GT.RMAX( 3 ) ) THEN
343 RMAX( 3 ) = VMAX
344 IF( NINFO( 3 ).EQ.0 )
345 $ LMAX( 3 ) = KNT
346 END IF
347 *
348 * Compute eigenvalue condition number only and compare
349 * Update Q
350 *
351 VMAX = ZERO
352 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
353 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
354 SEPTMP = -ONE
355 STMP = -ONE
356 CALL CTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
357 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
358 IF( INFO.NE.0 ) THEN
359 LMAX( 3 ) = KNT
360 NINFO( 3 ) = NINFO( 3 ) + 1
361 GO TO 200
362 END IF
363 IF( S.NE.STMP )
364 $ VMAX = ONE / EPS
365 IF( -ONE.NE.SEPTMP )
366 $ VMAX = ONE / EPS
367 DO 130 I = 1, N
368 DO 120 J = 1, N
369 IF( TTMP( I, J ).NE.T( I, J ) )
370 $ VMAX = ONE / EPS
371 IF( QTMP( I, J ).NE.Q( I, J ) )
372 $ VMAX = ONE / EPS
373 120 CONTINUE
374 130 CONTINUE
375 *
376 * Compute invariant subspace condition number only and compare
377 * Update Q
378 *
379 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
380 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
381 SEPTMP = -ONE
382 STMP = -ONE
383 CALL CTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
384 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
385 IF( INFO.NE.0 ) THEN
386 LMAX( 3 ) = KNT
387 NINFO( 3 ) = NINFO( 3 ) + 1
388 GO TO 200
389 END IF
390 IF( -ONE.NE.STMP )
391 $ VMAX = ONE / EPS
392 IF( SEP.NE.SEPTMP )
393 $ VMAX = ONE / EPS
394 DO 150 I = 1, N
395 DO 140 J = 1, N
396 IF( TTMP( I, J ).NE.T( I, J ) )
397 $ VMAX = ONE / EPS
398 IF( QTMP( I, J ).NE.Q( I, J ) )
399 $ VMAX = ONE / EPS
400 140 CONTINUE
401 150 CONTINUE
402 *
403 * Compute eigenvalue condition number only and compare
404 * Do not update Q
405 *
406 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
407 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
408 SEPTMP = -ONE
409 STMP = -ONE
410 CALL CTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
411 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
412 IF( INFO.NE.0 ) THEN
413 LMAX( 3 ) = KNT
414 NINFO( 3 ) = NINFO( 3 ) + 1
415 GO TO 200
416 END IF
417 IF( S.NE.STMP )
418 $ VMAX = ONE / EPS
419 IF( -ONE.NE.SEPTMP )
420 $ VMAX = ONE / EPS
421 DO 170 I = 1, N
422 DO 160 J = 1, N
423 IF( TTMP( I, J ).NE.T( I, J ) )
424 $ VMAX = ONE / EPS
425 IF( QTMP( I, J ).NE.QSAV( I, J ) )
426 $ VMAX = ONE / EPS
427 160 CONTINUE
428 170 CONTINUE
429 *
430 * Compute invariant subspace condition number only and compare
431 * Do not update Q
432 *
433 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
434 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
435 SEPTMP = -ONE
436 STMP = -ONE
437 CALL CTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
438 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
439 IF( INFO.NE.0 ) THEN
440 LMAX( 3 ) = KNT
441 NINFO( 3 ) = NINFO( 3 ) + 1
442 GO TO 200
443 END IF
444 IF( -ONE.NE.STMP )
445 $ VMAX = ONE / EPS
446 IF( SEP.NE.SEPTMP )
447 $ VMAX = ONE / EPS
448 DO 190 I = 1, N
449 DO 180 J = 1, N
450 IF( TTMP( I, J ).NE.T( I, J ) )
451 $ VMAX = ONE / EPS
452 IF( QTMP( I, J ).NE.QSAV( I, J ) )
453 $ VMAX = ONE / EPS
454 180 CONTINUE
455 190 CONTINUE
456 IF( VMAX.GT.RMAX( 1 ) ) THEN
457 RMAX( 1 ) = VMAX
458 IF( NINFO( 1 ).EQ.0 )
459 $ LMAX( 1 ) = KNT
460 END IF
461 200 CONTINUE
462 GO TO 10
463 *
464 * End of CGET38
465 *
466 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 REAL RMAX( 3 )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * CGET38 tests CTRSEN, 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) REAL array, dimension (3)
27 * Values of the largest test ratios.
28 * RMAX(1) = largest residuals from CHST01 or comparing
29 * different calls to CTRSEN
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 CGEHRD returns INFO nonzero on example i, LMAX(1)=i
40 * If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
41 * If CTRSEN returns INFO nonzero on example i, LMAX(3)=i
42 *
43 * NINFO (output) INTEGER array, dimension (3)
44 * NINFO(1) = No. of times CGEHRD returned INFO nonzero
45 * NINFO(2) = No. of times CHSEQR returned INFO nonzero
46 * NINFO(3) = No. of times CTRSEN 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 INTEGER LDT, LWORK
58 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
59 REAL ZERO, ONE, TWO
60 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
61 REAL EPSIN
62 PARAMETER ( EPSIN = 5.9605E-8 )
63 COMPLEX CZERO
64 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, INFO, ISCL, ISRT, ITMP, J, KMIN, M, N, NDIM
68 REAL BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
69 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VMAX, VMIN,
70 $ VMUL
71 * ..
72 * .. Local Arrays ..
73 LOGICAL SELECT( LDT )
74 INTEGER IPNT( LDT ), ISELEC( LDT )
75 REAL RESULT( 2 ), RWORK( LDT ), VAL( 3 ),
76 $ WSRT( LDT )
77 COMPLEX Q( LDT, LDT ), QSAV( LDT, LDT ),
78 $ QTMP( LDT, LDT ), T( LDT, LDT ),
79 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
80 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), W( LDT ),
81 $ WORK( LWORK ), WTMP( LDT )
82 * ..
83 * .. External Functions ..
84 REAL CLANGE, SLAMCH
85 EXTERNAL CLANGE, SLAMCH
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL CGEHRD, CHSEQR, CHST01, CLACPY, CSSCAL, CTRSEN,
89 $ CUNGHR, SLABAD
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC AIMAG, MAX, REAL, SQRT
93 * ..
94 * .. Executable Statements ..
95 *
96 EPS = SLAMCH( 'P' )
97 SMLNUM = SLAMCH( 'S' ) / EPS
98 BIGNUM = ONE / SMLNUM
99 CALL SLABAD( SMLNUM, BIGNUM )
100 *
101 * EPSIN = 2**(-24) = precision to which input data computed
102 *
103 EPS = MAX( EPS, EPSIN )
104 RMAX( 1 ) = ZERO
105 RMAX( 2 ) = ZERO
106 RMAX( 3 ) = ZERO
107 LMAX( 1 ) = 0
108 LMAX( 2 ) = 0
109 LMAX( 3 ) = 0
110 KNT = 0
111 NINFO( 1 ) = 0
112 NINFO( 2 ) = 0
113 NINFO( 3 ) = 0
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, ISRT
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 = CLANGE( 'M', N, N, TMP, LDT, RWORK )
133 DO 200 ISCL = 1, 3
134 *
135 * Scale input matrix
136 *
137 KNT = KNT + 1
138 CALL CLACPY( 'F', N, N, TMP, LDT, T, LDT )
139 VMUL = VAL( ISCL )
140 DO 30 I = 1, N
141 CALL CSSCAL( N, VMUL, T( 1, I ), 1 )
142 30 CONTINUE
143 IF( TNRM.EQ.ZERO )
144 $ VMUL = ONE
145 CALL CLACPY( 'F', N, N, T, LDT, TSAV, LDT )
146 *
147 * Compute Schur form
148 *
149 CALL CGEHRD( 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 200
155 END IF
156 *
157 * Generate unitary matrix
158 *
159 CALL CLACPY( 'L', N, N, T, LDT, Q, LDT )
160 CALL CUNGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
161 $ INFO )
162 *
163 * Compute Schur form
164 *
165 DO 50 J = 1, N - 2
166 DO 40 I = J + 2, N
167 T( I, J ) = CZERO
168 40 CONTINUE
169 50 CONTINUE
170 CALL CHSEQR( 'S', 'V', N, 1, N, T, LDT, W, Q, LDT, WORK, LWORK,
171 $ INFO )
172 IF( INFO.NE.0 ) THEN
173 LMAX( 2 ) = KNT
174 NINFO( 2 ) = NINFO( 2 ) + 1
175 GO TO 200
176 END IF
177 *
178 * Sort, select eigenvalues
179 *
180 DO 60 I = 1, N
181 IPNT( I ) = I
182 SELECT( I ) = .FALSE.
183 60 CONTINUE
184 IF( ISRT.EQ.0 ) THEN
185 DO 70 I = 1, N
186 WSRT( I ) = REAL( W( I ) )
187 70 CONTINUE
188 ELSE
189 DO 80 I = 1, N
190 WSRT( I ) = AIMAG( W( I ) )
191 80 CONTINUE
192 END IF
193 DO 100 I = 1, N - 1
194 KMIN = I
195 VMIN = WSRT( I )
196 DO 90 J = I + 1, N
197 IF( WSRT( J ).LT.VMIN ) THEN
198 KMIN = J
199 VMIN = WSRT( J )
200 END IF
201 90 CONTINUE
202 WSRT( KMIN ) = WSRT( I )
203 WSRT( I ) = VMIN
204 ITMP = IPNT( I )
205 IPNT( I ) = IPNT( KMIN )
206 IPNT( KMIN ) = ITMP
207 100 CONTINUE
208 DO 110 I = 1, NDIM
209 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE.
210 110 CONTINUE
211 *
212 * Compute condition numbers
213 *
214 CALL CLACPY( 'F', N, N, Q, LDT, QSAV, LDT )
215 CALL CLACPY( 'F', N, N, T, LDT, TSAV1, LDT )
216 CALL CTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WTMP, M, S,
217 $ SEP, WORK, LWORK, INFO )
218 IF( INFO.NE.0 ) THEN
219 LMAX( 3 ) = KNT
220 NINFO( 3 ) = NINFO( 3 ) + 1
221 GO TO 200
222 END IF
223 SEPTMP = SEP / VMUL
224 STMP = S
225 *
226 * Compute residuals
227 *
228 CALL CHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK,
229 $ RWORK, RESULT )
230 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) )
231 IF( VMAX.GT.RMAX( 1 ) ) THEN
232 RMAX( 1 ) = VMAX
233 IF( NINFO( 1 ).EQ.0 )
234 $ LMAX( 1 ) = KNT
235 END IF
236 *
237 * Compare condition number for eigenvalue cluster
238 * taking its condition number into account
239 *
240 V = MAX( TWO*REAL( N )*EPS*TNRM, SMLNUM )
241 IF( TNRM.EQ.ZERO )
242 $ V = ONE
243 IF( V.GT.SEPTMP ) THEN
244 TOL = ONE
245 ELSE
246 TOL = V / SEPTMP
247 END IF
248 IF( V.GT.SEPIN ) THEN
249 TOLIN = ONE
250 ELSE
251 TOLIN = V / SEPIN
252 END IF
253 TOL = MAX( TOL, SMLNUM / EPS )
254 TOLIN = MAX( TOLIN, SMLNUM / EPS )
255 IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN
256 VMAX = ONE / EPS
257 ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN
258 VMAX = ( SIN-TOLIN ) / ( STMP+TOL )
259 ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN
260 VMAX = ONE / EPS
261 ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN
262 VMAX = ( STMP-TOL ) / ( SIN+TOLIN )
263 ELSE
264 VMAX = ONE
265 END IF
266 IF( VMAX.GT.RMAX( 2 ) ) THEN
267 RMAX( 2 ) = VMAX
268 IF( NINFO( 2 ).EQ.0 )
269 $ LMAX( 2 ) = KNT
270 END IF
271 *
272 * Compare condition numbers for invariant subspace
273 * taking its condition number into account
274 *
275 IF( V.GT.SEPTMP*STMP ) THEN
276 TOL = SEPTMP
277 ELSE
278 TOL = V / STMP
279 END IF
280 IF( V.GT.SEPIN*SIN ) THEN
281 TOLIN = SEPIN
282 ELSE
283 TOLIN = V / SIN
284 END IF
285 TOL = MAX( TOL, SMLNUM / EPS )
286 TOLIN = MAX( TOLIN, SMLNUM / EPS )
287 IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN
288 VMAX = ONE / EPS
289 ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN
290 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL )
291 ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN
292 VMAX = ONE / EPS
293 ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN
294 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN )
295 ELSE
296 VMAX = ONE
297 END IF
298 IF( VMAX.GT.RMAX( 2 ) ) THEN
299 RMAX( 2 ) = VMAX
300 IF( NINFO( 2 ).EQ.0 )
301 $ LMAX( 2 ) = KNT
302 END IF
303 *
304 * Compare condition number for eigenvalue cluster
305 * without taking its condition number into account
306 *
307 IF( SIN.LE.REAL( 2*N )*EPS .AND. STMP.LE.REAL( 2*N )*EPS ) THEN
308 VMAX = ONE
309 ELSE IF( EPS*SIN.GT.STMP ) THEN
310 VMAX = ONE / EPS
311 ELSE IF( SIN.GT.STMP ) THEN
312 VMAX = SIN / STMP
313 ELSE IF( SIN.LT.EPS*STMP ) THEN
314 VMAX = ONE / EPS
315 ELSE IF( SIN.LT.STMP ) THEN
316 VMAX = STMP / SIN
317 ELSE
318 VMAX = ONE
319 END IF
320 IF( VMAX.GT.RMAX( 3 ) ) THEN
321 RMAX( 3 ) = VMAX
322 IF( NINFO( 3 ).EQ.0 )
323 $ LMAX( 3 ) = KNT
324 END IF
325 *
326 * Compare condition numbers for invariant subspace
327 * without taking its condition number into account
328 *
329 IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN
330 VMAX = ONE
331 ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN
332 VMAX = ONE / EPS
333 ELSE IF( SEPIN.GT.SEPTMP ) THEN
334 VMAX = SEPIN / SEPTMP
335 ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN
336 VMAX = ONE / EPS
337 ELSE IF( SEPIN.LT.SEPTMP ) THEN
338 VMAX = SEPTMP / SEPIN
339 ELSE
340 VMAX = ONE
341 END IF
342 IF( VMAX.GT.RMAX( 3 ) ) THEN
343 RMAX( 3 ) = VMAX
344 IF( NINFO( 3 ).EQ.0 )
345 $ LMAX( 3 ) = KNT
346 END IF
347 *
348 * Compute eigenvalue condition number only and compare
349 * Update Q
350 *
351 VMAX = ZERO
352 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
353 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
354 SEPTMP = -ONE
355 STMP = -ONE
356 CALL CTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
357 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
358 IF( INFO.NE.0 ) THEN
359 LMAX( 3 ) = KNT
360 NINFO( 3 ) = NINFO( 3 ) + 1
361 GO TO 200
362 END IF
363 IF( S.NE.STMP )
364 $ VMAX = ONE / EPS
365 IF( -ONE.NE.SEPTMP )
366 $ VMAX = ONE / EPS
367 DO 130 I = 1, N
368 DO 120 J = 1, N
369 IF( TTMP( I, J ).NE.T( I, J ) )
370 $ VMAX = ONE / EPS
371 IF( QTMP( I, J ).NE.Q( I, J ) )
372 $ VMAX = ONE / EPS
373 120 CONTINUE
374 130 CONTINUE
375 *
376 * Compute invariant subspace condition number only and compare
377 * Update Q
378 *
379 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
380 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
381 SEPTMP = -ONE
382 STMP = -ONE
383 CALL CTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
384 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
385 IF( INFO.NE.0 ) THEN
386 LMAX( 3 ) = KNT
387 NINFO( 3 ) = NINFO( 3 ) + 1
388 GO TO 200
389 END IF
390 IF( -ONE.NE.STMP )
391 $ VMAX = ONE / EPS
392 IF( SEP.NE.SEPTMP )
393 $ VMAX = ONE / EPS
394 DO 150 I = 1, N
395 DO 140 J = 1, N
396 IF( TTMP( I, J ).NE.T( I, J ) )
397 $ VMAX = ONE / EPS
398 IF( QTMP( I, J ).NE.Q( I, J ) )
399 $ VMAX = ONE / EPS
400 140 CONTINUE
401 150 CONTINUE
402 *
403 * Compute eigenvalue condition number only and compare
404 * Do not update Q
405 *
406 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
407 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
408 SEPTMP = -ONE
409 STMP = -ONE
410 CALL CTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
411 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
412 IF( INFO.NE.0 ) THEN
413 LMAX( 3 ) = KNT
414 NINFO( 3 ) = NINFO( 3 ) + 1
415 GO TO 200
416 END IF
417 IF( S.NE.STMP )
418 $ VMAX = ONE / EPS
419 IF( -ONE.NE.SEPTMP )
420 $ VMAX = ONE / EPS
421 DO 170 I = 1, N
422 DO 160 J = 1, N
423 IF( TTMP( I, J ).NE.T( I, J ) )
424 $ VMAX = ONE / EPS
425 IF( QTMP( I, J ).NE.QSAV( I, J ) )
426 $ VMAX = ONE / EPS
427 160 CONTINUE
428 170 CONTINUE
429 *
430 * Compute invariant subspace condition number only and compare
431 * Do not update Q
432 *
433 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT )
434 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT )
435 SEPTMP = -ONE
436 STMP = -ONE
437 CALL CTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP,
438 $ M, STMP, SEPTMP, WORK, LWORK, INFO )
439 IF( INFO.NE.0 ) THEN
440 LMAX( 3 ) = KNT
441 NINFO( 3 ) = NINFO( 3 ) + 1
442 GO TO 200
443 END IF
444 IF( -ONE.NE.STMP )
445 $ VMAX = ONE / EPS
446 IF( SEP.NE.SEPTMP )
447 $ VMAX = ONE / EPS
448 DO 190 I = 1, N
449 DO 180 J = 1, N
450 IF( TTMP( I, J ).NE.T( I, J ) )
451 $ VMAX = ONE / EPS
452 IF( QTMP( I, J ).NE.QSAV( I, J ) )
453 $ VMAX = ONE / EPS
454 180 CONTINUE
455 190 CONTINUE
456 IF( VMAX.GT.RMAX( 1 ) ) THEN
457 RMAX( 1 ) = VMAX
458 IF( NINFO( 1 ).EQ.0 )
459 $ LMAX( 1 ) = KNT
460 END IF
461 200 CONTINUE
462 GO TO 10
463 *
464 * End of CGET38
465 *
466 END