1 SUBROUTINE SGET37( 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 * SGET37 tests STRSNA, 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 STRSNA
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 SGEHRD returns INFO nonzero on example i, LMAX(1)=i
39 * If SHSEQR returns INFO nonzero on example i, LMAX(2)=i
40 * If STRSNA returns INFO nonzero on example i, LMAX(3)=i
41 *
42 * NINFO (output) INTEGER array, dimension (3)
43 * NINFO(1) = No. of times SGEHRD returned INFO nonzero
44 * NINFO(2) = No. of times SHSEQR returned INFO nonzero
45 * NINFO(3) = No. of times STRSNA 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, IFND, INFO, ISCL, J, KMIN, M, N
65 REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
66 $ VIMIN, VMAX, VMUL, VRMIN
67 * ..
68 * .. Local Arrays ..
69 LOGICAL SELECT( LDT )
70 INTEGER IWORK( 2*LDT ), LCMP( 3 )
71 REAL DUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
72 $ S( LDT ), SEP( LDT ), SEPIN( LDT ),
73 $ SEPTMP( LDT ), SIN( LDT ), STMP( LDT ),
74 $ T( LDT, LDT ), TMP( LDT, LDT ), VAL( 3 ),
75 $ WI( LDT ), WIIN( LDT ), WITMP( LDT ),
76 $ WORK( LWORK ), WR( LDT ), WRIN( LDT ),
77 $ WRTMP( LDT )
78 * ..
79 * .. External Functions ..
80 REAL SLAMCH, SLANGE
81 EXTERNAL SLAMCH, SLANGE
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL SCOPY, SGEHRD, SHSEQR, SLABAD, SLACPY, SSCAL,
85 $ STREVC, STRSNA
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC 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 *
111 VAL( 1 ) = SQRT( SMLNUM )
112 VAL( 2 ) = ONE
113 VAL( 3 ) = SQRT( BIGNUM )
114 *
115 * Read input data until N=0. Assume input eigenvalues are sorted
116 * lexicographically (increasing by real part, then decreasing by
117 * imaginary part)
118 *
119 10 CONTINUE
120 READ( NIN, FMT = * )N
121 IF( N.EQ.0 )
122 $ RETURN
123 DO 20 I = 1, N
124 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
125 20 CONTINUE
126 DO 30 I = 1, N
127 READ( NIN, FMT = * )WRIN( I ), WIIN( I ), SIN( I ), SEPIN( I )
128 30 CONTINUE
129 TNRM = SLANGE( 'M', N, N, TMP, LDT, WORK )
130 *
131 * Begin test
132 *
133 DO 240 ISCL = 1, 3
134 *
135 * Scale input matrix
136 *
137 KNT = KNT + 1
138 CALL SLACPY( 'F', N, N, TMP, LDT, T, LDT )
139 VMUL = VAL( ISCL )
140 DO 40 I = 1, N
141 CALL SSCAL( N, VMUL, T( 1, I ), 1 )
142 40 CONTINUE
143 IF( TNRM.EQ.ZERO )
144 $ VMUL = ONE
145 *
146 * Compute eigenvalues and eigenvectors
147 *
148 CALL SGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
149 $ INFO )
150 IF( INFO.NE.0 ) THEN
151 LMAX( 1 ) = KNT
152 NINFO( 1 ) = NINFO( 1 ) + 1
153 GO TO 240
154 END IF
155 DO 60 J = 1, N - 2
156 DO 50 I = J + 2, N
157 T( I, J ) = ZERO
158 50 CONTINUE
159 60 CONTINUE
160 *
161 * Compute Schur form
162 *
163 CALL SHSEQR( 'S', 'N', N, 1, N, T, LDT, WR, WI, DUM, 1, WORK,
164 $ LWORK, INFO )
165 IF( INFO.NE.0 ) THEN
166 LMAX( 2 ) = KNT
167 NINFO( 2 ) = NINFO( 2 ) + 1
168 GO TO 240
169 END IF
170 *
171 * Compute eigenvectors
172 *
173 CALL STREVC( 'Both', 'All', SELECT, N, T, LDT, LE, LDT, RE,
174 $ LDT, N, M, WORK, INFO )
175 *
176 * Compute condition numbers
177 *
178 CALL STRSNA( 'Both', 'All', SELECT, N, T, LDT, LE, LDT, RE,
179 $ LDT, S, SEP, N, M, WORK, N, IWORK, INFO )
180 IF( INFO.NE.0 ) THEN
181 LMAX( 3 ) = KNT
182 NINFO( 3 ) = NINFO( 3 ) + 1
183 GO TO 240
184 END IF
185 *
186 * Sort eigenvalues and condition numbers lexicographically
187 * to compare with inputs
188 *
189 CALL SCOPY( N, WR, 1, WRTMP, 1 )
190 CALL SCOPY( N, WI, 1, WITMP, 1 )
191 CALL SCOPY( N, S, 1, STMP, 1 )
192 CALL SCOPY( N, SEP, 1, SEPTMP, 1 )
193 CALL SSCAL( N, ONE / VMUL, SEPTMP, 1 )
194 DO 80 I = 1, N - 1
195 KMIN = I
196 VRMIN = WRTMP( I )
197 VIMIN = WITMP( I )
198 DO 70 J = I + 1, N
199 IF( WRTMP( J ).LT.VRMIN ) THEN
200 KMIN = J
201 VRMIN = WRTMP( J )
202 VIMIN = WITMP( J )
203 END IF
204 70 CONTINUE
205 WRTMP( KMIN ) = WRTMP( I )
206 WITMP( KMIN ) = WITMP( I )
207 WRTMP( I ) = VRMIN
208 WITMP( I ) = VIMIN
209 VRMIN = STMP( KMIN )
210 STMP( KMIN ) = STMP( I )
211 STMP( I ) = VRMIN
212 VRMIN = SEPTMP( KMIN )
213 SEPTMP( KMIN ) = SEPTMP( I )
214 SEPTMP( I ) = VRMIN
215 80 CONTINUE
216 *
217 * Compare condition numbers for eigenvalues
218 * taking their condition numbers into account
219 *
220 V = MAX( TWO*REAL( N )*EPS*TNRM, SMLNUM )
221 IF( TNRM.EQ.ZERO )
222 $ V = ONE
223 DO 90 I = 1, N
224 IF( V.GT.SEPTMP( I ) ) THEN
225 TOL = ONE
226 ELSE
227 TOL = V / SEPTMP( I )
228 END IF
229 IF( V.GT.SEPIN( I ) ) THEN
230 TOLIN = ONE
231 ELSE
232 TOLIN = V / SEPIN( I )
233 END IF
234 TOL = MAX( TOL, SMLNUM / EPS )
235 TOLIN = MAX( TOLIN, SMLNUM / EPS )
236 IF( EPS*( SIN( I )-TOLIN ).GT.STMP( I )+TOL ) THEN
237 VMAX = ONE / EPS
238 ELSE IF( SIN( I )-TOLIN.GT.STMP( I )+TOL ) THEN
239 VMAX = ( SIN( I )-TOLIN ) / ( STMP( I )+TOL )
240 ELSE IF( SIN( I )+TOLIN.LT.EPS*( STMP( I )-TOL ) ) THEN
241 VMAX = ONE / EPS
242 ELSE IF( SIN( I )+TOLIN.LT.STMP( I )-TOL ) THEN
243 VMAX = ( STMP( I )-TOL ) / ( SIN( I )+TOLIN )
244 ELSE
245 VMAX = ONE
246 END IF
247 IF( VMAX.GT.RMAX( 2 ) ) THEN
248 RMAX( 2 ) = VMAX
249 IF( NINFO( 2 ).EQ.0 )
250 $ LMAX( 2 ) = KNT
251 END IF
252 90 CONTINUE
253 *
254 * Compare condition numbers for eigenvectors
255 * taking their condition numbers into account
256 *
257 DO 100 I = 1, N
258 IF( V.GT.SEPTMP( I )*STMP( I ) ) THEN
259 TOL = SEPTMP( I )
260 ELSE
261 TOL = V / STMP( I )
262 END IF
263 IF( V.GT.SEPIN( I )*SIN( I ) ) THEN
264 TOLIN = SEPIN( I )
265 ELSE
266 TOLIN = V / SIN( I )
267 END IF
268 TOL = MAX( TOL, SMLNUM / EPS )
269 TOLIN = MAX( TOLIN, SMLNUM / EPS )
270 IF( EPS*( SEPIN( I )-TOLIN ).GT.SEPTMP( I )+TOL ) THEN
271 VMAX = ONE / EPS
272 ELSE IF( SEPIN( I )-TOLIN.GT.SEPTMP( I )+TOL ) THEN
273 VMAX = ( SEPIN( I )-TOLIN ) / ( SEPTMP( I )+TOL )
274 ELSE IF( SEPIN( I )+TOLIN.LT.EPS*( SEPTMP( I )-TOL ) ) THEN
275 VMAX = ONE / EPS
276 ELSE IF( SEPIN( I )+TOLIN.LT.SEPTMP( I )-TOL ) THEN
277 VMAX = ( SEPTMP( I )-TOL ) / ( SEPIN( I )+TOLIN )
278 ELSE
279 VMAX = ONE
280 END IF
281 IF( VMAX.GT.RMAX( 2 ) ) THEN
282 RMAX( 2 ) = VMAX
283 IF( NINFO( 2 ).EQ.0 )
284 $ LMAX( 2 ) = KNT
285 END IF
286 100 CONTINUE
287 *
288 * Compare condition numbers for eigenvalues
289 * without taking their condition numbers into account
290 *
291 DO 110 I = 1, N
292 IF( SIN( I ).LE.REAL( 2*N )*EPS .AND. STMP( I ).LE.
293 $ REAL( 2*N )*EPS ) THEN
294 VMAX = ONE
295 ELSE IF( EPS*SIN( I ).GT.STMP( I ) ) THEN
296 VMAX = ONE / EPS
297 ELSE IF( SIN( I ).GT.STMP( I ) ) THEN
298 VMAX = SIN( I ) / STMP( I )
299 ELSE IF( SIN( I ).LT.EPS*STMP( I ) ) THEN
300 VMAX = ONE / EPS
301 ELSE IF( SIN( I ).LT.STMP( I ) ) THEN
302 VMAX = STMP( I ) / SIN( I )
303 ELSE
304 VMAX = ONE
305 END IF
306 IF( VMAX.GT.RMAX( 3 ) ) THEN
307 RMAX( 3 ) = VMAX
308 IF( NINFO( 3 ).EQ.0 )
309 $ LMAX( 3 ) = KNT
310 END IF
311 110 CONTINUE
312 *
313 * Compare condition numbers for eigenvectors
314 * without taking their condition numbers into account
315 *
316 DO 120 I = 1, N
317 IF( SEPIN( I ).LE.V .AND. SEPTMP( I ).LE.V ) THEN
318 VMAX = ONE
319 ELSE IF( EPS*SEPIN( I ).GT.SEPTMP( I ) ) THEN
320 VMAX = ONE / EPS
321 ELSE IF( SEPIN( I ).GT.SEPTMP( I ) ) THEN
322 VMAX = SEPIN( I ) / SEPTMP( I )
323 ELSE IF( SEPIN( I ).LT.EPS*SEPTMP( I ) ) THEN
324 VMAX = ONE / EPS
325 ELSE IF( SEPIN( I ).LT.SEPTMP( I ) ) THEN
326 VMAX = SEPTMP( I ) / SEPIN( I )
327 ELSE
328 VMAX = ONE
329 END IF
330 IF( VMAX.GT.RMAX( 3 ) ) THEN
331 RMAX( 3 ) = VMAX
332 IF( NINFO( 3 ).EQ.0 )
333 $ LMAX( 3 ) = KNT
334 END IF
335 120 CONTINUE
336 *
337 * Compute eigenvalue condition numbers only and compare
338 *
339 VMAX = ZERO
340 DUM( 1 ) = -ONE
341 CALL SCOPY( N, DUM, 0, STMP, 1 )
342 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
343 CALL STRSNA( 'Eigcond', 'All', SELECT, N, T, LDT, LE, LDT, RE,
344 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
345 IF( INFO.NE.0 ) THEN
346 LMAX( 3 ) = KNT
347 NINFO( 3 ) = NINFO( 3 ) + 1
348 GO TO 240
349 END IF
350 DO 130 I = 1, N
351 IF( STMP( I ).NE.S( I ) )
352 $ VMAX = ONE / EPS
353 IF( SEPTMP( I ).NE.DUM( 1 ) )
354 $ VMAX = ONE / EPS
355 130 CONTINUE
356 *
357 * Compute eigenvector condition numbers only and compare
358 *
359 CALL SCOPY( N, DUM, 0, STMP, 1 )
360 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
361 CALL STRSNA( 'Veccond', 'All', SELECT, N, T, LDT, LE, LDT, RE,
362 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
363 IF( INFO.NE.0 ) THEN
364 LMAX( 3 ) = KNT
365 NINFO( 3 ) = NINFO( 3 ) + 1
366 GO TO 240
367 END IF
368 DO 140 I = 1, N
369 IF( STMP( I ).NE.DUM( 1 ) )
370 $ VMAX = ONE / EPS
371 IF( SEPTMP( I ).NE.SEP( I ) )
372 $ VMAX = ONE / EPS
373 140 CONTINUE
374 *
375 * Compute all condition numbers using SELECT and compare
376 *
377 DO 150 I = 1, N
378 SELECT( I ) = .TRUE.
379 150 CONTINUE
380 CALL SCOPY( N, DUM, 0, STMP, 1 )
381 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
382 CALL STRSNA( 'Bothcond', 'Some', SELECT, N, T, LDT, LE, LDT,
383 $ RE, LDT, STMP, SEPTMP, N, M, WORK, N, IWORK,
384 $ INFO )
385 IF( INFO.NE.0 ) THEN
386 LMAX( 3 ) = KNT
387 NINFO( 3 ) = NINFO( 3 ) + 1
388 GO TO 240
389 END IF
390 DO 160 I = 1, N
391 IF( SEPTMP( I ).NE.SEP( I ) )
392 $ VMAX = ONE / EPS
393 IF( STMP( I ).NE.S( I ) )
394 $ VMAX = ONE / EPS
395 160 CONTINUE
396 *
397 * Compute eigenvalue condition numbers using SELECT and compare
398 *
399 CALL SCOPY( N, DUM, 0, STMP, 1 )
400 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
401 CALL STRSNA( 'Eigcond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
402 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
403 IF( INFO.NE.0 ) THEN
404 LMAX( 3 ) = KNT
405 NINFO( 3 ) = NINFO( 3 ) + 1
406 GO TO 240
407 END IF
408 DO 170 I = 1, N
409 IF( STMP( I ).NE.S( I ) )
410 $ VMAX = ONE / EPS
411 IF( SEPTMP( I ).NE.DUM( 1 ) )
412 $ VMAX = ONE / EPS
413 170 CONTINUE
414 *
415 * Compute eigenvector condition numbers using SELECT and compare
416 *
417 CALL SCOPY( N, DUM, 0, STMP, 1 )
418 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
419 CALL STRSNA( 'Veccond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
420 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
421 IF( INFO.NE.0 ) THEN
422 LMAX( 3 ) = KNT
423 NINFO( 3 ) = NINFO( 3 ) + 1
424 GO TO 240
425 END IF
426 DO 180 I = 1, N
427 IF( STMP( I ).NE.DUM( 1 ) )
428 $ VMAX = ONE / EPS
429 IF( SEPTMP( I ).NE.SEP( I ) )
430 $ VMAX = ONE / EPS
431 180 CONTINUE
432 IF( VMAX.GT.RMAX( 1 ) ) THEN
433 RMAX( 1 ) = VMAX
434 IF( NINFO( 1 ).EQ.0 )
435 $ LMAX( 1 ) = KNT
436 END IF
437 *
438 * Select first real and first complex eigenvalue
439 *
440 IF( WI( 1 ).EQ.ZERO ) THEN
441 LCMP( 1 ) = 1
442 IFND = 0
443 DO 190 I = 2, N
444 IF( IFND.EQ.1 .OR. WI( I ).EQ.ZERO ) THEN
445 SELECT( I ) = .FALSE.
446 ELSE
447 IFND = 1
448 LCMP( 2 ) = I
449 LCMP( 3 ) = I + 1
450 CALL SCOPY( N, RE( 1, I ), 1, RE( 1, 2 ), 1 )
451 CALL SCOPY( N, RE( 1, I+1 ), 1, RE( 1, 3 ), 1 )
452 CALL SCOPY( N, LE( 1, I ), 1, LE( 1, 2 ), 1 )
453 CALL SCOPY( N, LE( 1, I+1 ), 1, LE( 1, 3 ), 1 )
454 END IF
455 190 CONTINUE
456 IF( IFND.EQ.0 ) THEN
457 ICMP = 1
458 ELSE
459 ICMP = 3
460 END IF
461 ELSE
462 LCMP( 1 ) = 1
463 LCMP( 2 ) = 2
464 IFND = 0
465 DO 200 I = 3, N
466 IF( IFND.EQ.1 .OR. WI( I ).NE.ZERO ) THEN
467 SELECT( I ) = .FALSE.
468 ELSE
469 LCMP( 3 ) = I
470 IFND = 1
471 CALL SCOPY( N, RE( 1, I ), 1, RE( 1, 3 ), 1 )
472 CALL SCOPY( N, LE( 1, I ), 1, LE( 1, 3 ), 1 )
473 END IF
474 200 CONTINUE
475 IF( IFND.EQ.0 ) THEN
476 ICMP = 2
477 ELSE
478 ICMP = 3
479 END IF
480 END IF
481 *
482 * Compute all selected condition numbers
483 *
484 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
485 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
486 CALL STRSNA( 'Bothcond', 'Some', SELECT, N, T, LDT, LE, LDT,
487 $ RE, LDT, STMP, SEPTMP, N, M, WORK, N, IWORK,
488 $ INFO )
489 IF( INFO.NE.0 ) THEN
490 LMAX( 3 ) = KNT
491 NINFO( 3 ) = NINFO( 3 ) + 1
492 GO TO 240
493 END IF
494 DO 210 I = 1, ICMP
495 J = LCMP( I )
496 IF( SEPTMP( I ).NE.SEP( J ) )
497 $ VMAX = ONE / EPS
498 IF( STMP( I ).NE.S( J ) )
499 $ VMAX = ONE / EPS
500 210 CONTINUE
501 *
502 * Compute selected eigenvalue condition numbers
503 *
504 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
505 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
506 CALL STRSNA( 'Eigcond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
507 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
508 IF( INFO.NE.0 ) THEN
509 LMAX( 3 ) = KNT
510 NINFO( 3 ) = NINFO( 3 ) + 1
511 GO TO 240
512 END IF
513 DO 220 I = 1, ICMP
514 J = LCMP( I )
515 IF( STMP( I ).NE.S( J ) )
516 $ VMAX = ONE / EPS
517 IF( SEPTMP( I ).NE.DUM( 1 ) )
518 $ VMAX = ONE / EPS
519 220 CONTINUE
520 *
521 * Compute selected eigenvector condition numbers
522 *
523 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
524 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
525 CALL STRSNA( 'Veccond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
526 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
527 IF( INFO.NE.0 ) THEN
528 LMAX( 3 ) = KNT
529 NINFO( 3 ) = NINFO( 3 ) + 1
530 GO TO 240
531 END IF
532 DO 230 I = 1, ICMP
533 J = LCMP( I )
534 IF( STMP( I ).NE.DUM( 1 ) )
535 $ VMAX = ONE / EPS
536 IF( SEPTMP( I ).NE.SEP( J ) )
537 $ VMAX = ONE / EPS
538 230 CONTINUE
539 IF( VMAX.GT.RMAX( 1 ) ) THEN
540 RMAX( 1 ) = VMAX
541 IF( NINFO( 1 ).EQ.0 )
542 $ LMAX( 1 ) = KNT
543 END IF
544 240 CONTINUE
545 GO TO 10
546 *
547 * End of SGET37
548 *
549 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 * SGET37 tests STRSNA, 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 STRSNA
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 SGEHRD returns INFO nonzero on example i, LMAX(1)=i
39 * If SHSEQR returns INFO nonzero on example i, LMAX(2)=i
40 * If STRSNA returns INFO nonzero on example i, LMAX(3)=i
41 *
42 * NINFO (output) INTEGER array, dimension (3)
43 * NINFO(1) = No. of times SGEHRD returned INFO nonzero
44 * NINFO(2) = No. of times SHSEQR returned INFO nonzero
45 * NINFO(3) = No. of times STRSNA 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, IFND, INFO, ISCL, J, KMIN, M, N
65 REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
66 $ VIMIN, VMAX, VMUL, VRMIN
67 * ..
68 * .. Local Arrays ..
69 LOGICAL SELECT( LDT )
70 INTEGER IWORK( 2*LDT ), LCMP( 3 )
71 REAL DUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
72 $ S( LDT ), SEP( LDT ), SEPIN( LDT ),
73 $ SEPTMP( LDT ), SIN( LDT ), STMP( LDT ),
74 $ T( LDT, LDT ), TMP( LDT, LDT ), VAL( 3 ),
75 $ WI( LDT ), WIIN( LDT ), WITMP( LDT ),
76 $ WORK( LWORK ), WR( LDT ), WRIN( LDT ),
77 $ WRTMP( LDT )
78 * ..
79 * .. External Functions ..
80 REAL SLAMCH, SLANGE
81 EXTERNAL SLAMCH, SLANGE
82 * ..
83 * .. External Subroutines ..
84 EXTERNAL SCOPY, SGEHRD, SHSEQR, SLABAD, SLACPY, SSCAL,
85 $ STREVC, STRSNA
86 * ..
87 * .. Intrinsic Functions ..
88 INTRINSIC 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 *
111 VAL( 1 ) = SQRT( SMLNUM )
112 VAL( 2 ) = ONE
113 VAL( 3 ) = SQRT( BIGNUM )
114 *
115 * Read input data until N=0. Assume input eigenvalues are sorted
116 * lexicographically (increasing by real part, then decreasing by
117 * imaginary part)
118 *
119 10 CONTINUE
120 READ( NIN, FMT = * )N
121 IF( N.EQ.0 )
122 $ RETURN
123 DO 20 I = 1, N
124 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
125 20 CONTINUE
126 DO 30 I = 1, N
127 READ( NIN, FMT = * )WRIN( I ), WIIN( I ), SIN( I ), SEPIN( I )
128 30 CONTINUE
129 TNRM = SLANGE( 'M', N, N, TMP, LDT, WORK )
130 *
131 * Begin test
132 *
133 DO 240 ISCL = 1, 3
134 *
135 * Scale input matrix
136 *
137 KNT = KNT + 1
138 CALL SLACPY( 'F', N, N, TMP, LDT, T, LDT )
139 VMUL = VAL( ISCL )
140 DO 40 I = 1, N
141 CALL SSCAL( N, VMUL, T( 1, I ), 1 )
142 40 CONTINUE
143 IF( TNRM.EQ.ZERO )
144 $ VMUL = ONE
145 *
146 * Compute eigenvalues and eigenvectors
147 *
148 CALL SGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N,
149 $ INFO )
150 IF( INFO.NE.0 ) THEN
151 LMAX( 1 ) = KNT
152 NINFO( 1 ) = NINFO( 1 ) + 1
153 GO TO 240
154 END IF
155 DO 60 J = 1, N - 2
156 DO 50 I = J + 2, N
157 T( I, J ) = ZERO
158 50 CONTINUE
159 60 CONTINUE
160 *
161 * Compute Schur form
162 *
163 CALL SHSEQR( 'S', 'N', N, 1, N, T, LDT, WR, WI, DUM, 1, WORK,
164 $ LWORK, INFO )
165 IF( INFO.NE.0 ) THEN
166 LMAX( 2 ) = KNT
167 NINFO( 2 ) = NINFO( 2 ) + 1
168 GO TO 240
169 END IF
170 *
171 * Compute eigenvectors
172 *
173 CALL STREVC( 'Both', 'All', SELECT, N, T, LDT, LE, LDT, RE,
174 $ LDT, N, M, WORK, INFO )
175 *
176 * Compute condition numbers
177 *
178 CALL STRSNA( 'Both', 'All', SELECT, N, T, LDT, LE, LDT, RE,
179 $ LDT, S, SEP, N, M, WORK, N, IWORK, INFO )
180 IF( INFO.NE.0 ) THEN
181 LMAX( 3 ) = KNT
182 NINFO( 3 ) = NINFO( 3 ) + 1
183 GO TO 240
184 END IF
185 *
186 * Sort eigenvalues and condition numbers lexicographically
187 * to compare with inputs
188 *
189 CALL SCOPY( N, WR, 1, WRTMP, 1 )
190 CALL SCOPY( N, WI, 1, WITMP, 1 )
191 CALL SCOPY( N, S, 1, STMP, 1 )
192 CALL SCOPY( N, SEP, 1, SEPTMP, 1 )
193 CALL SSCAL( N, ONE / VMUL, SEPTMP, 1 )
194 DO 80 I = 1, N - 1
195 KMIN = I
196 VRMIN = WRTMP( I )
197 VIMIN = WITMP( I )
198 DO 70 J = I + 1, N
199 IF( WRTMP( J ).LT.VRMIN ) THEN
200 KMIN = J
201 VRMIN = WRTMP( J )
202 VIMIN = WITMP( J )
203 END IF
204 70 CONTINUE
205 WRTMP( KMIN ) = WRTMP( I )
206 WITMP( KMIN ) = WITMP( I )
207 WRTMP( I ) = VRMIN
208 WITMP( I ) = VIMIN
209 VRMIN = STMP( KMIN )
210 STMP( KMIN ) = STMP( I )
211 STMP( I ) = VRMIN
212 VRMIN = SEPTMP( KMIN )
213 SEPTMP( KMIN ) = SEPTMP( I )
214 SEPTMP( I ) = VRMIN
215 80 CONTINUE
216 *
217 * Compare condition numbers for eigenvalues
218 * taking their condition numbers into account
219 *
220 V = MAX( TWO*REAL( N )*EPS*TNRM, SMLNUM )
221 IF( TNRM.EQ.ZERO )
222 $ V = ONE
223 DO 90 I = 1, N
224 IF( V.GT.SEPTMP( I ) ) THEN
225 TOL = ONE
226 ELSE
227 TOL = V / SEPTMP( I )
228 END IF
229 IF( V.GT.SEPIN( I ) ) THEN
230 TOLIN = ONE
231 ELSE
232 TOLIN = V / SEPIN( I )
233 END IF
234 TOL = MAX( TOL, SMLNUM / EPS )
235 TOLIN = MAX( TOLIN, SMLNUM / EPS )
236 IF( EPS*( SIN( I )-TOLIN ).GT.STMP( I )+TOL ) THEN
237 VMAX = ONE / EPS
238 ELSE IF( SIN( I )-TOLIN.GT.STMP( I )+TOL ) THEN
239 VMAX = ( SIN( I )-TOLIN ) / ( STMP( I )+TOL )
240 ELSE IF( SIN( I )+TOLIN.LT.EPS*( STMP( I )-TOL ) ) THEN
241 VMAX = ONE / EPS
242 ELSE IF( SIN( I )+TOLIN.LT.STMP( I )-TOL ) THEN
243 VMAX = ( STMP( I )-TOL ) / ( SIN( I )+TOLIN )
244 ELSE
245 VMAX = ONE
246 END IF
247 IF( VMAX.GT.RMAX( 2 ) ) THEN
248 RMAX( 2 ) = VMAX
249 IF( NINFO( 2 ).EQ.0 )
250 $ LMAX( 2 ) = KNT
251 END IF
252 90 CONTINUE
253 *
254 * Compare condition numbers for eigenvectors
255 * taking their condition numbers into account
256 *
257 DO 100 I = 1, N
258 IF( V.GT.SEPTMP( I )*STMP( I ) ) THEN
259 TOL = SEPTMP( I )
260 ELSE
261 TOL = V / STMP( I )
262 END IF
263 IF( V.GT.SEPIN( I )*SIN( I ) ) THEN
264 TOLIN = SEPIN( I )
265 ELSE
266 TOLIN = V / SIN( I )
267 END IF
268 TOL = MAX( TOL, SMLNUM / EPS )
269 TOLIN = MAX( TOLIN, SMLNUM / EPS )
270 IF( EPS*( SEPIN( I )-TOLIN ).GT.SEPTMP( I )+TOL ) THEN
271 VMAX = ONE / EPS
272 ELSE IF( SEPIN( I )-TOLIN.GT.SEPTMP( I )+TOL ) THEN
273 VMAX = ( SEPIN( I )-TOLIN ) / ( SEPTMP( I )+TOL )
274 ELSE IF( SEPIN( I )+TOLIN.LT.EPS*( SEPTMP( I )-TOL ) ) THEN
275 VMAX = ONE / EPS
276 ELSE IF( SEPIN( I )+TOLIN.LT.SEPTMP( I )-TOL ) THEN
277 VMAX = ( SEPTMP( I )-TOL ) / ( SEPIN( I )+TOLIN )
278 ELSE
279 VMAX = ONE
280 END IF
281 IF( VMAX.GT.RMAX( 2 ) ) THEN
282 RMAX( 2 ) = VMAX
283 IF( NINFO( 2 ).EQ.0 )
284 $ LMAX( 2 ) = KNT
285 END IF
286 100 CONTINUE
287 *
288 * Compare condition numbers for eigenvalues
289 * without taking their condition numbers into account
290 *
291 DO 110 I = 1, N
292 IF( SIN( I ).LE.REAL( 2*N )*EPS .AND. STMP( I ).LE.
293 $ REAL( 2*N )*EPS ) THEN
294 VMAX = ONE
295 ELSE IF( EPS*SIN( I ).GT.STMP( I ) ) THEN
296 VMAX = ONE / EPS
297 ELSE IF( SIN( I ).GT.STMP( I ) ) THEN
298 VMAX = SIN( I ) / STMP( I )
299 ELSE IF( SIN( I ).LT.EPS*STMP( I ) ) THEN
300 VMAX = ONE / EPS
301 ELSE IF( SIN( I ).LT.STMP( I ) ) THEN
302 VMAX = STMP( I ) / SIN( I )
303 ELSE
304 VMAX = ONE
305 END IF
306 IF( VMAX.GT.RMAX( 3 ) ) THEN
307 RMAX( 3 ) = VMAX
308 IF( NINFO( 3 ).EQ.0 )
309 $ LMAX( 3 ) = KNT
310 END IF
311 110 CONTINUE
312 *
313 * Compare condition numbers for eigenvectors
314 * without taking their condition numbers into account
315 *
316 DO 120 I = 1, N
317 IF( SEPIN( I ).LE.V .AND. SEPTMP( I ).LE.V ) THEN
318 VMAX = ONE
319 ELSE IF( EPS*SEPIN( I ).GT.SEPTMP( I ) ) THEN
320 VMAX = ONE / EPS
321 ELSE IF( SEPIN( I ).GT.SEPTMP( I ) ) THEN
322 VMAX = SEPIN( I ) / SEPTMP( I )
323 ELSE IF( SEPIN( I ).LT.EPS*SEPTMP( I ) ) THEN
324 VMAX = ONE / EPS
325 ELSE IF( SEPIN( I ).LT.SEPTMP( I ) ) THEN
326 VMAX = SEPTMP( I ) / SEPIN( I )
327 ELSE
328 VMAX = ONE
329 END IF
330 IF( VMAX.GT.RMAX( 3 ) ) THEN
331 RMAX( 3 ) = VMAX
332 IF( NINFO( 3 ).EQ.0 )
333 $ LMAX( 3 ) = KNT
334 END IF
335 120 CONTINUE
336 *
337 * Compute eigenvalue condition numbers only and compare
338 *
339 VMAX = ZERO
340 DUM( 1 ) = -ONE
341 CALL SCOPY( N, DUM, 0, STMP, 1 )
342 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
343 CALL STRSNA( 'Eigcond', 'All', SELECT, N, T, LDT, LE, LDT, RE,
344 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
345 IF( INFO.NE.0 ) THEN
346 LMAX( 3 ) = KNT
347 NINFO( 3 ) = NINFO( 3 ) + 1
348 GO TO 240
349 END IF
350 DO 130 I = 1, N
351 IF( STMP( I ).NE.S( I ) )
352 $ VMAX = ONE / EPS
353 IF( SEPTMP( I ).NE.DUM( 1 ) )
354 $ VMAX = ONE / EPS
355 130 CONTINUE
356 *
357 * Compute eigenvector condition numbers only and compare
358 *
359 CALL SCOPY( N, DUM, 0, STMP, 1 )
360 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
361 CALL STRSNA( 'Veccond', 'All', SELECT, N, T, LDT, LE, LDT, RE,
362 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
363 IF( INFO.NE.0 ) THEN
364 LMAX( 3 ) = KNT
365 NINFO( 3 ) = NINFO( 3 ) + 1
366 GO TO 240
367 END IF
368 DO 140 I = 1, N
369 IF( STMP( I ).NE.DUM( 1 ) )
370 $ VMAX = ONE / EPS
371 IF( SEPTMP( I ).NE.SEP( I ) )
372 $ VMAX = ONE / EPS
373 140 CONTINUE
374 *
375 * Compute all condition numbers using SELECT and compare
376 *
377 DO 150 I = 1, N
378 SELECT( I ) = .TRUE.
379 150 CONTINUE
380 CALL SCOPY( N, DUM, 0, STMP, 1 )
381 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
382 CALL STRSNA( 'Bothcond', 'Some', SELECT, N, T, LDT, LE, LDT,
383 $ RE, LDT, STMP, SEPTMP, N, M, WORK, N, IWORK,
384 $ INFO )
385 IF( INFO.NE.0 ) THEN
386 LMAX( 3 ) = KNT
387 NINFO( 3 ) = NINFO( 3 ) + 1
388 GO TO 240
389 END IF
390 DO 160 I = 1, N
391 IF( SEPTMP( I ).NE.SEP( I ) )
392 $ VMAX = ONE / EPS
393 IF( STMP( I ).NE.S( I ) )
394 $ VMAX = ONE / EPS
395 160 CONTINUE
396 *
397 * Compute eigenvalue condition numbers using SELECT and compare
398 *
399 CALL SCOPY( N, DUM, 0, STMP, 1 )
400 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
401 CALL STRSNA( 'Eigcond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
402 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
403 IF( INFO.NE.0 ) THEN
404 LMAX( 3 ) = KNT
405 NINFO( 3 ) = NINFO( 3 ) + 1
406 GO TO 240
407 END IF
408 DO 170 I = 1, N
409 IF( STMP( I ).NE.S( I ) )
410 $ VMAX = ONE / EPS
411 IF( SEPTMP( I ).NE.DUM( 1 ) )
412 $ VMAX = ONE / EPS
413 170 CONTINUE
414 *
415 * Compute eigenvector condition numbers using SELECT and compare
416 *
417 CALL SCOPY( N, DUM, 0, STMP, 1 )
418 CALL SCOPY( N, DUM, 0, SEPTMP, 1 )
419 CALL STRSNA( 'Veccond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
420 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
421 IF( INFO.NE.0 ) THEN
422 LMAX( 3 ) = KNT
423 NINFO( 3 ) = NINFO( 3 ) + 1
424 GO TO 240
425 END IF
426 DO 180 I = 1, N
427 IF( STMP( I ).NE.DUM( 1 ) )
428 $ VMAX = ONE / EPS
429 IF( SEPTMP( I ).NE.SEP( I ) )
430 $ VMAX = ONE / EPS
431 180 CONTINUE
432 IF( VMAX.GT.RMAX( 1 ) ) THEN
433 RMAX( 1 ) = VMAX
434 IF( NINFO( 1 ).EQ.0 )
435 $ LMAX( 1 ) = KNT
436 END IF
437 *
438 * Select first real and first complex eigenvalue
439 *
440 IF( WI( 1 ).EQ.ZERO ) THEN
441 LCMP( 1 ) = 1
442 IFND = 0
443 DO 190 I = 2, N
444 IF( IFND.EQ.1 .OR. WI( I ).EQ.ZERO ) THEN
445 SELECT( I ) = .FALSE.
446 ELSE
447 IFND = 1
448 LCMP( 2 ) = I
449 LCMP( 3 ) = I + 1
450 CALL SCOPY( N, RE( 1, I ), 1, RE( 1, 2 ), 1 )
451 CALL SCOPY( N, RE( 1, I+1 ), 1, RE( 1, 3 ), 1 )
452 CALL SCOPY( N, LE( 1, I ), 1, LE( 1, 2 ), 1 )
453 CALL SCOPY( N, LE( 1, I+1 ), 1, LE( 1, 3 ), 1 )
454 END IF
455 190 CONTINUE
456 IF( IFND.EQ.0 ) THEN
457 ICMP = 1
458 ELSE
459 ICMP = 3
460 END IF
461 ELSE
462 LCMP( 1 ) = 1
463 LCMP( 2 ) = 2
464 IFND = 0
465 DO 200 I = 3, N
466 IF( IFND.EQ.1 .OR. WI( I ).NE.ZERO ) THEN
467 SELECT( I ) = .FALSE.
468 ELSE
469 LCMP( 3 ) = I
470 IFND = 1
471 CALL SCOPY( N, RE( 1, I ), 1, RE( 1, 3 ), 1 )
472 CALL SCOPY( N, LE( 1, I ), 1, LE( 1, 3 ), 1 )
473 END IF
474 200 CONTINUE
475 IF( IFND.EQ.0 ) THEN
476 ICMP = 2
477 ELSE
478 ICMP = 3
479 END IF
480 END IF
481 *
482 * Compute all selected condition numbers
483 *
484 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
485 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
486 CALL STRSNA( 'Bothcond', 'Some', SELECT, N, T, LDT, LE, LDT,
487 $ RE, LDT, STMP, SEPTMP, N, M, WORK, N, IWORK,
488 $ INFO )
489 IF( INFO.NE.0 ) THEN
490 LMAX( 3 ) = KNT
491 NINFO( 3 ) = NINFO( 3 ) + 1
492 GO TO 240
493 END IF
494 DO 210 I = 1, ICMP
495 J = LCMP( I )
496 IF( SEPTMP( I ).NE.SEP( J ) )
497 $ VMAX = ONE / EPS
498 IF( STMP( I ).NE.S( J ) )
499 $ VMAX = ONE / EPS
500 210 CONTINUE
501 *
502 * Compute selected eigenvalue condition numbers
503 *
504 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
505 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
506 CALL STRSNA( 'Eigcond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
507 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
508 IF( INFO.NE.0 ) THEN
509 LMAX( 3 ) = KNT
510 NINFO( 3 ) = NINFO( 3 ) + 1
511 GO TO 240
512 END IF
513 DO 220 I = 1, ICMP
514 J = LCMP( I )
515 IF( STMP( I ).NE.S( J ) )
516 $ VMAX = ONE / EPS
517 IF( SEPTMP( I ).NE.DUM( 1 ) )
518 $ VMAX = ONE / EPS
519 220 CONTINUE
520 *
521 * Compute selected eigenvector condition numbers
522 *
523 CALL SCOPY( ICMP, DUM, 0, STMP, 1 )
524 CALL SCOPY( ICMP, DUM, 0, SEPTMP, 1 )
525 CALL STRSNA( 'Veccond', 'Some', SELECT, N, T, LDT, LE, LDT, RE,
526 $ LDT, STMP, SEPTMP, N, M, WORK, N, IWORK, INFO )
527 IF( INFO.NE.0 ) THEN
528 LMAX( 3 ) = KNT
529 NINFO( 3 ) = NINFO( 3 ) + 1
530 GO TO 240
531 END IF
532 DO 230 I = 1, ICMP
533 J = LCMP( I )
534 IF( STMP( I ).NE.DUM( 1 ) )
535 $ VMAX = ONE / EPS
536 IF( SEPTMP( I ).NE.SEP( J ) )
537 $ VMAX = ONE / EPS
538 230 CONTINUE
539 IF( VMAX.GT.RMAX( 1 ) ) THEN
540 RMAX( 1 ) = VMAX
541 IF( NINFO( 1 ).EQ.0 )
542 $ LMAX( 1 ) = KNT
543 END IF
544 240 CONTINUE
545 GO TO 10
546 *
547 * End of SGET37
548 *
549 END