1 SUBROUTINE SGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
2 $ H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
3 $ LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
4 $ RESULT, WORK, LWORK, IWORK, BWORK, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 LOGICAL COMP
12 INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
13 REAL RCDEIN, RCDVIN, THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL BWORK( * )
17 INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
18 REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ),
19 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
20 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
21 $ WR( * ), WRT( * ), WRTMP( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * SGET24 checks the nonsymmetric eigenvalue (Schur form) problem
28 * expert driver SGEESX.
29 *
30 * If COMP = .FALSE., the first 13 of the following tests will be
31 * be performed on the input matrix A, and also tests 14 and 15
32 * if LWORK is sufficiently large.
33 * If COMP = .TRUE., all 17 test will be performed.
34 *
35 * (1) 0 if T is in Schur form, 1/ulp otherwise
36 * (no sorting of eigenvalues)
37 *
38 * (2) | A - VS T VS' | / ( n |A| ulp )
39 *
40 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
41 * form (no sorting of eigenvalues).
42 *
43 * (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
44 *
45 * (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T
46 * 1/ulp otherwise
47 * (no sorting of eigenvalues)
48 *
49 * (5) 0 if T(with VS) = T(without VS),
50 * 1/ulp otherwise
51 * (no sorting of eigenvalues)
52 *
53 * (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
54 * 1/ulp otherwise
55 * (no sorting of eigenvalues)
56 *
57 * (7) 0 if T is in Schur form, 1/ulp otherwise
58 * (with sorting of eigenvalues)
59 *
60 * (8) | A - VS T VS' | / ( n |A| ulp )
61 *
62 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
63 * form (with sorting of eigenvalues).
64 *
65 * (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
66 *
67 * (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T
68 * 1/ulp otherwise
69 * If workspace sufficient, also compare WR, WI with and
70 * without reciprocal condition numbers
71 * (with sorting of eigenvalues)
72 *
73 * (11) 0 if T(with VS) = T(without VS),
74 * 1/ulp otherwise
75 * If workspace sufficient, also compare T with and without
76 * reciprocal condition numbers
77 * (with sorting of eigenvalues)
78 *
79 * (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
80 * 1/ulp otherwise
81 * If workspace sufficient, also compare VS with and without
82 * reciprocal condition numbers
83 * (with sorting of eigenvalues)
84 *
85 * (13) if sorting worked and SDIM is the number of
86 * eigenvalues which were SELECTed
87 * If workspace sufficient, also compare SDIM with and
88 * without reciprocal condition numbers
89 *
90 * (14) if RCONDE the same no matter if VS and/or RCONDV computed
91 *
92 * (15) if RCONDV the same no matter if VS and/or RCONDE computed
93 *
94 * (16) |RCONDE - RCDEIN| / cond(RCONDE)
95 *
96 * RCONDE is the reciprocal average eigenvalue condition number
97 * computed by SGEESX and RCDEIN (the precomputed true value)
98 * is supplied as input. cond(RCONDE) is the condition number
99 * of RCONDE, and takes errors in computing RCONDE into account,
100 * so that the resulting quantity should be O(ULP). cond(RCONDE)
101 * is essentially given by norm(A)/RCONDV.
102 *
103 * (17) |RCONDV - RCDVIN| / cond(RCONDV)
104 *
105 * RCONDV is the reciprocal right invariant subspace condition
106 * number computed by SGEESX and RCDVIN (the precomputed true
107 * value) is supplied as input. cond(RCONDV) is the condition
108 * number of RCONDV, and takes errors in computing RCONDV into
109 * account, so that the resulting quantity should be O(ULP).
110 * cond(RCONDV) is essentially given by norm(A)/RCONDE.
111 *
112 * Arguments
113 * =========
114 *
115 * COMP (input) LOGICAL
116 * COMP describes which input tests to perform:
117 * = .FALSE. if the computed condition numbers are not to
118 * be tested against RCDVIN and RCDEIN
119 * = .TRUE. if they are to be compared
120 *
121 * JTYPE (input) INTEGER
122 * Type of input matrix. Used to label output if error occurs.
123 *
124 * ISEED (input) INTEGER array, dimension (4)
125 * If COMP = .FALSE., the random number generator seed
126 * used to produce matrix.
127 * If COMP = .TRUE., ISEED(1) = the number of the example.
128 * Used to label output if error occurs.
129 *
130 * THRESH (input) REAL
131 * A test will count as "failed" if the "error", computed as
132 * described above, exceeds THRESH. Note that the error
133 * is scaled to be O(1), so THRESH should be a reasonably
134 * small multiple of 1, e.g., 10 or 100. In particular,
135 * it should not depend on the precision (single vs. double)
136 * or the size of the matrix. It must be at least zero.
137 *
138 * NOUNIT (input) INTEGER
139 * The FORTRAN unit number for printing out error messages
140 * (e.g., if a routine returns INFO not equal to 0.)
141 *
142 * N (input) INTEGER
143 * The dimension of A. N must be at least 0.
144 *
145 * A (input/output) REAL array, dimension (LDA, N)
146 * Used to hold the matrix whose eigenvalues are to be
147 * computed.
148 *
149 * LDA (input) INTEGER
150 * The leading dimension of A, and H. LDA must be at
151 * least 1 and at least N.
152 *
153 * H (workspace) REAL array, dimension (LDA, N)
154 * Another copy of the test matrix A, modified by SGEESX.
155 *
156 * HT (workspace) REAL array, dimension (LDA, N)
157 * Yet another copy of the test matrix A, modified by SGEESX.
158 *
159 * WR (workspace) REAL array, dimension (N)
160 * WI (workspace) REAL array, dimension (N)
161 * The real and imaginary parts of the eigenvalues of A.
162 * On exit, WR + WI*i are the eigenvalues of the matrix in A.
163 *
164 * WRT (workspace) REAL array, dimension (N)
165 * WIT (workspace) REAL array, dimension (N)
166 * Like WR, WI, these arrays contain the eigenvalues of A,
167 * but those computed when SGEESX only computes a partial
168 * eigendecomposition, i.e. not Schur vectors
169 *
170 * WRTMP (workspace) REAL array, dimension (N)
171 * WITMP (workspace) REAL array, dimension (N)
172 * Like WR, WI, these arrays contain the eigenvalues of A,
173 * but sorted by increasing real part.
174 *
175 * VS (workspace) REAL array, dimension (LDVS, N)
176 * VS holds the computed Schur vectors.
177 *
178 * LDVS (input) INTEGER
179 * Leading dimension of VS. Must be at least max(1, N).
180 *
181 * VS1 (workspace) REAL array, dimension (LDVS, N)
182 * VS1 holds another copy of the computed Schur vectors.
183 *
184 * RCDEIN (input) REAL
185 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
186 * condition number for the average of selected eigenvalues.
187 *
188 * RCDVIN (input) REAL
189 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
190 * condition number for the selected right invariant subspace.
191 *
192 * NSLCT (input) INTEGER
193 * When COMP = .TRUE. the number of selected eigenvalues
194 * corresponding to the precomputed values RCDEIN and RCDVIN.
195 *
196 * ISLCT (input) INTEGER array, dimension (NSLCT)
197 * When COMP = .TRUE. ISLCT selects the eigenvalues of the
198 * input matrix corresponding to the precomputed values RCDEIN
199 * and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
200 * eigenvalue with the J-th largest real part is selected.
201 * Not referenced if COMP = .FALSE.
202 *
203 * RESULT (output) REAL array, dimension (17)
204 * The values computed by the 17 tests described above.
205 * The values are currently limited to 1/ulp, to avoid
206 * overflow.
207 *
208 * WORK (workspace) REAL array, dimension (LWORK)
209 *
210 * LWORK (input) INTEGER
211 * The number of entries in WORK to be passed to SGEESX. This
212 * must be at least 3*N, and N+N**2 if tests 14--16 are to
213 * be performed.
214 *
215 * IWORK (workspace) INTEGER array, dimension (N*N)
216 *
217 * BWORK (workspace) LOGICAL array, dimension (N)
218 *
219 * INFO (output) INTEGER
220 * If 0, successful exit.
221 * If <0, input parameter -INFO had an incorrect value.
222 * If >0, SGEESX returned an error code, the absolute
223 * value of which is returned.
224 *
225 * =====================================================================
226 *
227 * .. Parameters ..
228 REAL ZERO, ONE
229 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
230 REAL EPSIN
231 PARAMETER ( EPSIN = 5.9605E-8 )
232 * ..
233 * .. Local Scalars ..
234 CHARACTER SORT
235 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
236 $ RSUB, SDIM, SDIM1
237 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
238 $ SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN,
239 $ VRMIN, WNORM
240 * ..
241 * .. Local Arrays ..
242 INTEGER IPNT( 20 )
243 * ..
244 * .. Arrays in Common ..
245 LOGICAL SELVAL( 20 )
246 REAL SELWI( 20 ), SELWR( 20 )
247 * ..
248 * .. Scalars in Common ..
249 INTEGER SELDIM, SELOPT
250 * ..
251 * .. Common blocks ..
252 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
253 * ..
254 * .. External Functions ..
255 LOGICAL SSLECT
256 REAL SLAMCH, SLANGE
257 EXTERNAL SSLECT, SLAMCH, SLANGE
258 * ..
259 * .. External Subroutines ..
260 EXTERNAL SCOPY, SGEESX, SGEMM, SLACPY, SORT01, XERBLA
261 * ..
262 * .. Intrinsic Functions ..
263 INTRINSIC ABS, MAX, MIN, REAL, SIGN, SQRT
264 * ..
265 * .. Executable Statements ..
266 *
267 * Check for errors
268 *
269 INFO = 0
270 IF( THRESH.LT.ZERO ) THEN
271 INFO = -3
272 ELSE IF( NOUNIT.LE.0 ) THEN
273 INFO = -5
274 ELSE IF( N.LT.0 ) THEN
275 INFO = -6
276 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
277 INFO = -8
278 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN
279 INFO = -18
280 ELSE IF( LWORK.LT.3*N ) THEN
281 INFO = -26
282 END IF
283 *
284 IF( INFO.NE.0 ) THEN
285 CALL XERBLA( 'SGET24', -INFO )
286 RETURN
287 END IF
288 *
289 * Quick return if nothing to do
290 *
291 DO 10 I = 1, 17
292 RESULT( I ) = -ONE
293 10 CONTINUE
294 *
295 IF( N.EQ.0 )
296 $ RETURN
297 *
298 * Important constants
299 *
300 SMLNUM = SLAMCH( 'Safe minimum' )
301 ULP = SLAMCH( 'Precision' )
302 ULPINV = ONE / ULP
303 *
304 * Perform tests (1)-(13)
305 *
306 SELOPT = 0
307 LIWORK = N*N
308 DO 120 ISORT = 0, 1
309 IF( ISORT.EQ.0 ) THEN
310 SORT = 'N'
311 RSUB = 0
312 ELSE
313 SORT = 'S'
314 RSUB = 6
315 END IF
316 *
317 * Compute Schur form and Schur vectors, and test them
318 *
319 CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
320 CALL SGEESX( 'V', SORT, SSLECT, 'N', N, H, LDA, SDIM, WR, WI,
321 $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK,
322 $ LIWORK, BWORK, IINFO )
323 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
324 RESULT( 1+RSUB ) = ULPINV
325 IF( JTYPE.NE.22 ) THEN
326 WRITE( NOUNIT, FMT = 9998 )'SGEESX1', IINFO, N, JTYPE,
327 $ ISEED
328 ELSE
329 WRITE( NOUNIT, FMT = 9999 )'SGEESX1', IINFO, N,
330 $ ISEED( 1 )
331 END IF
332 INFO = ABS( IINFO )
333 RETURN
334 END IF
335 IF( ISORT.EQ.0 ) THEN
336 CALL SCOPY( N, WR, 1, WRTMP, 1 )
337 CALL SCOPY( N, WI, 1, WITMP, 1 )
338 END IF
339 *
340 * Do Test (1) or Test (7)
341 *
342 RESULT( 1+RSUB ) = ZERO
343 DO 30 J = 1, N - 2
344 DO 20 I = J + 2, N
345 IF( H( I, J ).NE.ZERO )
346 $ RESULT( 1+RSUB ) = ULPINV
347 20 CONTINUE
348 30 CONTINUE
349 DO 40 I = 1, N - 2
350 IF( H( I+1, I ).NE.ZERO .AND. H( I+2, I+1 ).NE.ZERO )
351 $ RESULT( 1+RSUB ) = ULPINV
352 40 CONTINUE
353 DO 50 I = 1, N - 1
354 IF( H( I+1, I ).NE.ZERO ) THEN
355 IF( H( I, I ).NE.H( I+1, I+1 ) .OR. H( I, I+1 ).EQ.
356 $ ZERO .OR. SIGN( ONE, H( I+1, I ) ).EQ.
357 $ SIGN( ONE, H( I, I+1 ) ) )RESULT( 1+RSUB ) = ULPINV
358 END IF
359 50 CONTINUE
360 *
361 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
362 *
363 * Copy A to VS1, used as workspace
364 *
365 CALL SLACPY( ' ', N, N, A, LDA, VS1, LDVS )
366 *
367 * Compute Q*H and store in HT.
368 *
369 CALL SGEMM( 'No transpose', 'No transpose', N, N, N, ONE, VS,
370 $ LDVS, H, LDA, ZERO, HT, LDA )
371 *
372 * Compute A - Q*H*Q'
373 *
374 CALL SGEMM( 'No transpose', 'Transpose', N, N, N, -ONE, HT,
375 $ LDA, VS, LDVS, ONE, VS1, LDVS )
376 *
377 ANORM = MAX( SLANGE( '1', N, N, A, LDA, WORK ), SMLNUM )
378 WNORM = SLANGE( '1', N, N, VS1, LDVS, WORK )
379 *
380 IF( ANORM.GT.WNORM ) THEN
381 RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP )
382 ELSE
383 IF( ANORM.LT.ONE ) THEN
384 RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) /
385 $ ( N*ULP )
386 ELSE
387 RESULT( 2+RSUB ) = MIN( WNORM / ANORM, REAL( N ) ) /
388 $ ( N*ULP )
389 END IF
390 END IF
391 *
392 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
393 *
394 CALL SORT01( 'Columns', N, N, VS, LDVS, WORK, LWORK,
395 $ RESULT( 3+RSUB ) )
396 *
397 * Do Test (4) or Test (10)
398 *
399 RESULT( 4+RSUB ) = ZERO
400 DO 60 I = 1, N
401 IF( H( I, I ).NE.WR( I ) )
402 $ RESULT( 4+RSUB ) = ULPINV
403 60 CONTINUE
404 IF( N.GT.1 ) THEN
405 IF( H( 2, 1 ).EQ.ZERO .AND. WI( 1 ).NE.ZERO )
406 $ RESULT( 4+RSUB ) = ULPINV
407 IF( H( N, N-1 ).EQ.ZERO .AND. WI( N ).NE.ZERO )
408 $ RESULT( 4+RSUB ) = ULPINV
409 END IF
410 DO 70 I = 1, N - 1
411 IF( H( I+1, I ).NE.ZERO ) THEN
412 TMP = SQRT( ABS( H( I+1, I ) ) )*
413 $ SQRT( ABS( H( I, I+1 ) ) )
414 RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
415 $ ABS( WI( I )-TMP ) /
416 $ MAX( ULP*TMP, SMLNUM ) )
417 RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
418 $ ABS( WI( I+1 )+TMP ) /
419 $ MAX( ULP*TMP, SMLNUM ) )
420 ELSE IF( I.GT.1 ) THEN
421 IF( H( I+1, I ).EQ.ZERO .AND. H( I, I-1 ).EQ.ZERO .AND.
422 $ WI( I ).NE.ZERO )RESULT( 4+RSUB ) = ULPINV
423 END IF
424 70 CONTINUE
425 *
426 * Do Test (5) or Test (11)
427 *
428 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
429 CALL SGEESX( 'N', SORT, SSLECT, 'N', N, HT, LDA, SDIM, WRT,
430 $ WIT, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK,
431 $ LIWORK, BWORK, IINFO )
432 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
433 RESULT( 5+RSUB ) = ULPINV
434 IF( JTYPE.NE.22 ) THEN
435 WRITE( NOUNIT, FMT = 9998 )'SGEESX2', IINFO, N, JTYPE,
436 $ ISEED
437 ELSE
438 WRITE( NOUNIT, FMT = 9999 )'SGEESX2', IINFO, N,
439 $ ISEED( 1 )
440 END IF
441 INFO = ABS( IINFO )
442 GO TO 250
443 END IF
444 *
445 RESULT( 5+RSUB ) = ZERO
446 DO 90 J = 1, N
447 DO 80 I = 1, N
448 IF( H( I, J ).NE.HT( I, J ) )
449 $ RESULT( 5+RSUB ) = ULPINV
450 80 CONTINUE
451 90 CONTINUE
452 *
453 * Do Test (6) or Test (12)
454 *
455 RESULT( 6+RSUB ) = ZERO
456 DO 100 I = 1, N
457 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
458 $ RESULT( 6+RSUB ) = ULPINV
459 100 CONTINUE
460 *
461 * Do Test (13)
462 *
463 IF( ISORT.EQ.1 ) THEN
464 RESULT( 13 ) = ZERO
465 KNTEIG = 0
466 DO 110 I = 1, N
467 IF( SSLECT( WR( I ), WI( I ) ) .OR.
468 $ SSLECT( WR( I ), -WI( I ) ) )KNTEIG = KNTEIG + 1
469 IF( I.LT.N ) THEN
470 IF( ( SSLECT( WR( I+1 ), WI( I+1 ) ) .OR.
471 $ SSLECT( WR( I+1 ), -WI( I+1 ) ) ) .AND.
472 $ ( .NOT.( SSLECT( WR( I ),
473 $ WI( I ) ) .OR. SSLECT( WR( I ),
474 $ -WI( I ) ) ) ) .AND. IINFO.NE.N+2 )RESULT( 13 )
475 $ = ULPINV
476 END IF
477 110 CONTINUE
478 IF( SDIM.NE.KNTEIG )
479 $ RESULT( 13 ) = ULPINV
480 END IF
481 *
482 120 CONTINUE
483 *
484 * If there is enough workspace, perform tests (14) and (15)
485 * as well as (10) through (13)
486 *
487 IF( LWORK.GE.N+( N*N ) / 2 ) THEN
488 *
489 * Compute both RCONDE and RCONDV with VS
490 *
491 SORT = 'S'
492 RESULT( 14 ) = ZERO
493 RESULT( 15 ) = ZERO
494 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
495 CALL SGEESX( 'V', SORT, SSLECT, 'B', N, HT, LDA, SDIM1, WRT,
496 $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK,
497 $ IWORK, LIWORK, BWORK, IINFO )
498 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
499 RESULT( 14 ) = ULPINV
500 RESULT( 15 ) = ULPINV
501 IF( JTYPE.NE.22 ) THEN
502 WRITE( NOUNIT, FMT = 9998 )'SGEESX3', IINFO, N, JTYPE,
503 $ ISEED
504 ELSE
505 WRITE( NOUNIT, FMT = 9999 )'SGEESX3', IINFO, N,
506 $ ISEED( 1 )
507 END IF
508 INFO = ABS( IINFO )
509 GO TO 250
510 END IF
511 *
512 * Perform tests (10), (11), (12), and (13)
513 *
514 DO 140 I = 1, N
515 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
516 $ RESULT( 10 ) = ULPINV
517 DO 130 J = 1, N
518 IF( H( I, J ).NE.HT( I, J ) )
519 $ RESULT( 11 ) = ULPINV
520 IF( VS( I, J ).NE.VS1( I, J ) )
521 $ RESULT( 12 ) = ULPINV
522 130 CONTINUE
523 140 CONTINUE
524 IF( SDIM.NE.SDIM1 )
525 $ RESULT( 13 ) = ULPINV
526 *
527 * Compute both RCONDE and RCONDV without VS, and compare
528 *
529 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
530 CALL SGEESX( 'N', SORT, SSLECT, 'B', N, HT, LDA, SDIM1, WRT,
531 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
532 $ IWORK, LIWORK, BWORK, IINFO )
533 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
534 RESULT( 14 ) = ULPINV
535 RESULT( 15 ) = ULPINV
536 IF( JTYPE.NE.22 ) THEN
537 WRITE( NOUNIT, FMT = 9998 )'SGEESX4', IINFO, N, JTYPE,
538 $ ISEED
539 ELSE
540 WRITE( NOUNIT, FMT = 9999 )'SGEESX4', IINFO, N,
541 $ ISEED( 1 )
542 END IF
543 INFO = ABS( IINFO )
544 GO TO 250
545 END IF
546 *
547 * Perform tests (14) and (15)
548 *
549 IF( RCNDE1.NE.RCONDE )
550 $ RESULT( 14 ) = ULPINV
551 IF( RCNDV1.NE.RCONDV )
552 $ RESULT( 15 ) = ULPINV
553 *
554 * Perform tests (10), (11), (12), and (13)
555 *
556 DO 160 I = 1, N
557 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
558 $ RESULT( 10 ) = ULPINV
559 DO 150 J = 1, N
560 IF( H( I, J ).NE.HT( I, J ) )
561 $ RESULT( 11 ) = ULPINV
562 IF( VS( I, J ).NE.VS1( I, J ) )
563 $ RESULT( 12 ) = ULPINV
564 150 CONTINUE
565 160 CONTINUE
566 IF( SDIM.NE.SDIM1 )
567 $ RESULT( 13 ) = ULPINV
568 *
569 * Compute RCONDE with VS, and compare
570 *
571 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
572 CALL SGEESX( 'V', SORT, SSLECT, 'E', N, HT, LDA, SDIM1, WRT,
573 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
574 $ IWORK, LIWORK, BWORK, IINFO )
575 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
576 RESULT( 14 ) = ULPINV
577 IF( JTYPE.NE.22 ) THEN
578 WRITE( NOUNIT, FMT = 9998 )'SGEESX5', IINFO, N, JTYPE,
579 $ ISEED
580 ELSE
581 WRITE( NOUNIT, FMT = 9999 )'SGEESX5', IINFO, N,
582 $ ISEED( 1 )
583 END IF
584 INFO = ABS( IINFO )
585 GO TO 250
586 END IF
587 *
588 * Perform test (14)
589 *
590 IF( RCNDE1.NE.RCONDE )
591 $ RESULT( 14 ) = ULPINV
592 *
593 * Perform tests (10), (11), (12), and (13)
594 *
595 DO 180 I = 1, N
596 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
597 $ RESULT( 10 ) = ULPINV
598 DO 170 J = 1, N
599 IF( H( I, J ).NE.HT( I, J ) )
600 $ RESULT( 11 ) = ULPINV
601 IF( VS( I, J ).NE.VS1( I, J ) )
602 $ RESULT( 12 ) = ULPINV
603 170 CONTINUE
604 180 CONTINUE
605 IF( SDIM.NE.SDIM1 )
606 $ RESULT( 13 ) = ULPINV
607 *
608 * Compute RCONDE without VS, and compare
609 *
610 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
611 CALL SGEESX( 'N', SORT, SSLECT, 'E', N, HT, LDA, SDIM1, WRT,
612 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
613 $ IWORK, LIWORK, BWORK, IINFO )
614 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
615 RESULT( 14 ) = ULPINV
616 IF( JTYPE.NE.22 ) THEN
617 WRITE( NOUNIT, FMT = 9998 )'SGEESX6', IINFO, N, JTYPE,
618 $ ISEED
619 ELSE
620 WRITE( NOUNIT, FMT = 9999 )'SGEESX6', IINFO, N,
621 $ ISEED( 1 )
622 END IF
623 INFO = ABS( IINFO )
624 GO TO 250
625 END IF
626 *
627 * Perform test (14)
628 *
629 IF( RCNDE1.NE.RCONDE )
630 $ RESULT( 14 ) = ULPINV
631 *
632 * Perform tests (10), (11), (12), and (13)
633 *
634 DO 200 I = 1, N
635 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
636 $ RESULT( 10 ) = ULPINV
637 DO 190 J = 1, N
638 IF( H( I, J ).NE.HT( I, J ) )
639 $ RESULT( 11 ) = ULPINV
640 IF( VS( I, J ).NE.VS1( I, J ) )
641 $ RESULT( 12 ) = ULPINV
642 190 CONTINUE
643 200 CONTINUE
644 IF( SDIM.NE.SDIM1 )
645 $ RESULT( 13 ) = ULPINV
646 *
647 * Compute RCONDV with VS, and compare
648 *
649 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
650 CALL SGEESX( 'V', SORT, SSLECT, 'V', N, HT, LDA, SDIM1, WRT,
651 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
652 $ IWORK, LIWORK, BWORK, IINFO )
653 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
654 RESULT( 15 ) = ULPINV
655 IF( JTYPE.NE.22 ) THEN
656 WRITE( NOUNIT, FMT = 9998 )'SGEESX7', IINFO, N, JTYPE,
657 $ ISEED
658 ELSE
659 WRITE( NOUNIT, FMT = 9999 )'SGEESX7', IINFO, N,
660 $ ISEED( 1 )
661 END IF
662 INFO = ABS( IINFO )
663 GO TO 250
664 END IF
665 *
666 * Perform test (15)
667 *
668 IF( RCNDV1.NE.RCONDV )
669 $ RESULT( 15 ) = ULPINV
670 *
671 * Perform tests (10), (11), (12), and (13)
672 *
673 DO 220 I = 1, N
674 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
675 $ RESULT( 10 ) = ULPINV
676 DO 210 J = 1, N
677 IF( H( I, J ).NE.HT( I, J ) )
678 $ RESULT( 11 ) = ULPINV
679 IF( VS( I, J ).NE.VS1( I, J ) )
680 $ RESULT( 12 ) = ULPINV
681 210 CONTINUE
682 220 CONTINUE
683 IF( SDIM.NE.SDIM1 )
684 $ RESULT( 13 ) = ULPINV
685 *
686 * Compute RCONDV without VS, and compare
687 *
688 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
689 CALL SGEESX( 'N', SORT, SSLECT, 'V', N, HT, LDA, SDIM1, WRT,
690 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
691 $ IWORK, LIWORK, BWORK, IINFO )
692 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
693 RESULT( 15 ) = ULPINV
694 IF( JTYPE.NE.22 ) THEN
695 WRITE( NOUNIT, FMT = 9998 )'SGEESX8', IINFO, N, JTYPE,
696 $ ISEED
697 ELSE
698 WRITE( NOUNIT, FMT = 9999 )'SGEESX8', IINFO, N,
699 $ ISEED( 1 )
700 END IF
701 INFO = ABS( IINFO )
702 GO TO 250
703 END IF
704 *
705 * Perform test (15)
706 *
707 IF( RCNDV1.NE.RCONDV )
708 $ RESULT( 15 ) = ULPINV
709 *
710 * Perform tests (10), (11), (12), and (13)
711 *
712 DO 240 I = 1, N
713 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
714 $ RESULT( 10 ) = ULPINV
715 DO 230 J = 1, N
716 IF( H( I, J ).NE.HT( I, J ) )
717 $ RESULT( 11 ) = ULPINV
718 IF( VS( I, J ).NE.VS1( I, J ) )
719 $ RESULT( 12 ) = ULPINV
720 230 CONTINUE
721 240 CONTINUE
722 IF( SDIM.NE.SDIM1 )
723 $ RESULT( 13 ) = ULPINV
724 *
725 END IF
726 *
727 250 CONTINUE
728 *
729 * If there are precomputed reciprocal condition numbers, compare
730 * computed values with them.
731 *
732 IF( COMP ) THEN
733 *
734 * First set up SELOPT, SELDIM, SELVAL, SELWR, and SELWI so that
735 * the logical function SSLECT selects the eigenvalues specified
736 * by NSLCT and ISLCT.
737 *
738 SELDIM = N
739 SELOPT = 1
740 EPS = MAX( ULP, EPSIN )
741 DO 260 I = 1, N
742 IPNT( I ) = I
743 SELVAL( I ) = .FALSE.
744 SELWR( I ) = WRTMP( I )
745 SELWI( I ) = WITMP( I )
746 260 CONTINUE
747 DO 280 I = 1, N - 1
748 KMIN = I
749 VRMIN = WRTMP( I )
750 VIMIN = WITMP( I )
751 DO 270 J = I + 1, N
752 IF( WRTMP( J ).LT.VRMIN ) THEN
753 KMIN = J
754 VRMIN = WRTMP( J )
755 VIMIN = WITMP( J )
756 END IF
757 270 CONTINUE
758 WRTMP( KMIN ) = WRTMP( I )
759 WITMP( KMIN ) = WITMP( I )
760 WRTMP( I ) = VRMIN
761 WITMP( I ) = VIMIN
762 ITMP = IPNT( I )
763 IPNT( I ) = IPNT( KMIN )
764 IPNT( KMIN ) = ITMP
765 280 CONTINUE
766 DO 290 I = 1, NSLCT
767 SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE.
768 290 CONTINUE
769 *
770 * Compute condition numbers
771 *
772 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
773 CALL SGEESX( 'N', 'S', SSLECT, 'B', N, HT, LDA, SDIM1, WRT,
774 $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK,
775 $ IWORK, LIWORK, BWORK, IINFO )
776 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
777 RESULT( 16 ) = ULPINV
778 RESULT( 17 ) = ULPINV
779 WRITE( NOUNIT, FMT = 9999 )'SGEESX9', IINFO, N, ISEED( 1 )
780 INFO = ABS( IINFO )
781 GO TO 300
782 END IF
783 *
784 * Compare condition number for average of selected eigenvalues
785 * taking its condition number into account
786 *
787 ANORM = SLANGE( '1', N, N, A, LDA, WORK )
788 V = MAX( REAL( N )*EPS*ANORM, SMLNUM )
789 IF( ANORM.EQ.ZERO )
790 $ V = ONE
791 IF( V.GT.RCONDV ) THEN
792 TOL = ONE
793 ELSE
794 TOL = V / RCONDV
795 END IF
796 IF( V.GT.RCDVIN ) THEN
797 TOLIN = ONE
798 ELSE
799 TOLIN = V / RCDVIN
800 END IF
801 TOL = MAX( TOL, SMLNUM / EPS )
802 TOLIN = MAX( TOLIN, SMLNUM / EPS )
803 IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN
804 RESULT( 16 ) = ULPINV
805 ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN
806 RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
807 ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN
808 RESULT( 16 ) = ULPINV
809 ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN
810 RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN )
811 ELSE
812 RESULT( 16 ) = ONE
813 END IF
814 *
815 * Compare condition numbers for right invariant subspace
816 * taking its condition number into account
817 *
818 IF( V.GT.RCONDV*RCONDE ) THEN
819 TOL = RCONDV
820 ELSE
821 TOL = V / RCONDE
822 END IF
823 IF( V.GT.RCDVIN*RCDEIN ) THEN
824 TOLIN = RCDVIN
825 ELSE
826 TOLIN = V / RCDEIN
827 END IF
828 TOL = MAX( TOL, SMLNUM / EPS )
829 TOLIN = MAX( TOLIN, SMLNUM / EPS )
830 IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN
831 RESULT( 17 ) = ULPINV
832 ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN
833 RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
834 ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN
835 RESULT( 17 ) = ULPINV
836 ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN
837 RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN )
838 ELSE
839 RESULT( 17 ) = ONE
840 END IF
841 *
842 300 CONTINUE
843 *
844 END IF
845 *
846 9999 FORMAT( ' SGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
847 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
848 9998 FORMAT( ' SGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
849 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
850 *
851 RETURN
852 *
853 * End of SGET24
854 *
855 END
2 $ H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
3 $ LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
4 $ RESULT, WORK, LWORK, IWORK, BWORK, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 LOGICAL COMP
12 INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
13 REAL RCDEIN, RCDVIN, THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL BWORK( * )
17 INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
18 REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ),
19 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
20 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
21 $ WR( * ), WRT( * ), WRTMP( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * SGET24 checks the nonsymmetric eigenvalue (Schur form) problem
28 * expert driver SGEESX.
29 *
30 * If COMP = .FALSE., the first 13 of the following tests will be
31 * be performed on the input matrix A, and also tests 14 and 15
32 * if LWORK is sufficiently large.
33 * If COMP = .TRUE., all 17 test will be performed.
34 *
35 * (1) 0 if T is in Schur form, 1/ulp otherwise
36 * (no sorting of eigenvalues)
37 *
38 * (2) | A - VS T VS' | / ( n |A| ulp )
39 *
40 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
41 * form (no sorting of eigenvalues).
42 *
43 * (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
44 *
45 * (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T
46 * 1/ulp otherwise
47 * (no sorting of eigenvalues)
48 *
49 * (5) 0 if T(with VS) = T(without VS),
50 * 1/ulp otherwise
51 * (no sorting of eigenvalues)
52 *
53 * (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
54 * 1/ulp otherwise
55 * (no sorting of eigenvalues)
56 *
57 * (7) 0 if T is in Schur form, 1/ulp otherwise
58 * (with sorting of eigenvalues)
59 *
60 * (8) | A - VS T VS' | / ( n |A| ulp )
61 *
62 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
63 * form (with sorting of eigenvalues).
64 *
65 * (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
66 *
67 * (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T
68 * 1/ulp otherwise
69 * If workspace sufficient, also compare WR, WI with and
70 * without reciprocal condition numbers
71 * (with sorting of eigenvalues)
72 *
73 * (11) 0 if T(with VS) = T(without VS),
74 * 1/ulp otherwise
75 * If workspace sufficient, also compare T with and without
76 * reciprocal condition numbers
77 * (with sorting of eigenvalues)
78 *
79 * (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
80 * 1/ulp otherwise
81 * If workspace sufficient, also compare VS with and without
82 * reciprocal condition numbers
83 * (with sorting of eigenvalues)
84 *
85 * (13) if sorting worked and SDIM is the number of
86 * eigenvalues which were SELECTed
87 * If workspace sufficient, also compare SDIM with and
88 * without reciprocal condition numbers
89 *
90 * (14) if RCONDE the same no matter if VS and/or RCONDV computed
91 *
92 * (15) if RCONDV the same no matter if VS and/or RCONDE computed
93 *
94 * (16) |RCONDE - RCDEIN| / cond(RCONDE)
95 *
96 * RCONDE is the reciprocal average eigenvalue condition number
97 * computed by SGEESX and RCDEIN (the precomputed true value)
98 * is supplied as input. cond(RCONDE) is the condition number
99 * of RCONDE, and takes errors in computing RCONDE into account,
100 * so that the resulting quantity should be O(ULP). cond(RCONDE)
101 * is essentially given by norm(A)/RCONDV.
102 *
103 * (17) |RCONDV - RCDVIN| / cond(RCONDV)
104 *
105 * RCONDV is the reciprocal right invariant subspace condition
106 * number computed by SGEESX and RCDVIN (the precomputed true
107 * value) is supplied as input. cond(RCONDV) is the condition
108 * number of RCONDV, and takes errors in computing RCONDV into
109 * account, so that the resulting quantity should be O(ULP).
110 * cond(RCONDV) is essentially given by norm(A)/RCONDE.
111 *
112 * Arguments
113 * =========
114 *
115 * COMP (input) LOGICAL
116 * COMP describes which input tests to perform:
117 * = .FALSE. if the computed condition numbers are not to
118 * be tested against RCDVIN and RCDEIN
119 * = .TRUE. if they are to be compared
120 *
121 * JTYPE (input) INTEGER
122 * Type of input matrix. Used to label output if error occurs.
123 *
124 * ISEED (input) INTEGER array, dimension (4)
125 * If COMP = .FALSE., the random number generator seed
126 * used to produce matrix.
127 * If COMP = .TRUE., ISEED(1) = the number of the example.
128 * Used to label output if error occurs.
129 *
130 * THRESH (input) REAL
131 * A test will count as "failed" if the "error", computed as
132 * described above, exceeds THRESH. Note that the error
133 * is scaled to be O(1), so THRESH should be a reasonably
134 * small multiple of 1, e.g., 10 or 100. In particular,
135 * it should not depend on the precision (single vs. double)
136 * or the size of the matrix. It must be at least zero.
137 *
138 * NOUNIT (input) INTEGER
139 * The FORTRAN unit number for printing out error messages
140 * (e.g., if a routine returns INFO not equal to 0.)
141 *
142 * N (input) INTEGER
143 * The dimension of A. N must be at least 0.
144 *
145 * A (input/output) REAL array, dimension (LDA, N)
146 * Used to hold the matrix whose eigenvalues are to be
147 * computed.
148 *
149 * LDA (input) INTEGER
150 * The leading dimension of A, and H. LDA must be at
151 * least 1 and at least N.
152 *
153 * H (workspace) REAL array, dimension (LDA, N)
154 * Another copy of the test matrix A, modified by SGEESX.
155 *
156 * HT (workspace) REAL array, dimension (LDA, N)
157 * Yet another copy of the test matrix A, modified by SGEESX.
158 *
159 * WR (workspace) REAL array, dimension (N)
160 * WI (workspace) REAL array, dimension (N)
161 * The real and imaginary parts of the eigenvalues of A.
162 * On exit, WR + WI*i are the eigenvalues of the matrix in A.
163 *
164 * WRT (workspace) REAL array, dimension (N)
165 * WIT (workspace) REAL array, dimension (N)
166 * Like WR, WI, these arrays contain the eigenvalues of A,
167 * but those computed when SGEESX only computes a partial
168 * eigendecomposition, i.e. not Schur vectors
169 *
170 * WRTMP (workspace) REAL array, dimension (N)
171 * WITMP (workspace) REAL array, dimension (N)
172 * Like WR, WI, these arrays contain the eigenvalues of A,
173 * but sorted by increasing real part.
174 *
175 * VS (workspace) REAL array, dimension (LDVS, N)
176 * VS holds the computed Schur vectors.
177 *
178 * LDVS (input) INTEGER
179 * Leading dimension of VS. Must be at least max(1, N).
180 *
181 * VS1 (workspace) REAL array, dimension (LDVS, N)
182 * VS1 holds another copy of the computed Schur vectors.
183 *
184 * RCDEIN (input) REAL
185 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
186 * condition number for the average of selected eigenvalues.
187 *
188 * RCDVIN (input) REAL
189 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
190 * condition number for the selected right invariant subspace.
191 *
192 * NSLCT (input) INTEGER
193 * When COMP = .TRUE. the number of selected eigenvalues
194 * corresponding to the precomputed values RCDEIN and RCDVIN.
195 *
196 * ISLCT (input) INTEGER array, dimension (NSLCT)
197 * When COMP = .TRUE. ISLCT selects the eigenvalues of the
198 * input matrix corresponding to the precomputed values RCDEIN
199 * and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
200 * eigenvalue with the J-th largest real part is selected.
201 * Not referenced if COMP = .FALSE.
202 *
203 * RESULT (output) REAL array, dimension (17)
204 * The values computed by the 17 tests described above.
205 * The values are currently limited to 1/ulp, to avoid
206 * overflow.
207 *
208 * WORK (workspace) REAL array, dimension (LWORK)
209 *
210 * LWORK (input) INTEGER
211 * The number of entries in WORK to be passed to SGEESX. This
212 * must be at least 3*N, and N+N**2 if tests 14--16 are to
213 * be performed.
214 *
215 * IWORK (workspace) INTEGER array, dimension (N*N)
216 *
217 * BWORK (workspace) LOGICAL array, dimension (N)
218 *
219 * INFO (output) INTEGER
220 * If 0, successful exit.
221 * If <0, input parameter -INFO had an incorrect value.
222 * If >0, SGEESX returned an error code, the absolute
223 * value of which is returned.
224 *
225 * =====================================================================
226 *
227 * .. Parameters ..
228 REAL ZERO, ONE
229 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
230 REAL EPSIN
231 PARAMETER ( EPSIN = 5.9605E-8 )
232 * ..
233 * .. Local Scalars ..
234 CHARACTER SORT
235 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
236 $ RSUB, SDIM, SDIM1
237 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
238 $ SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN,
239 $ VRMIN, WNORM
240 * ..
241 * .. Local Arrays ..
242 INTEGER IPNT( 20 )
243 * ..
244 * .. Arrays in Common ..
245 LOGICAL SELVAL( 20 )
246 REAL SELWI( 20 ), SELWR( 20 )
247 * ..
248 * .. Scalars in Common ..
249 INTEGER SELDIM, SELOPT
250 * ..
251 * .. Common blocks ..
252 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
253 * ..
254 * .. External Functions ..
255 LOGICAL SSLECT
256 REAL SLAMCH, SLANGE
257 EXTERNAL SSLECT, SLAMCH, SLANGE
258 * ..
259 * .. External Subroutines ..
260 EXTERNAL SCOPY, SGEESX, SGEMM, SLACPY, SORT01, XERBLA
261 * ..
262 * .. Intrinsic Functions ..
263 INTRINSIC ABS, MAX, MIN, REAL, SIGN, SQRT
264 * ..
265 * .. Executable Statements ..
266 *
267 * Check for errors
268 *
269 INFO = 0
270 IF( THRESH.LT.ZERO ) THEN
271 INFO = -3
272 ELSE IF( NOUNIT.LE.0 ) THEN
273 INFO = -5
274 ELSE IF( N.LT.0 ) THEN
275 INFO = -6
276 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
277 INFO = -8
278 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN
279 INFO = -18
280 ELSE IF( LWORK.LT.3*N ) THEN
281 INFO = -26
282 END IF
283 *
284 IF( INFO.NE.0 ) THEN
285 CALL XERBLA( 'SGET24', -INFO )
286 RETURN
287 END IF
288 *
289 * Quick return if nothing to do
290 *
291 DO 10 I = 1, 17
292 RESULT( I ) = -ONE
293 10 CONTINUE
294 *
295 IF( N.EQ.0 )
296 $ RETURN
297 *
298 * Important constants
299 *
300 SMLNUM = SLAMCH( 'Safe minimum' )
301 ULP = SLAMCH( 'Precision' )
302 ULPINV = ONE / ULP
303 *
304 * Perform tests (1)-(13)
305 *
306 SELOPT = 0
307 LIWORK = N*N
308 DO 120 ISORT = 0, 1
309 IF( ISORT.EQ.0 ) THEN
310 SORT = 'N'
311 RSUB = 0
312 ELSE
313 SORT = 'S'
314 RSUB = 6
315 END IF
316 *
317 * Compute Schur form and Schur vectors, and test them
318 *
319 CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
320 CALL SGEESX( 'V', SORT, SSLECT, 'N', N, H, LDA, SDIM, WR, WI,
321 $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK,
322 $ LIWORK, BWORK, IINFO )
323 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
324 RESULT( 1+RSUB ) = ULPINV
325 IF( JTYPE.NE.22 ) THEN
326 WRITE( NOUNIT, FMT = 9998 )'SGEESX1', IINFO, N, JTYPE,
327 $ ISEED
328 ELSE
329 WRITE( NOUNIT, FMT = 9999 )'SGEESX1', IINFO, N,
330 $ ISEED( 1 )
331 END IF
332 INFO = ABS( IINFO )
333 RETURN
334 END IF
335 IF( ISORT.EQ.0 ) THEN
336 CALL SCOPY( N, WR, 1, WRTMP, 1 )
337 CALL SCOPY( N, WI, 1, WITMP, 1 )
338 END IF
339 *
340 * Do Test (1) or Test (7)
341 *
342 RESULT( 1+RSUB ) = ZERO
343 DO 30 J = 1, N - 2
344 DO 20 I = J + 2, N
345 IF( H( I, J ).NE.ZERO )
346 $ RESULT( 1+RSUB ) = ULPINV
347 20 CONTINUE
348 30 CONTINUE
349 DO 40 I = 1, N - 2
350 IF( H( I+1, I ).NE.ZERO .AND. H( I+2, I+1 ).NE.ZERO )
351 $ RESULT( 1+RSUB ) = ULPINV
352 40 CONTINUE
353 DO 50 I = 1, N - 1
354 IF( H( I+1, I ).NE.ZERO ) THEN
355 IF( H( I, I ).NE.H( I+1, I+1 ) .OR. H( I, I+1 ).EQ.
356 $ ZERO .OR. SIGN( ONE, H( I+1, I ) ).EQ.
357 $ SIGN( ONE, H( I, I+1 ) ) )RESULT( 1+RSUB ) = ULPINV
358 END IF
359 50 CONTINUE
360 *
361 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
362 *
363 * Copy A to VS1, used as workspace
364 *
365 CALL SLACPY( ' ', N, N, A, LDA, VS1, LDVS )
366 *
367 * Compute Q*H and store in HT.
368 *
369 CALL SGEMM( 'No transpose', 'No transpose', N, N, N, ONE, VS,
370 $ LDVS, H, LDA, ZERO, HT, LDA )
371 *
372 * Compute A - Q*H*Q'
373 *
374 CALL SGEMM( 'No transpose', 'Transpose', N, N, N, -ONE, HT,
375 $ LDA, VS, LDVS, ONE, VS1, LDVS )
376 *
377 ANORM = MAX( SLANGE( '1', N, N, A, LDA, WORK ), SMLNUM )
378 WNORM = SLANGE( '1', N, N, VS1, LDVS, WORK )
379 *
380 IF( ANORM.GT.WNORM ) THEN
381 RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP )
382 ELSE
383 IF( ANORM.LT.ONE ) THEN
384 RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) /
385 $ ( N*ULP )
386 ELSE
387 RESULT( 2+RSUB ) = MIN( WNORM / ANORM, REAL( N ) ) /
388 $ ( N*ULP )
389 END IF
390 END IF
391 *
392 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
393 *
394 CALL SORT01( 'Columns', N, N, VS, LDVS, WORK, LWORK,
395 $ RESULT( 3+RSUB ) )
396 *
397 * Do Test (4) or Test (10)
398 *
399 RESULT( 4+RSUB ) = ZERO
400 DO 60 I = 1, N
401 IF( H( I, I ).NE.WR( I ) )
402 $ RESULT( 4+RSUB ) = ULPINV
403 60 CONTINUE
404 IF( N.GT.1 ) THEN
405 IF( H( 2, 1 ).EQ.ZERO .AND. WI( 1 ).NE.ZERO )
406 $ RESULT( 4+RSUB ) = ULPINV
407 IF( H( N, N-1 ).EQ.ZERO .AND. WI( N ).NE.ZERO )
408 $ RESULT( 4+RSUB ) = ULPINV
409 END IF
410 DO 70 I = 1, N - 1
411 IF( H( I+1, I ).NE.ZERO ) THEN
412 TMP = SQRT( ABS( H( I+1, I ) ) )*
413 $ SQRT( ABS( H( I, I+1 ) ) )
414 RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
415 $ ABS( WI( I )-TMP ) /
416 $ MAX( ULP*TMP, SMLNUM ) )
417 RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
418 $ ABS( WI( I+1 )+TMP ) /
419 $ MAX( ULP*TMP, SMLNUM ) )
420 ELSE IF( I.GT.1 ) THEN
421 IF( H( I+1, I ).EQ.ZERO .AND. H( I, I-1 ).EQ.ZERO .AND.
422 $ WI( I ).NE.ZERO )RESULT( 4+RSUB ) = ULPINV
423 END IF
424 70 CONTINUE
425 *
426 * Do Test (5) or Test (11)
427 *
428 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
429 CALL SGEESX( 'N', SORT, SSLECT, 'N', N, HT, LDA, SDIM, WRT,
430 $ WIT, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK,
431 $ LIWORK, BWORK, IINFO )
432 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
433 RESULT( 5+RSUB ) = ULPINV
434 IF( JTYPE.NE.22 ) THEN
435 WRITE( NOUNIT, FMT = 9998 )'SGEESX2', IINFO, N, JTYPE,
436 $ ISEED
437 ELSE
438 WRITE( NOUNIT, FMT = 9999 )'SGEESX2', IINFO, N,
439 $ ISEED( 1 )
440 END IF
441 INFO = ABS( IINFO )
442 GO TO 250
443 END IF
444 *
445 RESULT( 5+RSUB ) = ZERO
446 DO 90 J = 1, N
447 DO 80 I = 1, N
448 IF( H( I, J ).NE.HT( I, J ) )
449 $ RESULT( 5+RSUB ) = ULPINV
450 80 CONTINUE
451 90 CONTINUE
452 *
453 * Do Test (6) or Test (12)
454 *
455 RESULT( 6+RSUB ) = ZERO
456 DO 100 I = 1, N
457 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
458 $ RESULT( 6+RSUB ) = ULPINV
459 100 CONTINUE
460 *
461 * Do Test (13)
462 *
463 IF( ISORT.EQ.1 ) THEN
464 RESULT( 13 ) = ZERO
465 KNTEIG = 0
466 DO 110 I = 1, N
467 IF( SSLECT( WR( I ), WI( I ) ) .OR.
468 $ SSLECT( WR( I ), -WI( I ) ) )KNTEIG = KNTEIG + 1
469 IF( I.LT.N ) THEN
470 IF( ( SSLECT( WR( I+1 ), WI( I+1 ) ) .OR.
471 $ SSLECT( WR( I+1 ), -WI( I+1 ) ) ) .AND.
472 $ ( .NOT.( SSLECT( WR( I ),
473 $ WI( I ) ) .OR. SSLECT( WR( I ),
474 $ -WI( I ) ) ) ) .AND. IINFO.NE.N+2 )RESULT( 13 )
475 $ = ULPINV
476 END IF
477 110 CONTINUE
478 IF( SDIM.NE.KNTEIG )
479 $ RESULT( 13 ) = ULPINV
480 END IF
481 *
482 120 CONTINUE
483 *
484 * If there is enough workspace, perform tests (14) and (15)
485 * as well as (10) through (13)
486 *
487 IF( LWORK.GE.N+( N*N ) / 2 ) THEN
488 *
489 * Compute both RCONDE and RCONDV with VS
490 *
491 SORT = 'S'
492 RESULT( 14 ) = ZERO
493 RESULT( 15 ) = ZERO
494 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
495 CALL SGEESX( 'V', SORT, SSLECT, 'B', N, HT, LDA, SDIM1, WRT,
496 $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK,
497 $ IWORK, LIWORK, BWORK, IINFO )
498 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
499 RESULT( 14 ) = ULPINV
500 RESULT( 15 ) = ULPINV
501 IF( JTYPE.NE.22 ) THEN
502 WRITE( NOUNIT, FMT = 9998 )'SGEESX3', IINFO, N, JTYPE,
503 $ ISEED
504 ELSE
505 WRITE( NOUNIT, FMT = 9999 )'SGEESX3', IINFO, N,
506 $ ISEED( 1 )
507 END IF
508 INFO = ABS( IINFO )
509 GO TO 250
510 END IF
511 *
512 * Perform tests (10), (11), (12), and (13)
513 *
514 DO 140 I = 1, N
515 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
516 $ RESULT( 10 ) = ULPINV
517 DO 130 J = 1, N
518 IF( H( I, J ).NE.HT( I, J ) )
519 $ RESULT( 11 ) = ULPINV
520 IF( VS( I, J ).NE.VS1( I, J ) )
521 $ RESULT( 12 ) = ULPINV
522 130 CONTINUE
523 140 CONTINUE
524 IF( SDIM.NE.SDIM1 )
525 $ RESULT( 13 ) = ULPINV
526 *
527 * Compute both RCONDE and RCONDV without VS, and compare
528 *
529 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
530 CALL SGEESX( 'N', SORT, SSLECT, 'B', N, HT, LDA, SDIM1, WRT,
531 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
532 $ IWORK, LIWORK, BWORK, IINFO )
533 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
534 RESULT( 14 ) = ULPINV
535 RESULT( 15 ) = ULPINV
536 IF( JTYPE.NE.22 ) THEN
537 WRITE( NOUNIT, FMT = 9998 )'SGEESX4', IINFO, N, JTYPE,
538 $ ISEED
539 ELSE
540 WRITE( NOUNIT, FMT = 9999 )'SGEESX4', IINFO, N,
541 $ ISEED( 1 )
542 END IF
543 INFO = ABS( IINFO )
544 GO TO 250
545 END IF
546 *
547 * Perform tests (14) and (15)
548 *
549 IF( RCNDE1.NE.RCONDE )
550 $ RESULT( 14 ) = ULPINV
551 IF( RCNDV1.NE.RCONDV )
552 $ RESULT( 15 ) = ULPINV
553 *
554 * Perform tests (10), (11), (12), and (13)
555 *
556 DO 160 I = 1, N
557 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
558 $ RESULT( 10 ) = ULPINV
559 DO 150 J = 1, N
560 IF( H( I, J ).NE.HT( I, J ) )
561 $ RESULT( 11 ) = ULPINV
562 IF( VS( I, J ).NE.VS1( I, J ) )
563 $ RESULT( 12 ) = ULPINV
564 150 CONTINUE
565 160 CONTINUE
566 IF( SDIM.NE.SDIM1 )
567 $ RESULT( 13 ) = ULPINV
568 *
569 * Compute RCONDE with VS, and compare
570 *
571 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
572 CALL SGEESX( 'V', SORT, SSLECT, 'E', N, HT, LDA, SDIM1, WRT,
573 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
574 $ IWORK, LIWORK, BWORK, IINFO )
575 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
576 RESULT( 14 ) = ULPINV
577 IF( JTYPE.NE.22 ) THEN
578 WRITE( NOUNIT, FMT = 9998 )'SGEESX5', IINFO, N, JTYPE,
579 $ ISEED
580 ELSE
581 WRITE( NOUNIT, FMT = 9999 )'SGEESX5', IINFO, N,
582 $ ISEED( 1 )
583 END IF
584 INFO = ABS( IINFO )
585 GO TO 250
586 END IF
587 *
588 * Perform test (14)
589 *
590 IF( RCNDE1.NE.RCONDE )
591 $ RESULT( 14 ) = ULPINV
592 *
593 * Perform tests (10), (11), (12), and (13)
594 *
595 DO 180 I = 1, N
596 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
597 $ RESULT( 10 ) = ULPINV
598 DO 170 J = 1, N
599 IF( H( I, J ).NE.HT( I, J ) )
600 $ RESULT( 11 ) = ULPINV
601 IF( VS( I, J ).NE.VS1( I, J ) )
602 $ RESULT( 12 ) = ULPINV
603 170 CONTINUE
604 180 CONTINUE
605 IF( SDIM.NE.SDIM1 )
606 $ RESULT( 13 ) = ULPINV
607 *
608 * Compute RCONDE without VS, and compare
609 *
610 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
611 CALL SGEESX( 'N', SORT, SSLECT, 'E', N, HT, LDA, SDIM1, WRT,
612 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
613 $ IWORK, LIWORK, BWORK, IINFO )
614 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
615 RESULT( 14 ) = ULPINV
616 IF( JTYPE.NE.22 ) THEN
617 WRITE( NOUNIT, FMT = 9998 )'SGEESX6', IINFO, N, JTYPE,
618 $ ISEED
619 ELSE
620 WRITE( NOUNIT, FMT = 9999 )'SGEESX6', IINFO, N,
621 $ ISEED( 1 )
622 END IF
623 INFO = ABS( IINFO )
624 GO TO 250
625 END IF
626 *
627 * Perform test (14)
628 *
629 IF( RCNDE1.NE.RCONDE )
630 $ RESULT( 14 ) = ULPINV
631 *
632 * Perform tests (10), (11), (12), and (13)
633 *
634 DO 200 I = 1, N
635 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
636 $ RESULT( 10 ) = ULPINV
637 DO 190 J = 1, N
638 IF( H( I, J ).NE.HT( I, J ) )
639 $ RESULT( 11 ) = ULPINV
640 IF( VS( I, J ).NE.VS1( I, J ) )
641 $ RESULT( 12 ) = ULPINV
642 190 CONTINUE
643 200 CONTINUE
644 IF( SDIM.NE.SDIM1 )
645 $ RESULT( 13 ) = ULPINV
646 *
647 * Compute RCONDV with VS, and compare
648 *
649 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
650 CALL SGEESX( 'V', SORT, SSLECT, 'V', N, HT, LDA, SDIM1, WRT,
651 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
652 $ IWORK, LIWORK, BWORK, IINFO )
653 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
654 RESULT( 15 ) = ULPINV
655 IF( JTYPE.NE.22 ) THEN
656 WRITE( NOUNIT, FMT = 9998 )'SGEESX7', IINFO, N, JTYPE,
657 $ ISEED
658 ELSE
659 WRITE( NOUNIT, FMT = 9999 )'SGEESX7', IINFO, N,
660 $ ISEED( 1 )
661 END IF
662 INFO = ABS( IINFO )
663 GO TO 250
664 END IF
665 *
666 * Perform test (15)
667 *
668 IF( RCNDV1.NE.RCONDV )
669 $ RESULT( 15 ) = ULPINV
670 *
671 * Perform tests (10), (11), (12), and (13)
672 *
673 DO 220 I = 1, N
674 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
675 $ RESULT( 10 ) = ULPINV
676 DO 210 J = 1, N
677 IF( H( I, J ).NE.HT( I, J ) )
678 $ RESULT( 11 ) = ULPINV
679 IF( VS( I, J ).NE.VS1( I, J ) )
680 $ RESULT( 12 ) = ULPINV
681 210 CONTINUE
682 220 CONTINUE
683 IF( SDIM.NE.SDIM1 )
684 $ RESULT( 13 ) = ULPINV
685 *
686 * Compute RCONDV without VS, and compare
687 *
688 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
689 CALL SGEESX( 'N', SORT, SSLECT, 'V', N, HT, LDA, SDIM1, WRT,
690 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
691 $ IWORK, LIWORK, BWORK, IINFO )
692 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
693 RESULT( 15 ) = ULPINV
694 IF( JTYPE.NE.22 ) THEN
695 WRITE( NOUNIT, FMT = 9998 )'SGEESX8', IINFO, N, JTYPE,
696 $ ISEED
697 ELSE
698 WRITE( NOUNIT, FMT = 9999 )'SGEESX8', IINFO, N,
699 $ ISEED( 1 )
700 END IF
701 INFO = ABS( IINFO )
702 GO TO 250
703 END IF
704 *
705 * Perform test (15)
706 *
707 IF( RCNDV1.NE.RCONDV )
708 $ RESULT( 15 ) = ULPINV
709 *
710 * Perform tests (10), (11), (12), and (13)
711 *
712 DO 240 I = 1, N
713 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
714 $ RESULT( 10 ) = ULPINV
715 DO 230 J = 1, N
716 IF( H( I, J ).NE.HT( I, J ) )
717 $ RESULT( 11 ) = ULPINV
718 IF( VS( I, J ).NE.VS1( I, J ) )
719 $ RESULT( 12 ) = ULPINV
720 230 CONTINUE
721 240 CONTINUE
722 IF( SDIM.NE.SDIM1 )
723 $ RESULT( 13 ) = ULPINV
724 *
725 END IF
726 *
727 250 CONTINUE
728 *
729 * If there are precomputed reciprocal condition numbers, compare
730 * computed values with them.
731 *
732 IF( COMP ) THEN
733 *
734 * First set up SELOPT, SELDIM, SELVAL, SELWR, and SELWI so that
735 * the logical function SSLECT selects the eigenvalues specified
736 * by NSLCT and ISLCT.
737 *
738 SELDIM = N
739 SELOPT = 1
740 EPS = MAX( ULP, EPSIN )
741 DO 260 I = 1, N
742 IPNT( I ) = I
743 SELVAL( I ) = .FALSE.
744 SELWR( I ) = WRTMP( I )
745 SELWI( I ) = WITMP( I )
746 260 CONTINUE
747 DO 280 I = 1, N - 1
748 KMIN = I
749 VRMIN = WRTMP( I )
750 VIMIN = WITMP( I )
751 DO 270 J = I + 1, N
752 IF( WRTMP( J ).LT.VRMIN ) THEN
753 KMIN = J
754 VRMIN = WRTMP( J )
755 VIMIN = WITMP( J )
756 END IF
757 270 CONTINUE
758 WRTMP( KMIN ) = WRTMP( I )
759 WITMP( KMIN ) = WITMP( I )
760 WRTMP( I ) = VRMIN
761 WITMP( I ) = VIMIN
762 ITMP = IPNT( I )
763 IPNT( I ) = IPNT( KMIN )
764 IPNT( KMIN ) = ITMP
765 280 CONTINUE
766 DO 290 I = 1, NSLCT
767 SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE.
768 290 CONTINUE
769 *
770 * Compute condition numbers
771 *
772 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
773 CALL SGEESX( 'N', 'S', SSLECT, 'B', N, HT, LDA, SDIM1, WRT,
774 $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK,
775 $ IWORK, LIWORK, BWORK, IINFO )
776 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
777 RESULT( 16 ) = ULPINV
778 RESULT( 17 ) = ULPINV
779 WRITE( NOUNIT, FMT = 9999 )'SGEESX9', IINFO, N, ISEED( 1 )
780 INFO = ABS( IINFO )
781 GO TO 300
782 END IF
783 *
784 * Compare condition number for average of selected eigenvalues
785 * taking its condition number into account
786 *
787 ANORM = SLANGE( '1', N, N, A, LDA, WORK )
788 V = MAX( REAL( N )*EPS*ANORM, SMLNUM )
789 IF( ANORM.EQ.ZERO )
790 $ V = ONE
791 IF( V.GT.RCONDV ) THEN
792 TOL = ONE
793 ELSE
794 TOL = V / RCONDV
795 END IF
796 IF( V.GT.RCDVIN ) THEN
797 TOLIN = ONE
798 ELSE
799 TOLIN = V / RCDVIN
800 END IF
801 TOL = MAX( TOL, SMLNUM / EPS )
802 TOLIN = MAX( TOLIN, SMLNUM / EPS )
803 IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN
804 RESULT( 16 ) = ULPINV
805 ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN
806 RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
807 ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN
808 RESULT( 16 ) = ULPINV
809 ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN
810 RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN )
811 ELSE
812 RESULT( 16 ) = ONE
813 END IF
814 *
815 * Compare condition numbers for right invariant subspace
816 * taking its condition number into account
817 *
818 IF( V.GT.RCONDV*RCONDE ) THEN
819 TOL = RCONDV
820 ELSE
821 TOL = V / RCONDE
822 END IF
823 IF( V.GT.RCDVIN*RCDEIN ) THEN
824 TOLIN = RCDVIN
825 ELSE
826 TOLIN = V / RCDEIN
827 END IF
828 TOL = MAX( TOL, SMLNUM / EPS )
829 TOLIN = MAX( TOLIN, SMLNUM / EPS )
830 IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN
831 RESULT( 17 ) = ULPINV
832 ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN
833 RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
834 ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN
835 RESULT( 17 ) = ULPINV
836 ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN
837 RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN )
838 ELSE
839 RESULT( 17 ) = ONE
840 END IF
841 *
842 300 CONTINUE
843 *
844 END IF
845 *
846 9999 FORMAT( ' SGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
847 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
848 9998 FORMAT( ' SGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
849 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
850 *
851 RETURN
852 *
853 * End of SGET24
854 *
855 END