1 SUBROUTINE CGET37( 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 * CGET37 tests CTRSNA, a routine for estimating condition numbers of
19 * eigenvalues and/or right eigenvectors of a matrix.
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 * Value of the largest test ratio.
28 * RMAX(1) = largest ratio comparing different calls to CTRSNA
29 * RMAX(2) = largest error in reciprocal condition
30 * numbers taking their conditioning into account
31 * RMAX(3) = largest error in reciprocal condition
32 * numbers not taking their conditioning into
33 * account (may be larger than RMAX(2))
34 *
35 * LMAX (output) INTEGER array, dimension (3)
36 * LMAX(i) is example number where largest test ratio
37 * RMAX(i) is achieved. Also:
38 * If CGEHRD returns INFO nonzero on example i, LMAX(1)=i
39 * If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
40 * If CTRSNA returns INFO nonzero on example i, LMAX(3)=i
41 *
42 * NINFO (output) INTEGER array, dimension (3)
43 * NINFO(1) = No. of times CGEHRD returned INFO nonzero
44 * NINFO(2) = No. of times CHSEQR returned INFO nonzero
45 * NINFO(3) = No. of times CTRSNA returned INFO nonzero
46 *
47 * KNT (output) INTEGER
48 * Total number of examples tested.
49 *
50 * NIN (input) INTEGER
51 * Input logical unit number
52 *
53 * =====================================================================
54 *
55 * .. Parameters ..
56 REAL ZERO, ONE, TWO
57 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
58 REAL EPSIN
59 PARAMETER ( EPSIN = 5.9605E-8 )
60 INTEGER LDT, LWORK
61 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
62 * ..
63 * .. Local Scalars ..
64 INTEGER I, ICMP, INFO, ISCL, ISRT, J, KMIN, M, N
65 REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
66 $ VCMIN, VMAX, VMIN, VMUL
67 * ..
68 * .. Local Arrays ..
69 LOGICAL SELECT( LDT )
70 INTEGER LCMP( 3 )
71 REAL DUM( 1 ), RWORK( 2*LDT ), S( LDT ), SEP( LDT ),
72 $ SEPIN( LDT ), SEPTMP( LDT ), SIN( LDT ),
73 $ STMP( LDT ), VAL( 3 ), WIIN( LDT ),
74 $ WRIN( LDT ), WSRT( LDT )
75 COMPLEX CDUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
76 $ T( LDT, LDT ), TMP( LDT, LDT ), W( LDT ),
77 $ WORK( LWORK ), WTMP( LDT )
78 * ..
79 * .. External Functions ..
80 REAL CLANGE, SLAMCH
81 EXTERNAL CLANGE, SLAMCH
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL CCOPY, CGEHRD, CHSEQR, CLACPY, CSSCAL, CTREVC,
85 $ CTRSNA, SCOPY, SLABAD, SSCAL
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC AIMAG, MAX, REAL, SQRT
89 * ..
90 * .. Executable Statements ..
91 *
92 EPS = SLAMCH( 'P' )
93 SMLNUM = SLAMCH( 'S' ) / EPS
94 BIGNUM = ONE / SMLNUM
95 CALL SLABAD( SMLNUM, BIGNUM )
96 *
97 * EPSIN = 2**(-24) = precision to which input data computed
98 *
99 EPS = MAX( EPS, EPSIN )
100 RMAX( 1 ) = ZERO
101 RMAX( 2 ) = ZERO
102 RMAX( 3 ) = ZERO
103 LMAX( 1 ) = 0
104 LMAX( 2 ) = 0
105 LMAX( 3 ) = 0
106 KNT = 0
107 NINFO( 1 ) = 0
108 NINFO( 2 ) = 0
109 NINFO( 3 ) = 0
110 VAL( 1 ) = SQRT( SMLNUM )
111 VAL( 2 ) = ONE
112 VAL( 3 ) = SQRT( BIGNUM )
113 *
114 * Read input data until N=0. Assume input eigenvalues are sorted
115 * lexicographically (increasing by real part if ISRT = 0,
116 * increasing by imaginary part if ISRT = 1)
117 *
118 10 CONTINUE
119 READ( NIN, FMT = * )N, ISRT
120 IF( N.EQ.0 )
121 $ RETURN
122 DO 20 I = 1, N
123 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
124 20 CONTINUE
125 DO 30 I = 1, N
126 READ( NIN, FMT = * )WRIN( I ), WIIN( I ), SIN( I ), SEPIN( I )
127 30 CONTINUE
128 TNRM = CLANGE( 'M', N, N, TMP, LDT, RWORK )
129 DO 260 ISCL = 1, 3
130 *
131 * Scale input matrix
132 *
133 KNT = KNT + 1
134 CALL CLACPY( 'F', N, N, TMP, LDT, T, LDT )
135 VMUL = VAL( ISCL )
136 DO 40 I = 1, N
137 CALL CSSCAL( N, VMUL, T( 1, I ), 1 )
138 40 CONTINUE
139 IF( TNRM.EQ.ZERO )
140 $ VMUL = ONE
141 *
142 * Compute eigenvalues and eigenvectors
143 *
144 CALL CGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
145 $ INFO )
146 IF( INFO.NE.0 ) THEN
147 LMAX( 1 ) = KNT
148 NINFO( 1 ) = NINFO( 1 ) + 1
149 GO TO 260
150 END IF
151 DO 60 J = 1, N - 2
152 DO 50 I = J + 2, N
153 T( I, J ) = ZERO
154 50 CONTINUE
155 60 CONTINUE
156 *
157 * Compute Schur form
158 *
159 CALL CHSEQR( 'S', 'N', N, 1, N, T, LDT, W, CDUM, 1, WORK,
160 $ LWORK, INFO )
161 IF( INFO.NE.0 ) THEN
162 LMAX( 2 ) = KNT
163 NINFO( 2 ) = NINFO( 2 ) + 1
164 GO TO 260
165 END IF
166 *
167 * Compute eigenvectors
168 *
169 DO 70 I = 1, N
170 SELECT( I ) = .TRUE.
171 70 CONTINUE
172 CALL CTREVC( 'B', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT, N,
173 $ M, WORK, RWORK, INFO )
174 *
175 * Compute condition numbers
176 *
177 CALL CTRSNA( 'B', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT, S,
178 $ SEP, N, M, WORK, N, RWORK, INFO )
179 IF( INFO.NE.0 ) THEN
180 LMAX( 3 ) = KNT
181 NINFO( 3 ) = NINFO( 3 ) + 1
182 GO TO 260
183 END IF
184 *
185 * Sort eigenvalues and condition numbers lexicographically
186 * to compare with inputs
187 *
188 CALL CCOPY( N, W, 1, WTMP, 1 )
189 IF( ISRT.EQ.0 ) THEN
190 *
191 * Sort by increasing real part
192 *
193 DO 80 I = 1, N
194 WSRT( I ) = REAL( W( I ) )
195 80 CONTINUE
196 ELSE
197 *
198 * Sort by increasing imaginary part
199 *
200 DO 90 I = 1, N
201 WSRT( I ) = AIMAG( W( I ) )
202 90 CONTINUE
203 END IF
204 CALL SCOPY( N, S, 1, STMP, 1 )
205 CALL SCOPY( N, SEP, 1, SEPTMP, 1 )
206 CALL SSCAL( N, ONE / VMUL, SEPTMP, 1 )
207 DO 110 I = 1, N - 1
208 KMIN = I
209 VMIN = WSRT( I )
210 DO 100 J = I + 1, N
211 IF( WSRT( J ).LT.VMIN ) THEN
212 KMIN = J
213 VMIN = WSRT( J )
214 END IF
215 100 CONTINUE
216 WSRT( KMIN ) = WSRT( I )
217 WSRT( I ) = VMIN
218 VCMIN = WTMP( I )
219 WTMP( I ) = W( KMIN )
220 WTMP( KMIN ) = VCMIN
221 VMIN = STMP( KMIN )
222 STMP( KMIN ) = STMP( I )
223 STMP( I ) = VMIN
224 VMIN = SEPTMP( KMIN )
225 SEPTMP( KMIN ) = SEPTMP( I )
226 SEPTMP( I ) = VMIN
227 110 CONTINUE
228 *
229 * Compare condition numbers for eigenvalues
230 * taking their condition numbers into account
231 *
232 V = MAX( TWO*REAL( N )*EPS*TNRM, SMLNUM )
233 IF( TNRM.EQ.ZERO )
234 $ V = ONE
235 DO 120 I = 1, N
236 IF( V.GT.SEPTMP( I ) ) THEN
237 TOL = ONE
238 ELSE
239 TOL = V / SEPTMP( I )
240 END IF
241 IF( V.GT.SEPIN( I ) ) THEN
242 TOLIN = ONE
243 ELSE
244 TOLIN = V / SEPIN( I )
245 END IF
246 TOL = MAX( TOL, SMLNUM / EPS )
247 TOLIN = MAX( TOLIN, SMLNUM / EPS )
248 IF( EPS*( SIN( I )-TOLIN ).GT.STMP( I )+TOL ) THEN
249 VMAX = ONE / EPS
250 ELSE IF( SIN( I )-TOLIN.GT.STMP( I )+TOL ) THEN
251 VMAX = ( SIN( I )-TOLIN ) / ( STMP( I )+TOL )
252 ELSE IF( SIN( I )+TOLIN.LT.EPS*( STMP( I )-TOL ) ) THEN
253 VMAX = ONE / EPS
254 ELSE IF( SIN( I )+TOLIN.LT.STMP( I )-TOL ) THEN
255 VMAX = ( STMP( I )-TOL ) / ( SIN( I )+TOLIN )
256 ELSE
257 VMAX = ONE
258 END IF
259 IF( VMAX.GT.RMAX( 2 ) ) THEN
260 RMAX( 2 ) = VMAX
261 IF( NINFO( 2 ).EQ.0 )
262 $ LMAX( 2 ) = KNT
263 END IF
264 120 CONTINUE
265 *
266 * Compare condition numbers for eigenvectors
267 * taking their condition numbers into account
268 *
269 DO 130 I = 1, N
270 IF( V.GT.SEPTMP( I )*STMP( I ) ) THEN
271 TOL = SEPTMP( I )
272 ELSE
273 TOL = V / STMP( I )
274 END IF
275 IF( V.GT.SEPIN( I )*SIN( I ) ) THEN
276 TOLIN = SEPIN( I )
277 ELSE
278 TOLIN = V / SIN( I )
279 END IF
280 TOL = MAX( TOL, SMLNUM / EPS )
281 TOLIN = MAX( TOLIN, SMLNUM / EPS )
282 IF( EPS*( SEPIN( I )-TOLIN ).GT.SEPTMP( I )+TOL ) THEN
283 VMAX = ONE / EPS
284 ELSE IF( SEPIN( I )-TOLIN.GT.SEPTMP( I )+TOL ) THEN
285 VMAX = ( SEPIN( I )-TOLIN ) / ( SEPTMP( I )+TOL )
286 ELSE IF( SEPIN( I )+TOLIN.LT.EPS*( SEPTMP( I )-TOL ) ) THEN
287 VMAX = ONE / EPS
288 ELSE IF( SEPIN( I )+TOLIN.LT.SEPTMP( I )-TOL ) THEN
289 VMAX = ( SEPTMP( I )-TOL ) / ( SEPIN( I )+TOLIN )
290 ELSE
291 VMAX = ONE
292 END IF
293 IF( VMAX.GT.RMAX( 2 ) ) THEN
294 RMAX( 2 ) = VMAX
295 IF( NINFO( 2 ).EQ.0 )
296 $ LMAX( 2 ) = KNT
297 END IF
298 130 CONTINUE
299 *
300 * Compare condition numbers for eigenvalues
301 * without taking their condition numbers into account
302 *
303 DO 140 I = 1, N
304 IF( SIN( I ).LE.REAL( 2*N )*EPS .AND. STMP( I ).LE.
305 $ REAL( 2*N )*EPS ) THEN
306 VMAX = ONE
307 ELSE IF( EPS*SIN( I ).GT.STMP( I ) ) THEN
308 VMAX = ONE / EPS
309 ELSE IF( SIN( I ).GT.STMP( I ) ) THEN
310 VMAX = SIN( I ) / STMP( I )
311 ELSE IF( SIN( I ).LT.EPS*STMP( I ) ) THEN
312 VMAX = ONE / EPS
313 ELSE IF( SIN( I ).LT.STMP( I ) ) THEN
314 VMAX = STMP( I ) / SIN( I )
315 ELSE
316 VMAX = ONE
317 END IF
318 IF( VMAX.GT.RMAX( 3 ) ) THEN
319 RMAX( 3 ) = VMAX
320 IF( NINFO( 3 ).EQ.0 )
321 $ LMAX( 3 ) = KNT
322 END IF
323 140 CONTINUE
324 *
325 * Compare condition numbers for eigenvectors
326 * without taking their condition numbers into account
327 *
328 DO 150 I = 1, N
329 IF( SEPIN( I ).LE.V .AND. SEPTMP( I ).LE.V ) THEN
330 VMAX = ONE
331 ELSE IF( EPS*SEPIN( I ).GT.SEPTMP( I ) ) THEN
332 VMAX = ONE / EPS
333 ELSE IF( SEPIN( I ).GT.SEPTMP( I ) ) THEN
334 VMAX = SEPIN( I ) / SEPTMP( I )
335 ELSE IF( SEPIN( I ).LT.EPS*SEPTMP( I ) ) THEN
336 VMAX = ONE / EPS
337 ELSE IF( SEPIN( I ).LT.SEPTMP( I ) ) THEN
338 VMAX = SEPTMP( I ) / SEPIN( I )
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 150 CONTINUE
348 *
349 * Compute eigenvalue condition numbers only and compare
350 *
351 VMAX = ZERO
352 DUM( 1 ) = -ONE
353 CALL SCOPY( N, DUM, 0, STMP, 1 )
354 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
355 CALL CTRSNA( 'E', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT,
356 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
357 IF( INFO.NE.0 ) THEN
358 LMAX( 3 ) = KNT
359 NINFO( 3 ) = NINFO( 3 ) + 1
360 GO TO 260
361 END IF
362 DO 160 I = 1, N
363 IF( STMP( I ).NE.S( I ) )
364 $ VMAX = ONE / EPS
365 IF( SEPTMP( I ).NE.DUM( 1 ) )
366 $ VMAX = ONE / EPS
367 160 CONTINUE
368 *
369 * Compute eigenvector condition numbers only and compare
370 *
371 CALL SCOPY( N, DUM, 0, STMP, 1 )
372 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
373 CALL CTRSNA( 'V', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT,
374 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
375 IF( INFO.NE.0 ) THEN
376 LMAX( 3 ) = KNT
377 NINFO( 3 ) = NINFO( 3 ) + 1
378 GO TO 260
379 END IF
380 DO 170 I = 1, N
381 IF( STMP( I ).NE.DUM( 1 ) )
382 $ VMAX = ONE / EPS
383 IF( SEPTMP( I ).NE.SEP( I ) )
384 $ VMAX = ONE / EPS
385 170 CONTINUE
386 *
387 * Compute all condition numbers using SELECT and compare
388 *
389 DO 180 I = 1, N
390 SELECT( I ) = .TRUE.
391 180 CONTINUE
392 CALL SCOPY( N, DUM, 0, STMP, 1 )
393 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
394 CALL CTRSNA( 'B', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
395 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
396 IF( INFO.NE.0 ) THEN
397 LMAX( 3 ) = KNT
398 NINFO( 3 ) = NINFO( 3 ) + 1
399 GO TO 260
400 END IF
401 DO 190 I = 1, N
402 IF( SEPTMP( I ).NE.SEP( I ) )
403 $ VMAX = ONE / EPS
404 IF( STMP( I ).NE.S( I ) )
405 $ VMAX = ONE / EPS
406 190 CONTINUE
407 *
408 * Compute eigenvalue condition numbers using SELECT and compare
409 *
410 CALL SCOPY( N, DUM, 0, STMP, 1 )
411 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
412 CALL CTRSNA( 'E', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
413 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
414 IF( INFO.NE.0 ) THEN
415 LMAX( 3 ) = KNT
416 NINFO( 3 ) = NINFO( 3 ) + 1
417 GO TO 260
418 END IF
419 DO 200 I = 1, N
420 IF( STMP( I ).NE.S( I ) )
421 $ VMAX = ONE / EPS
422 IF( SEPTMP( I ).NE.DUM( 1 ) )
423 $ VMAX = ONE / EPS
424 200 CONTINUE
425 *
426 * Compute eigenvector condition numbers using SELECT and compare
427 *
428 CALL SCOPY( N, DUM, 0, STMP, 1 )
429 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
430 CALL CTRSNA( 'V', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
431 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
432 IF( INFO.NE.0 ) THEN
433 LMAX( 3 ) = KNT
434 NINFO( 3 ) = NINFO( 3 ) + 1
435 GO TO 260
436 END IF
437 DO 210 I = 1, N
438 IF( STMP( I ).NE.DUM( 1 ) )
439 $ VMAX = ONE / EPS
440 IF( SEPTMP( I ).NE.SEP( I ) )
441 $ VMAX = ONE / EPS
442 210 CONTINUE
443 IF( VMAX.GT.RMAX( 1 ) ) THEN
444 RMAX( 1 ) = VMAX
445 IF( NINFO( 1 ).EQ.0 )
446 $ LMAX( 1 ) = KNT
447 END IF
448 *
449 * Select second and next to last eigenvalues
450 *
451 DO 220 I = 1, N
452 SELECT( I ) = .FALSE.
453 220 CONTINUE
454 ICMP = 0
455 IF( N.GT.1 ) THEN
456 ICMP = 1
457 LCMP( 1 ) = 2
458 SELECT( 2 ) = .TRUE.
459 CALL CCOPY( N, RE( 1, 2 ), 1, RE( 1, 1 ), 1 )
460 CALL CCOPY( N, LE( 1, 2 ), 1, LE( 1, 1 ), 1 )
461 END IF
462 IF( N.GT.3 ) THEN
463 ICMP = 2
464 LCMP( 2 ) = N - 1
465 SELECT( N-1 ) = .TRUE.
466 CALL CCOPY( N, RE( 1, N-1 ), 1, RE( 1, 2 ), 1 )
467 CALL CCOPY( N, LE( 1, N-1 ), 1, LE( 1, 2 ), 1 )
468 END IF
469 *
470 * Compute all selected condition numbers
471 *
472 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
473 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
474 CALL CTRSNA( 'B', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
475 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
476 IF( INFO.NE.0 ) THEN
477 LMAX( 3 ) = KNT
478 NINFO( 3 ) = NINFO( 3 ) + 1
479 GO TO 260
480 END IF
481 DO 230 I = 1, ICMP
482 J = LCMP( I )
483 IF( SEPTMP( I ).NE.SEP( J ) )
484 $ VMAX = ONE / EPS
485 IF( STMP( I ).NE.S( J ) )
486 $ VMAX = ONE / EPS
487 230 CONTINUE
488 *
489 * Compute selected eigenvalue condition numbers
490 *
491 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
492 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
493 CALL CTRSNA( 'E', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
494 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
495 IF( INFO.NE.0 ) THEN
496 LMAX( 3 ) = KNT
497 NINFO( 3 ) = NINFO( 3 ) + 1
498 GO TO 260
499 END IF
500 DO 240 I = 1, ICMP
501 J = LCMP( I )
502 IF( STMP( I ).NE.S( J ) )
503 $ VMAX = ONE / EPS
504 IF( SEPTMP( I ).NE.DUM( 1 ) )
505 $ VMAX = ONE / EPS
506 240 CONTINUE
507 *
508 * Compute selected eigenvector condition numbers
509 *
510 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
511 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
512 CALL CTRSNA( 'V', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
513 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
514 IF( INFO.NE.0 ) THEN
515 LMAX( 3 ) = KNT
516 NINFO( 3 ) = NINFO( 3 ) + 1
517 GO TO 260
518 END IF
519 DO 250 I = 1, ICMP
520 J = LCMP( I )
521 IF( STMP( I ).NE.DUM( 1 ) )
522 $ VMAX = ONE / EPS
523 IF( SEPTMP( I ).NE.SEP( J ) )
524 $ VMAX = ONE / EPS
525 250 CONTINUE
526 IF( VMAX.GT.RMAX( 1 ) ) THEN
527 RMAX( 1 ) = VMAX
528 IF( NINFO( 1 ).EQ.0 )
529 $ LMAX( 1 ) = KNT
530 END IF
531 260 CONTINUE
532 GO TO 10
533 *
534 * End of CGET37
535 *
536 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 * CGET37 tests CTRSNA, a routine for estimating condition numbers of
19 * eigenvalues and/or right eigenvectors of a matrix.
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 * Value of the largest test ratio.
28 * RMAX(1) = largest ratio comparing different calls to CTRSNA
29 * RMAX(2) = largest error in reciprocal condition
30 * numbers taking their conditioning into account
31 * RMAX(3) = largest error in reciprocal condition
32 * numbers not taking their conditioning into
33 * account (may be larger than RMAX(2))
34 *
35 * LMAX (output) INTEGER array, dimension (3)
36 * LMAX(i) is example number where largest test ratio
37 * RMAX(i) is achieved. Also:
38 * If CGEHRD returns INFO nonzero on example i, LMAX(1)=i
39 * If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
40 * If CTRSNA returns INFO nonzero on example i, LMAX(3)=i
41 *
42 * NINFO (output) INTEGER array, dimension (3)
43 * NINFO(1) = No. of times CGEHRD returned INFO nonzero
44 * NINFO(2) = No. of times CHSEQR returned INFO nonzero
45 * NINFO(3) = No. of times CTRSNA returned INFO nonzero
46 *
47 * KNT (output) INTEGER
48 * Total number of examples tested.
49 *
50 * NIN (input) INTEGER
51 * Input logical unit number
52 *
53 * =====================================================================
54 *
55 * .. Parameters ..
56 REAL ZERO, ONE, TWO
57 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
58 REAL EPSIN
59 PARAMETER ( EPSIN = 5.9605E-8 )
60 INTEGER LDT, LWORK
61 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) )
62 * ..
63 * .. Local Scalars ..
64 INTEGER I, ICMP, INFO, ISCL, ISRT, J, KMIN, M, N
65 REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
66 $ VCMIN, VMAX, VMIN, VMUL
67 * ..
68 * .. Local Arrays ..
69 LOGICAL SELECT( LDT )
70 INTEGER LCMP( 3 )
71 REAL DUM( 1 ), RWORK( 2*LDT ), S( LDT ), SEP( LDT ),
72 $ SEPIN( LDT ), SEPTMP( LDT ), SIN( LDT ),
73 $ STMP( LDT ), VAL( 3 ), WIIN( LDT ),
74 $ WRIN( LDT ), WSRT( LDT )
75 COMPLEX CDUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
76 $ T( LDT, LDT ), TMP( LDT, LDT ), W( LDT ),
77 $ WORK( LWORK ), WTMP( LDT )
78 * ..
79 * .. External Functions ..
80 REAL CLANGE, SLAMCH
81 EXTERNAL CLANGE, SLAMCH
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL CCOPY, CGEHRD, CHSEQR, CLACPY, CSSCAL, CTREVC,
85 $ CTRSNA, SCOPY, SLABAD, SSCAL
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC AIMAG, MAX, REAL, SQRT
89 * ..
90 * .. Executable Statements ..
91 *
92 EPS = SLAMCH( 'P' )
93 SMLNUM = SLAMCH( 'S' ) / EPS
94 BIGNUM = ONE / SMLNUM
95 CALL SLABAD( SMLNUM, BIGNUM )
96 *
97 * EPSIN = 2**(-24) = precision to which input data computed
98 *
99 EPS = MAX( EPS, EPSIN )
100 RMAX( 1 ) = ZERO
101 RMAX( 2 ) = ZERO
102 RMAX( 3 ) = ZERO
103 LMAX( 1 ) = 0
104 LMAX( 2 ) = 0
105 LMAX( 3 ) = 0
106 KNT = 0
107 NINFO( 1 ) = 0
108 NINFO( 2 ) = 0
109 NINFO( 3 ) = 0
110 VAL( 1 ) = SQRT( SMLNUM )
111 VAL( 2 ) = ONE
112 VAL( 3 ) = SQRT( BIGNUM )
113 *
114 * Read input data until N=0. Assume input eigenvalues are sorted
115 * lexicographically (increasing by real part if ISRT = 0,
116 * increasing by imaginary part if ISRT = 1)
117 *
118 10 CONTINUE
119 READ( NIN, FMT = * )N, ISRT
120 IF( N.EQ.0 )
121 $ RETURN
122 DO 20 I = 1, N
123 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
124 20 CONTINUE
125 DO 30 I = 1, N
126 READ( NIN, FMT = * )WRIN( I ), WIIN( I ), SIN( I ), SEPIN( I )
127 30 CONTINUE
128 TNRM = CLANGE( 'M', N, N, TMP, LDT, RWORK )
129 DO 260 ISCL = 1, 3
130 *
131 * Scale input matrix
132 *
133 KNT = KNT + 1
134 CALL CLACPY( 'F', N, N, TMP, LDT, T, LDT )
135 VMUL = VAL( ISCL )
136 DO 40 I = 1, N
137 CALL CSSCAL( N, VMUL, T( 1, I ), 1 )
138 40 CONTINUE
139 IF( TNRM.EQ.ZERO )
140 $ VMUL = ONE
141 *
142 * Compute eigenvalues and eigenvectors
143 *
144 CALL CGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
145 $ INFO )
146 IF( INFO.NE.0 ) THEN
147 LMAX( 1 ) = KNT
148 NINFO( 1 ) = NINFO( 1 ) + 1
149 GO TO 260
150 END IF
151 DO 60 J = 1, N - 2
152 DO 50 I = J + 2, N
153 T( I, J ) = ZERO
154 50 CONTINUE
155 60 CONTINUE
156 *
157 * Compute Schur form
158 *
159 CALL CHSEQR( 'S', 'N', N, 1, N, T, LDT, W, CDUM, 1, WORK,
160 $ LWORK, INFO )
161 IF( INFO.NE.0 ) THEN
162 LMAX( 2 ) = KNT
163 NINFO( 2 ) = NINFO( 2 ) + 1
164 GO TO 260
165 END IF
166 *
167 * Compute eigenvectors
168 *
169 DO 70 I = 1, N
170 SELECT( I ) = .TRUE.
171 70 CONTINUE
172 CALL CTREVC( 'B', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT, N,
173 $ M, WORK, RWORK, INFO )
174 *
175 * Compute condition numbers
176 *
177 CALL CTRSNA( 'B', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT, S,
178 $ SEP, N, M, WORK, N, RWORK, INFO )
179 IF( INFO.NE.0 ) THEN
180 LMAX( 3 ) = KNT
181 NINFO( 3 ) = NINFO( 3 ) + 1
182 GO TO 260
183 END IF
184 *
185 * Sort eigenvalues and condition numbers lexicographically
186 * to compare with inputs
187 *
188 CALL CCOPY( N, W, 1, WTMP, 1 )
189 IF( ISRT.EQ.0 ) THEN
190 *
191 * Sort by increasing real part
192 *
193 DO 80 I = 1, N
194 WSRT( I ) = REAL( W( I ) )
195 80 CONTINUE
196 ELSE
197 *
198 * Sort by increasing imaginary part
199 *
200 DO 90 I = 1, N
201 WSRT( I ) = AIMAG( W( I ) )
202 90 CONTINUE
203 END IF
204 CALL SCOPY( N, S, 1, STMP, 1 )
205 CALL SCOPY( N, SEP, 1, SEPTMP, 1 )
206 CALL SSCAL( N, ONE / VMUL, SEPTMP, 1 )
207 DO 110 I = 1, N - 1
208 KMIN = I
209 VMIN = WSRT( I )
210 DO 100 J = I + 1, N
211 IF( WSRT( J ).LT.VMIN ) THEN
212 KMIN = J
213 VMIN = WSRT( J )
214 END IF
215 100 CONTINUE
216 WSRT( KMIN ) = WSRT( I )
217 WSRT( I ) = VMIN
218 VCMIN = WTMP( I )
219 WTMP( I ) = W( KMIN )
220 WTMP( KMIN ) = VCMIN
221 VMIN = STMP( KMIN )
222 STMP( KMIN ) = STMP( I )
223 STMP( I ) = VMIN
224 VMIN = SEPTMP( KMIN )
225 SEPTMP( KMIN ) = SEPTMP( I )
226 SEPTMP( I ) = VMIN
227 110 CONTINUE
228 *
229 * Compare condition numbers for eigenvalues
230 * taking their condition numbers into account
231 *
232 V = MAX( TWO*REAL( N )*EPS*TNRM, SMLNUM )
233 IF( TNRM.EQ.ZERO )
234 $ V = ONE
235 DO 120 I = 1, N
236 IF( V.GT.SEPTMP( I ) ) THEN
237 TOL = ONE
238 ELSE
239 TOL = V / SEPTMP( I )
240 END IF
241 IF( V.GT.SEPIN( I ) ) THEN
242 TOLIN = ONE
243 ELSE
244 TOLIN = V / SEPIN( I )
245 END IF
246 TOL = MAX( TOL, SMLNUM / EPS )
247 TOLIN = MAX( TOLIN, SMLNUM / EPS )
248 IF( EPS*( SIN( I )-TOLIN ).GT.STMP( I )+TOL ) THEN
249 VMAX = ONE / EPS
250 ELSE IF( SIN( I )-TOLIN.GT.STMP( I )+TOL ) THEN
251 VMAX = ( SIN( I )-TOLIN ) / ( STMP( I )+TOL )
252 ELSE IF( SIN( I )+TOLIN.LT.EPS*( STMP( I )-TOL ) ) THEN
253 VMAX = ONE / EPS
254 ELSE IF( SIN( I )+TOLIN.LT.STMP( I )-TOL ) THEN
255 VMAX = ( STMP( I )-TOL ) / ( SIN( I )+TOLIN )
256 ELSE
257 VMAX = ONE
258 END IF
259 IF( VMAX.GT.RMAX( 2 ) ) THEN
260 RMAX( 2 ) = VMAX
261 IF( NINFO( 2 ).EQ.0 )
262 $ LMAX( 2 ) = KNT
263 END IF
264 120 CONTINUE
265 *
266 * Compare condition numbers for eigenvectors
267 * taking their condition numbers into account
268 *
269 DO 130 I = 1, N
270 IF( V.GT.SEPTMP( I )*STMP( I ) ) THEN
271 TOL = SEPTMP( I )
272 ELSE
273 TOL = V / STMP( I )
274 END IF
275 IF( V.GT.SEPIN( I )*SIN( I ) ) THEN
276 TOLIN = SEPIN( I )
277 ELSE
278 TOLIN = V / SIN( I )
279 END IF
280 TOL = MAX( TOL, SMLNUM / EPS )
281 TOLIN = MAX( TOLIN, SMLNUM / EPS )
282 IF( EPS*( SEPIN( I )-TOLIN ).GT.SEPTMP( I )+TOL ) THEN
283 VMAX = ONE / EPS
284 ELSE IF( SEPIN( I )-TOLIN.GT.SEPTMP( I )+TOL ) THEN
285 VMAX = ( SEPIN( I )-TOLIN ) / ( SEPTMP( I )+TOL )
286 ELSE IF( SEPIN( I )+TOLIN.LT.EPS*( SEPTMP( I )-TOL ) ) THEN
287 VMAX = ONE / EPS
288 ELSE IF( SEPIN( I )+TOLIN.LT.SEPTMP( I )-TOL ) THEN
289 VMAX = ( SEPTMP( I )-TOL ) / ( SEPIN( I )+TOLIN )
290 ELSE
291 VMAX = ONE
292 END IF
293 IF( VMAX.GT.RMAX( 2 ) ) THEN
294 RMAX( 2 ) = VMAX
295 IF( NINFO( 2 ).EQ.0 )
296 $ LMAX( 2 ) = KNT
297 END IF
298 130 CONTINUE
299 *
300 * Compare condition numbers for eigenvalues
301 * without taking their condition numbers into account
302 *
303 DO 140 I = 1, N
304 IF( SIN( I ).LE.REAL( 2*N )*EPS .AND. STMP( I ).LE.
305 $ REAL( 2*N )*EPS ) THEN
306 VMAX = ONE
307 ELSE IF( EPS*SIN( I ).GT.STMP( I ) ) THEN
308 VMAX = ONE / EPS
309 ELSE IF( SIN( I ).GT.STMP( I ) ) THEN
310 VMAX = SIN( I ) / STMP( I )
311 ELSE IF( SIN( I ).LT.EPS*STMP( I ) ) THEN
312 VMAX = ONE / EPS
313 ELSE IF( SIN( I ).LT.STMP( I ) ) THEN
314 VMAX = STMP( I ) / SIN( I )
315 ELSE
316 VMAX = ONE
317 END IF
318 IF( VMAX.GT.RMAX( 3 ) ) THEN
319 RMAX( 3 ) = VMAX
320 IF( NINFO( 3 ).EQ.0 )
321 $ LMAX( 3 ) = KNT
322 END IF
323 140 CONTINUE
324 *
325 * Compare condition numbers for eigenvectors
326 * without taking their condition numbers into account
327 *
328 DO 150 I = 1, N
329 IF( SEPIN( I ).LE.V .AND. SEPTMP( I ).LE.V ) THEN
330 VMAX = ONE
331 ELSE IF( EPS*SEPIN( I ).GT.SEPTMP( I ) ) THEN
332 VMAX = ONE / EPS
333 ELSE IF( SEPIN( I ).GT.SEPTMP( I ) ) THEN
334 VMAX = SEPIN( I ) / SEPTMP( I )
335 ELSE IF( SEPIN( I ).LT.EPS*SEPTMP( I ) ) THEN
336 VMAX = ONE / EPS
337 ELSE IF( SEPIN( I ).LT.SEPTMP( I ) ) THEN
338 VMAX = SEPTMP( I ) / SEPIN( I )
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 150 CONTINUE
348 *
349 * Compute eigenvalue condition numbers only and compare
350 *
351 VMAX = ZERO
352 DUM( 1 ) = -ONE
353 CALL SCOPY( N, DUM, 0, STMP, 1 )
354 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
355 CALL CTRSNA( 'E', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT,
356 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
357 IF( INFO.NE.0 ) THEN
358 LMAX( 3 ) = KNT
359 NINFO( 3 ) = NINFO( 3 ) + 1
360 GO TO 260
361 END IF
362 DO 160 I = 1, N
363 IF( STMP( I ).NE.S( I ) )
364 $ VMAX = ONE / EPS
365 IF( SEPTMP( I ).NE.DUM( 1 ) )
366 $ VMAX = ONE / EPS
367 160 CONTINUE
368 *
369 * Compute eigenvector condition numbers only and compare
370 *
371 CALL SCOPY( N, DUM, 0, STMP, 1 )
372 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
373 CALL CTRSNA( 'V', 'A', SELECT, N, T, LDT, LE, LDT, RE, LDT,
374 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
375 IF( INFO.NE.0 ) THEN
376 LMAX( 3 ) = KNT
377 NINFO( 3 ) = NINFO( 3 ) + 1
378 GO TO 260
379 END IF
380 DO 170 I = 1, N
381 IF( STMP( I ).NE.DUM( 1 ) )
382 $ VMAX = ONE / EPS
383 IF( SEPTMP( I ).NE.SEP( I ) )
384 $ VMAX = ONE / EPS
385 170 CONTINUE
386 *
387 * Compute all condition numbers using SELECT and compare
388 *
389 DO 180 I = 1, N
390 SELECT( I ) = .TRUE.
391 180 CONTINUE
392 CALL SCOPY( N, DUM, 0, STMP, 1 )
393 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
394 CALL CTRSNA( 'B', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
395 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
396 IF( INFO.NE.0 ) THEN
397 LMAX( 3 ) = KNT
398 NINFO( 3 ) = NINFO( 3 ) + 1
399 GO TO 260
400 END IF
401 DO 190 I = 1, N
402 IF( SEPTMP( I ).NE.SEP( I ) )
403 $ VMAX = ONE / EPS
404 IF( STMP( I ).NE.S( I ) )
405 $ VMAX = ONE / EPS
406 190 CONTINUE
407 *
408 * Compute eigenvalue condition numbers using SELECT and compare
409 *
410 CALL SCOPY( N, DUM, 0, STMP, 1 )
411 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
412 CALL CTRSNA( 'E', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
413 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
414 IF( INFO.NE.0 ) THEN
415 LMAX( 3 ) = KNT
416 NINFO( 3 ) = NINFO( 3 ) + 1
417 GO TO 260
418 END IF
419 DO 200 I = 1, N
420 IF( STMP( I ).NE.S( I ) )
421 $ VMAX = ONE / EPS
422 IF( SEPTMP( I ).NE.DUM( 1 ) )
423 $ VMAX = ONE / EPS
424 200 CONTINUE
425 *
426 * Compute eigenvector condition numbers using SELECT and compare
427 *
428 CALL SCOPY( N, DUM, 0, STMP, 1 )
429 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
430 CALL CTRSNA( 'V', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
431 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
432 IF( INFO.NE.0 ) THEN
433 LMAX( 3 ) = KNT
434 NINFO( 3 ) = NINFO( 3 ) + 1
435 GO TO 260
436 END IF
437 DO 210 I = 1, N
438 IF( STMP( I ).NE.DUM( 1 ) )
439 $ VMAX = ONE / EPS
440 IF( SEPTMP( I ).NE.SEP( I ) )
441 $ VMAX = ONE / EPS
442 210 CONTINUE
443 IF( VMAX.GT.RMAX( 1 ) ) THEN
444 RMAX( 1 ) = VMAX
445 IF( NINFO( 1 ).EQ.0 )
446 $ LMAX( 1 ) = KNT
447 END IF
448 *
449 * Select second and next to last eigenvalues
450 *
451 DO 220 I = 1, N
452 SELECT( I ) = .FALSE.
453 220 CONTINUE
454 ICMP = 0
455 IF( N.GT.1 ) THEN
456 ICMP = 1
457 LCMP( 1 ) = 2
458 SELECT( 2 ) = .TRUE.
459 CALL CCOPY( N, RE( 1, 2 ), 1, RE( 1, 1 ), 1 )
460 CALL CCOPY( N, LE( 1, 2 ), 1, LE( 1, 1 ), 1 )
461 END IF
462 IF( N.GT.3 ) THEN
463 ICMP = 2
464 LCMP( 2 ) = N - 1
465 SELECT( N-1 ) = .TRUE.
466 CALL CCOPY( N, RE( 1, N-1 ), 1, RE( 1, 2 ), 1 )
467 CALL CCOPY( N, LE( 1, N-1 ), 1, LE( 1, 2 ), 1 )
468 END IF
469 *
470 * Compute all selected condition numbers
471 *
472 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
473 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
474 CALL CTRSNA( 'B', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
475 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
476 IF( INFO.NE.0 ) THEN
477 LMAX( 3 ) = KNT
478 NINFO( 3 ) = NINFO( 3 ) + 1
479 GO TO 260
480 END IF
481 DO 230 I = 1, ICMP
482 J = LCMP( I )
483 IF( SEPTMP( I ).NE.SEP( J ) )
484 $ VMAX = ONE / EPS
485 IF( STMP( I ).NE.S( J ) )
486 $ VMAX = ONE / EPS
487 230 CONTINUE
488 *
489 * Compute selected eigenvalue condition numbers
490 *
491 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
492 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
493 CALL CTRSNA( 'E', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
494 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
495 IF( INFO.NE.0 ) THEN
496 LMAX( 3 ) = KNT
497 NINFO( 3 ) = NINFO( 3 ) + 1
498 GO TO 260
499 END IF
500 DO 240 I = 1, ICMP
501 J = LCMP( I )
502 IF( STMP( I ).NE.S( J ) )
503 $ VMAX = ONE / EPS
504 IF( SEPTMP( I ).NE.DUM( 1 ) )
505 $ VMAX = ONE / EPS
506 240 CONTINUE
507 *
508 * Compute selected eigenvector condition numbers
509 *
510 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
511 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
512 CALL CTRSNA( 'V', 'S', SELECT, N, T, LDT, LE, LDT, RE, LDT,
513 $ STMP, SEPTMP, N, M, WORK, N, RWORK, INFO )
514 IF( INFO.NE.0 ) THEN
515 LMAX( 3 ) = KNT
516 NINFO( 3 ) = NINFO( 3 ) + 1
517 GO TO 260
518 END IF
519 DO 250 I = 1, ICMP
520 J = LCMP( I )
521 IF( STMP( I ).NE.DUM( 1 ) )
522 $ VMAX = ONE / EPS
523 IF( SEPTMP( I ).NE.SEP( J ) )
524 $ VMAX = ONE / EPS
525 250 CONTINUE
526 IF( VMAX.GT.RMAX( 1 ) ) THEN
527 RMAX( 1 ) = VMAX
528 IF( NINFO( 1 ).EQ.0 )
529 $ LMAX( 1 ) = KNT
530 END IF
531 260 CONTINUE
532 GO TO 10
533 *
534 * End of CGET37
535 *
536 END