1 SUBROUTINE SDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT,
3 $ WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK,
4 $ 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 INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
12 $ NTYPES
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL BWORK( * ), DOTYPE( * )
17 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
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 * SDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
28 * expert driver SGEESX.
29 *
30 * SDRVSX uses both test matrices generated randomly depending on
31 * data supplied in the calling sequence, as well as on data
32 * read from an input file and including precomputed condition
33 * numbers to which it compares the ones it computes.
34 *
35 * When SDRVSX is called, a number of matrix "sizes" ("n's") and a
36 * number of matrix "types" are specified. For each size ("n")
37 * and each type of matrix, one matrix will be generated and used
38 * to test the nonsymmetric eigenroutines. For each matrix, 15
39 * tests will be performed:
40 *
41 * (1) 0 if T is in Schur form, 1/ulp otherwise
42 * (no sorting of eigenvalues)
43 *
44 * (2) | A - VS T VS' | / ( n |A| ulp )
45 *
46 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
47 * form (no sorting of eigenvalues).
48 *
49 * (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
50 *
51 * (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T
52 * 1/ulp otherwise
53 * (no sorting of eigenvalues)
54 *
55 * (5) 0 if T(with VS) = T(without VS),
56 * 1/ulp otherwise
57 * (no sorting of eigenvalues)
58 *
59 * (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
60 * 1/ulp otherwise
61 * (no sorting of eigenvalues)
62 *
63 * (7) 0 if T is in Schur form, 1/ulp otherwise
64 * (with sorting of eigenvalues)
65 *
66 * (8) | A - VS T VS' | / ( n |A| ulp )
67 *
68 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
69 * form (with sorting of eigenvalues).
70 *
71 * (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
72 *
73 * (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T
74 * 1/ulp otherwise
75 * If workspace sufficient, also compare WR, WI with and
76 * without reciprocal condition numbers
77 * (with sorting of eigenvalues)
78 *
79 * (11) 0 if T(with VS) = T(without VS),
80 * 1/ulp otherwise
81 * If workspace sufficient, also compare T with and without
82 * reciprocal condition numbers
83 * (with sorting of eigenvalues)
84 *
85 * (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
86 * 1/ulp otherwise
87 * If workspace sufficient, also compare VS with and without
88 * reciprocal condition numbers
89 * (with sorting of eigenvalues)
90 *
91 * (13) if sorting worked and SDIM is the number of
92 * eigenvalues which were SELECTed
93 * If workspace sufficient, also compare SDIM with and
94 * without reciprocal condition numbers
95 *
96 * (14) if RCONDE the same no matter if VS and/or RCONDV computed
97 *
98 * (15) if RCONDV the same no matter if VS and/or RCONDE computed
99 *
100 * The "sizes" are specified by an array NN(1:NSIZES); the value of
101 * each element NN(j) specifies one size.
102 * The "types" are specified by a logical array DOTYPE( 1:NTYPES );
103 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
104 * Currently, the list of possible types is:
105 *
106 * (1) The zero matrix.
107 * (2) The identity matrix.
108 * (3) A (transposed) Jordan block, with 1's on the diagonal.
109 *
110 * (4) A diagonal matrix with evenly spaced entries
111 * 1, ..., ULP and random signs.
112 * (ULP = (first number larger than 1) - 1 )
113 * (5) A diagonal matrix with geometrically spaced entries
114 * 1, ..., ULP and random signs.
115 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
116 * and random signs.
117 *
118 * (7) Same as (4), but multiplied by a constant near
119 * the overflow threshold
120 * (8) Same as (4), but multiplied by a constant near
121 * the underflow threshold
122 *
123 * (9) A matrix of the form U' T U, where U is orthogonal and
124 * T has evenly spaced entries 1, ..., ULP with random signs
125 * on the diagonal and random O(1) entries in the upper
126 * triangle.
127 *
128 * (10) A matrix of the form U' T U, where U is orthogonal and
129 * T has geometrically spaced entries 1, ..., ULP with random
130 * signs on the diagonal and random O(1) entries in the upper
131 * triangle.
132 *
133 * (11) A matrix of the form U' T U, where U is orthogonal and
134 * T has "clustered" entries 1, ULP,..., ULP with random
135 * signs on the diagonal and random O(1) entries in the upper
136 * triangle.
137 *
138 * (12) A matrix of the form U' T U, where U is orthogonal and
139 * T has real or complex conjugate paired eigenvalues randomly
140 * chosen from ( ULP, 1 ) and random O(1) entries in the upper
141 * triangle.
142 *
143 * (13) A matrix of the form X' T X, where X has condition
144 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
145 * with random signs on the diagonal and random O(1) entries
146 * in the upper triangle.
147 *
148 * (14) A matrix of the form X' T X, where X has condition
149 * SQRT( ULP ) and T has geometrically spaced entries
150 * 1, ..., ULP with random signs on the diagonal and random
151 * O(1) entries in the upper triangle.
152 *
153 * (15) A matrix of the form X' T X, where X has condition
154 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
155 * with random signs on the diagonal and random O(1) entries
156 * in the upper triangle.
157 *
158 * (16) A matrix of the form X' T X, where X has condition
159 * SQRT( ULP ) and T has real or complex conjugate paired
160 * eigenvalues randomly chosen from ( ULP, 1 ) and random
161 * O(1) entries in the upper triangle.
162 *
163 * (17) Same as (16), but multiplied by a constant
164 * near the overflow threshold
165 * (18) Same as (16), but multiplied by a constant
166 * near the underflow threshold
167 *
168 * (19) Nonsymmetric matrix with random entries chosen from (-1,1).
169 * If N is at least 4, all entries in first two rows and last
170 * row, and first column and last two columns are zero.
171 * (20) Same as (19), but multiplied by a constant
172 * near the overflow threshold
173 * (21) Same as (19), but multiplied by a constant
174 * near the underflow threshold
175 *
176 * In addition, an input file will be read from logical unit number
177 * NIUNIT. The file contains matrices along with precomputed
178 * eigenvalues and reciprocal condition numbers for the eigenvalue
179 * average and right invariant subspace. For these matrices, in
180 * addition to tests (1) to (15) we will compute the following two
181 * tests:
182 *
183 * (16) |RCONDE - RCDEIN| / cond(RCONDE)
184 *
185 * RCONDE is the reciprocal average eigenvalue condition number
186 * computed by SGEESX and RCDEIN (the precomputed true value)
187 * is supplied as input. cond(RCONDE) is the condition number
188 * of RCONDE, and takes errors in computing RCONDE into account,
189 * so that the resulting quantity should be O(ULP). cond(RCONDE)
190 * is essentially given by norm(A)/RCONDV.
191 *
192 * (17) |RCONDV - RCDVIN| / cond(RCONDV)
193 *
194 * RCONDV is the reciprocal right invariant subspace condition
195 * number computed by SGEESX and RCDVIN (the precomputed true
196 * value) is supplied as input. cond(RCONDV) is the condition
197 * number of RCONDV, and takes errors in computing RCONDV into
198 * account, so that the resulting quantity should be O(ULP).
199 * cond(RCONDV) is essentially given by norm(A)/RCONDE.
200 *
201 * Arguments
202 * =========
203 *
204 * NSIZES (input) INTEGER
205 * The number of sizes of matrices to use. NSIZES must be at
206 * least zero. If it is zero, no randomly generated matrices
207 * are tested, but any test matrices read from NIUNIT will be
208 * tested.
209 *
210 * NN (input) INTEGER array, dimension (NSIZES)
211 * An array containing the sizes to be used for the matrices.
212 * Zero values will be skipped. The values must be at least
213 * zero.
214 *
215 * NTYPES (input) INTEGER
216 * The number of elements in DOTYPE. NTYPES must be at least
217 * zero. If it is zero, no randomly generated test matrices
218 * are tested, but and test matrices read from NIUNIT will be
219 * tested. If it is MAXTYP+1 and NSIZES is 1, then an
220 * additional type, MAXTYP+1 is defined, which is to use
221 * whatever matrix is in A. This is only useful if
222 * DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
223 *
224 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
225 * If DOTYPE(j) is .TRUE., then for each size in NN a
226 * matrix of that size and of type j will be generated.
227 * If NTYPES is smaller than the maximum number of types
228 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
229 * MAXTYP will not be generated. If NTYPES is larger
230 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
231 * will be ignored.
232 *
233 * ISEED (input/output) INTEGER array, dimension (4)
234 * On entry ISEED specifies the seed of the random number
235 * generator. The array elements should be between 0 and 4095;
236 * if not they will be reduced mod 4096. Also, ISEED(4) must
237 * be odd. The random number generator uses a linear
238 * congruential sequence limited to small integers, and so
239 * should produce machine independent random numbers. The
240 * values of ISEED are changed on exit, and can be used in the
241 * next call to SDRVSX to continue the same random number
242 * sequence.
243 *
244 * THRESH (input) REAL
245 * A test will count as "failed" if the "error", computed as
246 * described above, exceeds THRESH. Note that the error
247 * is scaled to be O(1), so THRESH should be a reasonably
248 * small multiple of 1, e.g., 10 or 100. In particular,
249 * it should not depend on the precision (single vs. double)
250 * or the size of the matrix. It must be at least zero.
251 *
252 * NIUNIT (input) INTEGER
253 * The FORTRAN unit number for reading in the data file of
254 * problems to solve.
255 *
256 * NOUNIT (input) INTEGER
257 * The FORTRAN unit number for printing out error messages
258 * (e.g., if a routine returns INFO not equal to 0.)
259 *
260 * A (workspace) REAL array, dimension (LDA, max(NN))
261 * Used to hold the matrix whose eigenvalues are to be
262 * computed. On exit, A contains the last matrix actually used.
263 *
264 * LDA (input) INTEGER
265 * The leading dimension of A, and H. LDA must be at
266 * least 1 and at least max( NN ).
267 *
268 * H (workspace) REAL array, dimension (LDA, max(NN))
269 * Another copy of the test matrix A, modified by SGEESX.
270 *
271 * HT (workspace) REAL array, dimension (LDA, max(NN))
272 * Yet another copy of the test matrix A, modified by SGEESX.
273 *
274 * WR (workspace) REAL array, dimension (max(NN))
275 * WI (workspace) REAL array, dimension (max(NN))
276 * The real and imaginary parts of the eigenvalues of A.
277 * On exit, WR + WI*i are the eigenvalues of the matrix in A.
278 *
279 * WRT (workspace) REAL array, dimension (max(NN))
280 * WIT (workspace) REAL array, dimension (max(NN))
281 * Like WR, WI, these arrays contain the eigenvalues of A,
282 * but those computed when SGEESX only computes a partial
283 * eigendecomposition, i.e. not Schur vectors
284 *
285 * WRTMP (workspace) REAL array, dimension (max(NN))
286 * WITMP (workspace) REAL array, dimension (max(NN))
287 * More temporary storage for eigenvalues.
288 *
289 * VS (workspace) REAL array, dimension (LDVS, max(NN))
290 * VS holds the computed Schur vectors.
291 *
292 * LDVS (input) INTEGER
293 * Leading dimension of VS. Must be at least max(1,max(NN)).
294 *
295 * VS1 (workspace) REAL array, dimension (LDVS, max(NN))
296 * VS1 holds another copy of the computed Schur vectors.
297 *
298 * RESULT (output) REAL array, dimension (17)
299 * The values computed by the 17 tests described above.
300 * The values are currently limited to 1/ulp, to avoid overflow.
301 *
302 * WORK (workspace) REAL array, dimension (LWORK)
303 *
304 * LWORK (input) INTEGER
305 * The number of entries in WORK. This must be at least
306 * max(3*NN(j),2*NN(j)**2) for all j.
307 *
308 * IWORK (workspace) INTEGER array, dimension (max(NN)*max(NN))
309 *
310 * INFO (output) INTEGER
311 * If 0, successful exit.
312 * <0, input parameter -INFO is incorrect
313 * >0, SLATMR, SLATMS, SLATME or SGET24 returned an error
314 * code and INFO is its absolute value
315 *
316 *-----------------------------------------------------------------------
317 *
318 * Some Local Variables and Parameters:
319 * ---- ----- --------- --- ----------
320 * ZERO, ONE Real 0 and 1.
321 * MAXTYP The number of types defined.
322 * NMAX Largest value in NN.
323 * NERRS The number of tests which have exceeded THRESH
324 * COND, CONDS,
325 * IMODE Values to be passed to the matrix generators.
326 * ANORM Norm of A; passed to matrix generators.
327 *
328 * OVFL, UNFL Overflow and underflow thresholds.
329 * ULP, ULPINV Finest relative precision and its inverse.
330 * RTULP, RTULPI Square roots of the previous 4 values.
331 * The following four arrays decode JTYPE:
332 * KTYPE(j) The general type (1-10) for type "j".
333 * KMODE(j) The MODE value to be passed to the matrix
334 * generator for type "j".
335 * KMAGN(j) The order of magnitude ( O(1),
336 * O(overflow^(1/2) ), O(underflow^(1/2) )
337 * KCONDS(j) Selectw whether CONDS is to be 1 or
338 * 1/sqrt(ulp). (0 means irrelevant.)
339 *
340 * =====================================================================
341 *
342 * .. Parameters ..
343 REAL ZERO, ONE
344 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
345 INTEGER MAXTYP
346 PARAMETER ( MAXTYP = 21 )
347 * ..
348 * .. Local Scalars ..
349 LOGICAL BADNN
350 CHARACTER*3 PATH
351 INTEGER I, IINFO, IMODE, ITYPE, IWK, J, JCOL, JSIZE,
352 $ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX,
353 $ NNWORK, NSLCT, NTEST, NTESTF, NTESTT
354 REAL ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
355 $ RTULP, RTULPI, ULP, ULPINV, UNFL
356 * ..
357 * .. Local Arrays ..
358 CHARACTER ADUMMA( 1 )
359 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
360 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ),
361 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
362 * ..
363 * .. Arrays in Common ..
364 LOGICAL SELVAL( 20 )
365 REAL SELWI( 20 ), SELWR( 20 )
366 * ..
367 * .. Scalars in Common ..
368 INTEGER SELDIM, SELOPT
369 * ..
370 * .. Common blocks ..
371 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
372 * ..
373 * .. External Functions ..
374 REAL SLAMCH
375 EXTERNAL SLAMCH
376 * ..
377 * .. External Subroutines ..
378 EXTERNAL SGET24, SLABAD, SLASUM, SLATME, SLATMR, SLATMS,
379 $ SLASET, XERBLA
380 * ..
381 * .. Intrinsic Functions ..
382 INTRINSIC ABS, MAX, MIN, SQRT
383 * ..
384 * .. Data statements ..
385 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
386 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
387 $ 3, 1, 2, 3 /
388 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
389 $ 1, 5, 5, 5, 4, 3, 1 /
390 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
391 * ..
392 * .. Executable Statements ..
393 *
394 PATH( 1: 1 ) = 'Single precision'
395 PATH( 2: 3 ) = 'SX'
396 *
397 * Check for errors
398 *
399 NTESTT = 0
400 NTESTF = 0
401 INFO = 0
402 *
403 * Important constants
404 *
405 BADNN = .FALSE.
406 *
407 * 12 is the largest dimension in the input file of precomputed
408 * problems
409 *
410 NMAX = 12
411 DO 10 J = 1, NSIZES
412 NMAX = MAX( NMAX, NN( J ) )
413 IF( NN( J ).LT.0 )
414 $ BADNN = .TRUE.
415 10 CONTINUE
416 *
417 * Check for errors
418 *
419 IF( NSIZES.LT.0 ) THEN
420 INFO = -1
421 ELSE IF( BADNN ) THEN
422 INFO = -2
423 ELSE IF( NTYPES.LT.0 ) THEN
424 INFO = -3
425 ELSE IF( THRESH.LT.ZERO ) THEN
426 INFO = -6
427 ELSE IF( NIUNIT.LE.0 ) THEN
428 INFO = -7
429 ELSE IF( NOUNIT.LE.0 ) THEN
430 INFO = -8
431 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
432 INFO = -10
433 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN
434 INFO = -20
435 ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN
436 INFO = -24
437 END IF
438 *
439 IF( INFO.NE.0 ) THEN
440 CALL XERBLA( 'SDRVSX', -INFO )
441 RETURN
442 END IF
443 *
444 * If nothing to do check on NIUNIT
445 *
446 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
447 $ GO TO 150
448 *
449 * More Important constants
450 *
451 UNFL = SLAMCH( 'Safe minimum' )
452 OVFL = ONE / UNFL
453 CALL SLABAD( UNFL, OVFL )
454 ULP = SLAMCH( 'Precision' )
455 ULPINV = ONE / ULP
456 RTULP = SQRT( ULP )
457 RTULPI = ONE / RTULP
458 *
459 * Loop over sizes, types
460 *
461 NERRS = 0
462 *
463 DO 140 JSIZE = 1, NSIZES
464 N = NN( JSIZE )
465 IF( NSIZES.NE.1 ) THEN
466 MTYPES = MIN( MAXTYP, NTYPES )
467 ELSE
468 MTYPES = MIN( MAXTYP+1, NTYPES )
469 END IF
470 *
471 DO 130 JTYPE = 1, MTYPES
472 IF( .NOT.DOTYPE( JTYPE ) )
473 $ GO TO 130
474 *
475 * Save ISEED in case of an error.
476 *
477 DO 20 J = 1, 4
478 IOLDSD( J ) = ISEED( J )
479 20 CONTINUE
480 *
481 * Compute "A"
482 *
483 * Control parameters:
484 *
485 * KMAGN KCONDS KMODE KTYPE
486 * =1 O(1) 1 clustered 1 zero
487 * =2 large large clustered 2 identity
488 * =3 small exponential Jordan
489 * =4 arithmetic diagonal, (w/ eigenvalues)
490 * =5 random log symmetric, w/ eigenvalues
491 * =6 random general, w/ eigenvalues
492 * =7 random diagonal
493 * =8 random symmetric
494 * =9 random general
495 * =10 random triangular
496 *
497 IF( MTYPES.GT.MAXTYP )
498 $ GO TO 90
499 *
500 ITYPE = KTYPE( JTYPE )
501 IMODE = KMODE( JTYPE )
502 *
503 * Compute norm
504 *
505 GO TO ( 30, 40, 50 )KMAGN( JTYPE )
506 *
507 30 CONTINUE
508 ANORM = ONE
509 GO TO 60
510 *
511 40 CONTINUE
512 ANORM = OVFL*ULP
513 GO TO 60
514 *
515 50 CONTINUE
516 ANORM = UNFL*ULPINV
517 GO TO 60
518 *
519 60 CONTINUE
520 *
521 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
522 IINFO = 0
523 COND = ULPINV
524 *
525 * Special Matrices -- Identity & Jordan block
526 *
527 * Zero
528 *
529 IF( ITYPE.EQ.1 ) THEN
530 IINFO = 0
531 *
532 ELSE IF( ITYPE.EQ.2 ) THEN
533 *
534 * Identity
535 *
536 DO 70 JCOL = 1, N
537 A( JCOL, JCOL ) = ANORM
538 70 CONTINUE
539 *
540 ELSE IF( ITYPE.EQ.3 ) THEN
541 *
542 * Jordan Block
543 *
544 DO 80 JCOL = 1, N
545 A( JCOL, JCOL ) = ANORM
546 IF( JCOL.GT.1 )
547 $ A( JCOL, JCOL-1 ) = ONE
548 80 CONTINUE
549 *
550 ELSE IF( ITYPE.EQ.4 ) THEN
551 *
552 * Diagonal Matrix, [Eigen]values Specified
553 *
554 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
555 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
556 $ IINFO )
557 *
558 ELSE IF( ITYPE.EQ.5 ) THEN
559 *
560 * Symmetric, eigenvalues specified
561 *
562 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
563 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
564 $ IINFO )
565 *
566 ELSE IF( ITYPE.EQ.6 ) THEN
567 *
568 * General, eigenvalues specified
569 *
570 IF( KCONDS( JTYPE ).EQ.1 ) THEN
571 CONDS = ONE
572 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
573 CONDS = RTULPI
574 ELSE
575 CONDS = ZERO
576 END IF
577 *
578 ADUMMA( 1 ) = ' '
579 CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE,
580 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4,
581 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ),
582 $ IINFO )
583 *
584 ELSE IF( ITYPE.EQ.7 ) THEN
585 *
586 * Diagonal, random eigenvalues
587 *
588 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
589 $ 'T', 'N', WORK( N+1 ), 1, ONE,
590 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
591 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
592 *
593 ELSE IF( ITYPE.EQ.8 ) THEN
594 *
595 * Symmetric, random eigenvalues
596 *
597 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
598 $ 'T', 'N', WORK( N+1 ), 1, ONE,
599 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
600 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
601 *
602 ELSE IF( ITYPE.EQ.9 ) THEN
603 *
604 * General, random eigenvalues
605 *
606 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
607 $ 'T', 'N', WORK( N+1 ), 1, ONE,
608 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
609 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
610 IF( N.GE.4 ) THEN
611 CALL SLASET( 'Full', 2, N, ZERO, ZERO, A, LDA )
612 CALL SLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ),
613 $ LDA )
614 CALL SLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ),
615 $ LDA )
616 CALL SLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ),
617 $ LDA )
618 END IF
619 *
620 ELSE IF( ITYPE.EQ.10 ) THEN
621 *
622 * Triangular, random eigenvalues
623 *
624 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
625 $ 'T', 'N', WORK( N+1 ), 1, ONE,
626 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
627 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
628 *
629 ELSE
630 *
631 IINFO = 1
632 END IF
633 *
634 IF( IINFO.NE.0 ) THEN
635 WRITE( NOUNIT, FMT = 9991 )'Generator', IINFO, N, JTYPE,
636 $ IOLDSD
637 INFO = ABS( IINFO )
638 RETURN
639 END IF
640 *
641 90 CONTINUE
642 *
643 * Test for minimal and generous workspace
644 *
645 DO 120 IWK = 1, 2
646 IF( IWK.EQ.1 ) THEN
647 NNWORK = 3*N
648 ELSE
649 NNWORK = MAX( 3*N, 2*N*N )
650 END IF
651 NNWORK = MAX( NNWORK, 1 )
652 *
653 CALL SGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N,
654 $ A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP,
655 $ WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT,
656 $ ISLCT, RESULT, WORK, NNWORK, IWORK, BWORK,
657 $ INFO )
658 *
659 * Check for RESULT(j) > THRESH
660 *
661 NTEST = 0
662 NFAIL = 0
663 DO 100 J = 1, 15
664 IF( RESULT( J ).GE.ZERO )
665 $ NTEST = NTEST + 1
666 IF( RESULT( J ).GE.THRESH )
667 $ NFAIL = NFAIL + 1
668 100 CONTINUE
669 *
670 IF( NFAIL.GT.0 )
671 $ NTESTF = NTESTF + 1
672 IF( NTESTF.EQ.1 ) THEN
673 WRITE( NOUNIT, FMT = 9999 )PATH
674 WRITE( NOUNIT, FMT = 9998 )
675 WRITE( NOUNIT, FMT = 9997 )
676 WRITE( NOUNIT, FMT = 9996 )
677 WRITE( NOUNIT, FMT = 9995 )THRESH
678 WRITE( NOUNIT, FMT = 9994 )
679 NTESTF = 2
680 END IF
681 *
682 DO 110 J = 1, 15
683 IF( RESULT( J ).GE.THRESH ) THEN
684 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
685 $ J, RESULT( J )
686 END IF
687 110 CONTINUE
688 *
689 NERRS = NERRS + NFAIL
690 NTESTT = NTESTT + NTEST
691 *
692 120 CONTINUE
693 130 CONTINUE
694 140 CONTINUE
695 *
696 150 CONTINUE
697 *
698 * Read in data from file to check accuracy of condition estimation
699 * Read input data until N=0
700 *
701 JTYPE = 0
702 160 CONTINUE
703 READ( NIUNIT, FMT = *, END = 200 )N, NSLCT
704 IF( N.EQ.0 )
705 $ GO TO 200
706 JTYPE = JTYPE + 1
707 ISEED( 1 ) = JTYPE
708 IF( NSLCT.GT.0 )
709 $ READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT )
710 DO 170 I = 1, N
711 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
712 170 CONTINUE
713 READ( NIUNIT, FMT = * )RCDEIN, RCDVIN
714 *
715 CALL SGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT,
716 $ WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1,
717 $ RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK,
718 $ IWORK, BWORK, INFO )
719 *
720 * Check for RESULT(j) > THRESH
721 *
722 NTEST = 0
723 NFAIL = 0
724 DO 180 J = 1, 17
725 IF( RESULT( J ).GE.ZERO )
726 $ NTEST = NTEST + 1
727 IF( RESULT( J ).GE.THRESH )
728 $ NFAIL = NFAIL + 1
729 180 CONTINUE
730 *
731 IF( NFAIL.GT.0 )
732 $ NTESTF = NTESTF + 1
733 IF( NTESTF.EQ.1 ) THEN
734 WRITE( NOUNIT, FMT = 9999 )PATH
735 WRITE( NOUNIT, FMT = 9998 )
736 WRITE( NOUNIT, FMT = 9997 )
737 WRITE( NOUNIT, FMT = 9996 )
738 WRITE( NOUNIT, FMT = 9995 )THRESH
739 WRITE( NOUNIT, FMT = 9994 )
740 NTESTF = 2
741 END IF
742 DO 190 J = 1, 17
743 IF( RESULT( J ).GE.THRESH ) THEN
744 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J )
745 END IF
746 190 CONTINUE
747 *
748 NERRS = NERRS + NFAIL
749 NTESTT = NTESTT + NTEST
750 GO TO 160
751 200 CONTINUE
752 *
753 * Summary
754 *
755 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
756 *
757 9999 FORMAT( / 1X, A3, ' -- Real Schur Form Decomposition Expert ',
758 $ 'Driver', / ' Matrix types (see SDRVSX for details):' )
759 *
760 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
761 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
762 $ / ' 2=Identity matrix. ', ' 6=Diagona',
763 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
764 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
765 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
766 $ 'mall, evenly spaced.' )
767 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
768 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
769 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
770 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
771 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
772 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
773 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
774 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
775 $ ' complx ' )
776 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
777 $ 'with small random entries.', / ' 20=Matrix with large ran',
778 $ 'dom entries. ', / )
779 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
780 $ / ' ( A denotes A on input and T denotes A on output)',
781 $ / / ' 1 = 0 if T in Schur form (no sort), ',
782 $ ' 1/ulp otherwise', /
783 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
784 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
785 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
786 $ ' 1/ulp otherwise', /
787 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
788 $ ' 1/ulp otherwise', /
789 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
790 $ ', 1/ulp otherwise' )
791 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
792 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
793 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
794 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
795 $ ' 1/ulp otherwise', /
796 $ ' 11 = 0 if T same no matter what else computed (sort),',
797 $ ' 1/ulp otherwise', /
798 $ ' 12 = 0 if WR, WI same no matter what else computed ',
799 $ '(sort), 1/ulp otherwise', /
800 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise',
801 $ / ' 14 = 0 if RCONDE same no matter what else computed,',
802 $ ' 1/ulp otherwise', /
803 $ ' 15 = 0 if RCONDv same no matter what else computed,',
804 $ ' 1/ulp otherwise', /
805 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
806 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
807 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
808 $ ' type ', I2, ', test(', I2, ')=', G10.3 )
809 9992 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=',
810 $ G10.3 )
811 9991 FORMAT( ' SDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
812 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
813 *
814 RETURN
815 *
816 * End of SDRVSX
817 *
818 END
2 $ NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT,
3 $ WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK,
4 $ 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 INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
12 $ NTYPES
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL BWORK( * ), DOTYPE( * )
17 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
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 * SDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
28 * expert driver SGEESX.
29 *
30 * SDRVSX uses both test matrices generated randomly depending on
31 * data supplied in the calling sequence, as well as on data
32 * read from an input file and including precomputed condition
33 * numbers to which it compares the ones it computes.
34 *
35 * When SDRVSX is called, a number of matrix "sizes" ("n's") and a
36 * number of matrix "types" are specified. For each size ("n")
37 * and each type of matrix, one matrix will be generated and used
38 * to test the nonsymmetric eigenroutines. For each matrix, 15
39 * tests will be performed:
40 *
41 * (1) 0 if T is in Schur form, 1/ulp otherwise
42 * (no sorting of eigenvalues)
43 *
44 * (2) | A - VS T VS' | / ( n |A| ulp )
45 *
46 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
47 * form (no sorting of eigenvalues).
48 *
49 * (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
50 *
51 * (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T
52 * 1/ulp otherwise
53 * (no sorting of eigenvalues)
54 *
55 * (5) 0 if T(with VS) = T(without VS),
56 * 1/ulp otherwise
57 * (no sorting of eigenvalues)
58 *
59 * (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
60 * 1/ulp otherwise
61 * (no sorting of eigenvalues)
62 *
63 * (7) 0 if T is in Schur form, 1/ulp otherwise
64 * (with sorting of eigenvalues)
65 *
66 * (8) | A - VS T VS' | / ( n |A| ulp )
67 *
68 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
69 * form (with sorting of eigenvalues).
70 *
71 * (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
72 *
73 * (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T
74 * 1/ulp otherwise
75 * If workspace sufficient, also compare WR, WI with and
76 * without reciprocal condition numbers
77 * (with sorting of eigenvalues)
78 *
79 * (11) 0 if T(with VS) = T(without VS),
80 * 1/ulp otherwise
81 * If workspace sufficient, also compare T with and without
82 * reciprocal condition numbers
83 * (with sorting of eigenvalues)
84 *
85 * (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
86 * 1/ulp otherwise
87 * If workspace sufficient, also compare VS with and without
88 * reciprocal condition numbers
89 * (with sorting of eigenvalues)
90 *
91 * (13) if sorting worked and SDIM is the number of
92 * eigenvalues which were SELECTed
93 * If workspace sufficient, also compare SDIM with and
94 * without reciprocal condition numbers
95 *
96 * (14) if RCONDE the same no matter if VS and/or RCONDV computed
97 *
98 * (15) if RCONDV the same no matter if VS and/or RCONDE computed
99 *
100 * The "sizes" are specified by an array NN(1:NSIZES); the value of
101 * each element NN(j) specifies one size.
102 * The "types" are specified by a logical array DOTYPE( 1:NTYPES );
103 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
104 * Currently, the list of possible types is:
105 *
106 * (1) The zero matrix.
107 * (2) The identity matrix.
108 * (3) A (transposed) Jordan block, with 1's on the diagonal.
109 *
110 * (4) A diagonal matrix with evenly spaced entries
111 * 1, ..., ULP and random signs.
112 * (ULP = (first number larger than 1) - 1 )
113 * (5) A diagonal matrix with geometrically spaced entries
114 * 1, ..., ULP and random signs.
115 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
116 * and random signs.
117 *
118 * (7) Same as (4), but multiplied by a constant near
119 * the overflow threshold
120 * (8) Same as (4), but multiplied by a constant near
121 * the underflow threshold
122 *
123 * (9) A matrix of the form U' T U, where U is orthogonal and
124 * T has evenly spaced entries 1, ..., ULP with random signs
125 * on the diagonal and random O(1) entries in the upper
126 * triangle.
127 *
128 * (10) A matrix of the form U' T U, where U is orthogonal and
129 * T has geometrically spaced entries 1, ..., ULP with random
130 * signs on the diagonal and random O(1) entries in the upper
131 * triangle.
132 *
133 * (11) A matrix of the form U' T U, where U is orthogonal and
134 * T has "clustered" entries 1, ULP,..., ULP with random
135 * signs on the diagonal and random O(1) entries in the upper
136 * triangle.
137 *
138 * (12) A matrix of the form U' T U, where U is orthogonal and
139 * T has real or complex conjugate paired eigenvalues randomly
140 * chosen from ( ULP, 1 ) and random O(1) entries in the upper
141 * triangle.
142 *
143 * (13) A matrix of the form X' T X, where X has condition
144 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
145 * with random signs on the diagonal and random O(1) entries
146 * in the upper triangle.
147 *
148 * (14) A matrix of the form X' T X, where X has condition
149 * SQRT( ULP ) and T has geometrically spaced entries
150 * 1, ..., ULP with random signs on the diagonal and random
151 * O(1) entries in the upper triangle.
152 *
153 * (15) A matrix of the form X' T X, where X has condition
154 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
155 * with random signs on the diagonal and random O(1) entries
156 * in the upper triangle.
157 *
158 * (16) A matrix of the form X' T X, where X has condition
159 * SQRT( ULP ) and T has real or complex conjugate paired
160 * eigenvalues randomly chosen from ( ULP, 1 ) and random
161 * O(1) entries in the upper triangle.
162 *
163 * (17) Same as (16), but multiplied by a constant
164 * near the overflow threshold
165 * (18) Same as (16), but multiplied by a constant
166 * near the underflow threshold
167 *
168 * (19) Nonsymmetric matrix with random entries chosen from (-1,1).
169 * If N is at least 4, all entries in first two rows and last
170 * row, and first column and last two columns are zero.
171 * (20) Same as (19), but multiplied by a constant
172 * near the overflow threshold
173 * (21) Same as (19), but multiplied by a constant
174 * near the underflow threshold
175 *
176 * In addition, an input file will be read from logical unit number
177 * NIUNIT. The file contains matrices along with precomputed
178 * eigenvalues and reciprocal condition numbers for the eigenvalue
179 * average and right invariant subspace. For these matrices, in
180 * addition to tests (1) to (15) we will compute the following two
181 * tests:
182 *
183 * (16) |RCONDE - RCDEIN| / cond(RCONDE)
184 *
185 * RCONDE is the reciprocal average eigenvalue condition number
186 * computed by SGEESX and RCDEIN (the precomputed true value)
187 * is supplied as input. cond(RCONDE) is the condition number
188 * of RCONDE, and takes errors in computing RCONDE into account,
189 * so that the resulting quantity should be O(ULP). cond(RCONDE)
190 * is essentially given by norm(A)/RCONDV.
191 *
192 * (17) |RCONDV - RCDVIN| / cond(RCONDV)
193 *
194 * RCONDV is the reciprocal right invariant subspace condition
195 * number computed by SGEESX and RCDVIN (the precomputed true
196 * value) is supplied as input. cond(RCONDV) is the condition
197 * number of RCONDV, and takes errors in computing RCONDV into
198 * account, so that the resulting quantity should be O(ULP).
199 * cond(RCONDV) is essentially given by norm(A)/RCONDE.
200 *
201 * Arguments
202 * =========
203 *
204 * NSIZES (input) INTEGER
205 * The number of sizes of matrices to use. NSIZES must be at
206 * least zero. If it is zero, no randomly generated matrices
207 * are tested, but any test matrices read from NIUNIT will be
208 * tested.
209 *
210 * NN (input) INTEGER array, dimension (NSIZES)
211 * An array containing the sizes to be used for the matrices.
212 * Zero values will be skipped. The values must be at least
213 * zero.
214 *
215 * NTYPES (input) INTEGER
216 * The number of elements in DOTYPE. NTYPES must be at least
217 * zero. If it is zero, no randomly generated test matrices
218 * are tested, but and test matrices read from NIUNIT will be
219 * tested. If it is MAXTYP+1 and NSIZES is 1, then an
220 * additional type, MAXTYP+1 is defined, which is to use
221 * whatever matrix is in A. This is only useful if
222 * DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
223 *
224 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
225 * If DOTYPE(j) is .TRUE., then for each size in NN a
226 * matrix of that size and of type j will be generated.
227 * If NTYPES is smaller than the maximum number of types
228 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
229 * MAXTYP will not be generated. If NTYPES is larger
230 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
231 * will be ignored.
232 *
233 * ISEED (input/output) INTEGER array, dimension (4)
234 * On entry ISEED specifies the seed of the random number
235 * generator. The array elements should be between 0 and 4095;
236 * if not they will be reduced mod 4096. Also, ISEED(4) must
237 * be odd. The random number generator uses a linear
238 * congruential sequence limited to small integers, and so
239 * should produce machine independent random numbers. The
240 * values of ISEED are changed on exit, and can be used in the
241 * next call to SDRVSX to continue the same random number
242 * sequence.
243 *
244 * THRESH (input) REAL
245 * A test will count as "failed" if the "error", computed as
246 * described above, exceeds THRESH. Note that the error
247 * is scaled to be O(1), so THRESH should be a reasonably
248 * small multiple of 1, e.g., 10 or 100. In particular,
249 * it should not depend on the precision (single vs. double)
250 * or the size of the matrix. It must be at least zero.
251 *
252 * NIUNIT (input) INTEGER
253 * The FORTRAN unit number for reading in the data file of
254 * problems to solve.
255 *
256 * NOUNIT (input) INTEGER
257 * The FORTRAN unit number for printing out error messages
258 * (e.g., if a routine returns INFO not equal to 0.)
259 *
260 * A (workspace) REAL array, dimension (LDA, max(NN))
261 * Used to hold the matrix whose eigenvalues are to be
262 * computed. On exit, A contains the last matrix actually used.
263 *
264 * LDA (input) INTEGER
265 * The leading dimension of A, and H. LDA must be at
266 * least 1 and at least max( NN ).
267 *
268 * H (workspace) REAL array, dimension (LDA, max(NN))
269 * Another copy of the test matrix A, modified by SGEESX.
270 *
271 * HT (workspace) REAL array, dimension (LDA, max(NN))
272 * Yet another copy of the test matrix A, modified by SGEESX.
273 *
274 * WR (workspace) REAL array, dimension (max(NN))
275 * WI (workspace) REAL array, dimension (max(NN))
276 * The real and imaginary parts of the eigenvalues of A.
277 * On exit, WR + WI*i are the eigenvalues of the matrix in A.
278 *
279 * WRT (workspace) REAL array, dimension (max(NN))
280 * WIT (workspace) REAL array, dimension (max(NN))
281 * Like WR, WI, these arrays contain the eigenvalues of A,
282 * but those computed when SGEESX only computes a partial
283 * eigendecomposition, i.e. not Schur vectors
284 *
285 * WRTMP (workspace) REAL array, dimension (max(NN))
286 * WITMP (workspace) REAL array, dimension (max(NN))
287 * More temporary storage for eigenvalues.
288 *
289 * VS (workspace) REAL array, dimension (LDVS, max(NN))
290 * VS holds the computed Schur vectors.
291 *
292 * LDVS (input) INTEGER
293 * Leading dimension of VS. Must be at least max(1,max(NN)).
294 *
295 * VS1 (workspace) REAL array, dimension (LDVS, max(NN))
296 * VS1 holds another copy of the computed Schur vectors.
297 *
298 * RESULT (output) REAL array, dimension (17)
299 * The values computed by the 17 tests described above.
300 * The values are currently limited to 1/ulp, to avoid overflow.
301 *
302 * WORK (workspace) REAL array, dimension (LWORK)
303 *
304 * LWORK (input) INTEGER
305 * The number of entries in WORK. This must be at least
306 * max(3*NN(j),2*NN(j)**2) for all j.
307 *
308 * IWORK (workspace) INTEGER array, dimension (max(NN)*max(NN))
309 *
310 * INFO (output) INTEGER
311 * If 0, successful exit.
312 * <0, input parameter -INFO is incorrect
313 * >0, SLATMR, SLATMS, SLATME or SGET24 returned an error
314 * code and INFO is its absolute value
315 *
316 *-----------------------------------------------------------------------
317 *
318 * Some Local Variables and Parameters:
319 * ---- ----- --------- --- ----------
320 * ZERO, ONE Real 0 and 1.
321 * MAXTYP The number of types defined.
322 * NMAX Largest value in NN.
323 * NERRS The number of tests which have exceeded THRESH
324 * COND, CONDS,
325 * IMODE Values to be passed to the matrix generators.
326 * ANORM Norm of A; passed to matrix generators.
327 *
328 * OVFL, UNFL Overflow and underflow thresholds.
329 * ULP, ULPINV Finest relative precision and its inverse.
330 * RTULP, RTULPI Square roots of the previous 4 values.
331 * The following four arrays decode JTYPE:
332 * KTYPE(j) The general type (1-10) for type "j".
333 * KMODE(j) The MODE value to be passed to the matrix
334 * generator for type "j".
335 * KMAGN(j) The order of magnitude ( O(1),
336 * O(overflow^(1/2) ), O(underflow^(1/2) )
337 * KCONDS(j) Selectw whether CONDS is to be 1 or
338 * 1/sqrt(ulp). (0 means irrelevant.)
339 *
340 * =====================================================================
341 *
342 * .. Parameters ..
343 REAL ZERO, ONE
344 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
345 INTEGER MAXTYP
346 PARAMETER ( MAXTYP = 21 )
347 * ..
348 * .. Local Scalars ..
349 LOGICAL BADNN
350 CHARACTER*3 PATH
351 INTEGER I, IINFO, IMODE, ITYPE, IWK, J, JCOL, JSIZE,
352 $ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX,
353 $ NNWORK, NSLCT, NTEST, NTESTF, NTESTT
354 REAL ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
355 $ RTULP, RTULPI, ULP, ULPINV, UNFL
356 * ..
357 * .. Local Arrays ..
358 CHARACTER ADUMMA( 1 )
359 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
360 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ),
361 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
362 * ..
363 * .. Arrays in Common ..
364 LOGICAL SELVAL( 20 )
365 REAL SELWI( 20 ), SELWR( 20 )
366 * ..
367 * .. Scalars in Common ..
368 INTEGER SELDIM, SELOPT
369 * ..
370 * .. Common blocks ..
371 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
372 * ..
373 * .. External Functions ..
374 REAL SLAMCH
375 EXTERNAL SLAMCH
376 * ..
377 * .. External Subroutines ..
378 EXTERNAL SGET24, SLABAD, SLASUM, SLATME, SLATMR, SLATMS,
379 $ SLASET, XERBLA
380 * ..
381 * .. Intrinsic Functions ..
382 INTRINSIC ABS, MAX, MIN, SQRT
383 * ..
384 * .. Data statements ..
385 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
386 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
387 $ 3, 1, 2, 3 /
388 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
389 $ 1, 5, 5, 5, 4, 3, 1 /
390 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
391 * ..
392 * .. Executable Statements ..
393 *
394 PATH( 1: 1 ) = 'Single precision'
395 PATH( 2: 3 ) = 'SX'
396 *
397 * Check for errors
398 *
399 NTESTT = 0
400 NTESTF = 0
401 INFO = 0
402 *
403 * Important constants
404 *
405 BADNN = .FALSE.
406 *
407 * 12 is the largest dimension in the input file of precomputed
408 * problems
409 *
410 NMAX = 12
411 DO 10 J = 1, NSIZES
412 NMAX = MAX( NMAX, NN( J ) )
413 IF( NN( J ).LT.0 )
414 $ BADNN = .TRUE.
415 10 CONTINUE
416 *
417 * Check for errors
418 *
419 IF( NSIZES.LT.0 ) THEN
420 INFO = -1
421 ELSE IF( BADNN ) THEN
422 INFO = -2
423 ELSE IF( NTYPES.LT.0 ) THEN
424 INFO = -3
425 ELSE IF( THRESH.LT.ZERO ) THEN
426 INFO = -6
427 ELSE IF( NIUNIT.LE.0 ) THEN
428 INFO = -7
429 ELSE IF( NOUNIT.LE.0 ) THEN
430 INFO = -8
431 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
432 INFO = -10
433 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN
434 INFO = -20
435 ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN
436 INFO = -24
437 END IF
438 *
439 IF( INFO.NE.0 ) THEN
440 CALL XERBLA( 'SDRVSX', -INFO )
441 RETURN
442 END IF
443 *
444 * If nothing to do check on NIUNIT
445 *
446 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
447 $ GO TO 150
448 *
449 * More Important constants
450 *
451 UNFL = SLAMCH( 'Safe minimum' )
452 OVFL = ONE / UNFL
453 CALL SLABAD( UNFL, OVFL )
454 ULP = SLAMCH( 'Precision' )
455 ULPINV = ONE / ULP
456 RTULP = SQRT( ULP )
457 RTULPI = ONE / RTULP
458 *
459 * Loop over sizes, types
460 *
461 NERRS = 0
462 *
463 DO 140 JSIZE = 1, NSIZES
464 N = NN( JSIZE )
465 IF( NSIZES.NE.1 ) THEN
466 MTYPES = MIN( MAXTYP, NTYPES )
467 ELSE
468 MTYPES = MIN( MAXTYP+1, NTYPES )
469 END IF
470 *
471 DO 130 JTYPE = 1, MTYPES
472 IF( .NOT.DOTYPE( JTYPE ) )
473 $ GO TO 130
474 *
475 * Save ISEED in case of an error.
476 *
477 DO 20 J = 1, 4
478 IOLDSD( J ) = ISEED( J )
479 20 CONTINUE
480 *
481 * Compute "A"
482 *
483 * Control parameters:
484 *
485 * KMAGN KCONDS KMODE KTYPE
486 * =1 O(1) 1 clustered 1 zero
487 * =2 large large clustered 2 identity
488 * =3 small exponential Jordan
489 * =4 arithmetic diagonal, (w/ eigenvalues)
490 * =5 random log symmetric, w/ eigenvalues
491 * =6 random general, w/ eigenvalues
492 * =7 random diagonal
493 * =8 random symmetric
494 * =9 random general
495 * =10 random triangular
496 *
497 IF( MTYPES.GT.MAXTYP )
498 $ GO TO 90
499 *
500 ITYPE = KTYPE( JTYPE )
501 IMODE = KMODE( JTYPE )
502 *
503 * Compute norm
504 *
505 GO TO ( 30, 40, 50 )KMAGN( JTYPE )
506 *
507 30 CONTINUE
508 ANORM = ONE
509 GO TO 60
510 *
511 40 CONTINUE
512 ANORM = OVFL*ULP
513 GO TO 60
514 *
515 50 CONTINUE
516 ANORM = UNFL*ULPINV
517 GO TO 60
518 *
519 60 CONTINUE
520 *
521 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
522 IINFO = 0
523 COND = ULPINV
524 *
525 * Special Matrices -- Identity & Jordan block
526 *
527 * Zero
528 *
529 IF( ITYPE.EQ.1 ) THEN
530 IINFO = 0
531 *
532 ELSE IF( ITYPE.EQ.2 ) THEN
533 *
534 * Identity
535 *
536 DO 70 JCOL = 1, N
537 A( JCOL, JCOL ) = ANORM
538 70 CONTINUE
539 *
540 ELSE IF( ITYPE.EQ.3 ) THEN
541 *
542 * Jordan Block
543 *
544 DO 80 JCOL = 1, N
545 A( JCOL, JCOL ) = ANORM
546 IF( JCOL.GT.1 )
547 $ A( JCOL, JCOL-1 ) = ONE
548 80 CONTINUE
549 *
550 ELSE IF( ITYPE.EQ.4 ) THEN
551 *
552 * Diagonal Matrix, [Eigen]values Specified
553 *
554 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
555 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
556 $ IINFO )
557 *
558 ELSE IF( ITYPE.EQ.5 ) THEN
559 *
560 * Symmetric, eigenvalues specified
561 *
562 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
563 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
564 $ IINFO )
565 *
566 ELSE IF( ITYPE.EQ.6 ) THEN
567 *
568 * General, eigenvalues specified
569 *
570 IF( KCONDS( JTYPE ).EQ.1 ) THEN
571 CONDS = ONE
572 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
573 CONDS = RTULPI
574 ELSE
575 CONDS = ZERO
576 END IF
577 *
578 ADUMMA( 1 ) = ' '
579 CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE,
580 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4,
581 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ),
582 $ IINFO )
583 *
584 ELSE IF( ITYPE.EQ.7 ) THEN
585 *
586 * Diagonal, random eigenvalues
587 *
588 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
589 $ 'T', 'N', WORK( N+1 ), 1, ONE,
590 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
591 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
592 *
593 ELSE IF( ITYPE.EQ.8 ) THEN
594 *
595 * Symmetric, random eigenvalues
596 *
597 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
598 $ 'T', 'N', WORK( N+1 ), 1, ONE,
599 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
600 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
601 *
602 ELSE IF( ITYPE.EQ.9 ) THEN
603 *
604 * General, random eigenvalues
605 *
606 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
607 $ 'T', 'N', WORK( N+1 ), 1, ONE,
608 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
609 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
610 IF( N.GE.4 ) THEN
611 CALL SLASET( 'Full', 2, N, ZERO, ZERO, A, LDA )
612 CALL SLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ),
613 $ LDA )
614 CALL SLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ),
615 $ LDA )
616 CALL SLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ),
617 $ LDA )
618 END IF
619 *
620 ELSE IF( ITYPE.EQ.10 ) THEN
621 *
622 * Triangular, random eigenvalues
623 *
624 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
625 $ 'T', 'N', WORK( N+1 ), 1, ONE,
626 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
627 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
628 *
629 ELSE
630 *
631 IINFO = 1
632 END IF
633 *
634 IF( IINFO.NE.0 ) THEN
635 WRITE( NOUNIT, FMT = 9991 )'Generator', IINFO, N, JTYPE,
636 $ IOLDSD
637 INFO = ABS( IINFO )
638 RETURN
639 END IF
640 *
641 90 CONTINUE
642 *
643 * Test for minimal and generous workspace
644 *
645 DO 120 IWK = 1, 2
646 IF( IWK.EQ.1 ) THEN
647 NNWORK = 3*N
648 ELSE
649 NNWORK = MAX( 3*N, 2*N*N )
650 END IF
651 NNWORK = MAX( NNWORK, 1 )
652 *
653 CALL SGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N,
654 $ A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP,
655 $ WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT,
656 $ ISLCT, RESULT, WORK, NNWORK, IWORK, BWORK,
657 $ INFO )
658 *
659 * Check for RESULT(j) > THRESH
660 *
661 NTEST = 0
662 NFAIL = 0
663 DO 100 J = 1, 15
664 IF( RESULT( J ).GE.ZERO )
665 $ NTEST = NTEST + 1
666 IF( RESULT( J ).GE.THRESH )
667 $ NFAIL = NFAIL + 1
668 100 CONTINUE
669 *
670 IF( NFAIL.GT.0 )
671 $ NTESTF = NTESTF + 1
672 IF( NTESTF.EQ.1 ) THEN
673 WRITE( NOUNIT, FMT = 9999 )PATH
674 WRITE( NOUNIT, FMT = 9998 )
675 WRITE( NOUNIT, FMT = 9997 )
676 WRITE( NOUNIT, FMT = 9996 )
677 WRITE( NOUNIT, FMT = 9995 )THRESH
678 WRITE( NOUNIT, FMT = 9994 )
679 NTESTF = 2
680 END IF
681 *
682 DO 110 J = 1, 15
683 IF( RESULT( J ).GE.THRESH ) THEN
684 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
685 $ J, RESULT( J )
686 END IF
687 110 CONTINUE
688 *
689 NERRS = NERRS + NFAIL
690 NTESTT = NTESTT + NTEST
691 *
692 120 CONTINUE
693 130 CONTINUE
694 140 CONTINUE
695 *
696 150 CONTINUE
697 *
698 * Read in data from file to check accuracy of condition estimation
699 * Read input data until N=0
700 *
701 JTYPE = 0
702 160 CONTINUE
703 READ( NIUNIT, FMT = *, END = 200 )N, NSLCT
704 IF( N.EQ.0 )
705 $ GO TO 200
706 JTYPE = JTYPE + 1
707 ISEED( 1 ) = JTYPE
708 IF( NSLCT.GT.0 )
709 $ READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT )
710 DO 170 I = 1, N
711 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
712 170 CONTINUE
713 READ( NIUNIT, FMT = * )RCDEIN, RCDVIN
714 *
715 CALL SGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT,
716 $ WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1,
717 $ RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK,
718 $ IWORK, BWORK, INFO )
719 *
720 * Check for RESULT(j) > THRESH
721 *
722 NTEST = 0
723 NFAIL = 0
724 DO 180 J = 1, 17
725 IF( RESULT( J ).GE.ZERO )
726 $ NTEST = NTEST + 1
727 IF( RESULT( J ).GE.THRESH )
728 $ NFAIL = NFAIL + 1
729 180 CONTINUE
730 *
731 IF( NFAIL.GT.0 )
732 $ NTESTF = NTESTF + 1
733 IF( NTESTF.EQ.1 ) THEN
734 WRITE( NOUNIT, FMT = 9999 )PATH
735 WRITE( NOUNIT, FMT = 9998 )
736 WRITE( NOUNIT, FMT = 9997 )
737 WRITE( NOUNIT, FMT = 9996 )
738 WRITE( NOUNIT, FMT = 9995 )THRESH
739 WRITE( NOUNIT, FMT = 9994 )
740 NTESTF = 2
741 END IF
742 DO 190 J = 1, 17
743 IF( RESULT( J ).GE.THRESH ) THEN
744 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J )
745 END IF
746 190 CONTINUE
747 *
748 NERRS = NERRS + NFAIL
749 NTESTT = NTESTT + NTEST
750 GO TO 160
751 200 CONTINUE
752 *
753 * Summary
754 *
755 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
756 *
757 9999 FORMAT( / 1X, A3, ' -- Real Schur Form Decomposition Expert ',
758 $ 'Driver', / ' Matrix types (see SDRVSX for details):' )
759 *
760 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
761 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
762 $ / ' 2=Identity matrix. ', ' 6=Diagona',
763 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
764 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
765 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
766 $ 'mall, evenly spaced.' )
767 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
768 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
769 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
770 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
771 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
772 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
773 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
774 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
775 $ ' complx ' )
776 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
777 $ 'with small random entries.', / ' 20=Matrix with large ran',
778 $ 'dom entries. ', / )
779 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
780 $ / ' ( A denotes A on input and T denotes A on output)',
781 $ / / ' 1 = 0 if T in Schur form (no sort), ',
782 $ ' 1/ulp otherwise', /
783 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
784 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
785 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
786 $ ' 1/ulp otherwise', /
787 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
788 $ ' 1/ulp otherwise', /
789 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
790 $ ', 1/ulp otherwise' )
791 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
792 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
793 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
794 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
795 $ ' 1/ulp otherwise', /
796 $ ' 11 = 0 if T same no matter what else computed (sort),',
797 $ ' 1/ulp otherwise', /
798 $ ' 12 = 0 if WR, WI same no matter what else computed ',
799 $ '(sort), 1/ulp otherwise', /
800 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise',
801 $ / ' 14 = 0 if RCONDE same no matter what else computed,',
802 $ ' 1/ulp otherwise', /
803 $ ' 15 = 0 if RCONDv same no matter what else computed,',
804 $ ' 1/ulp otherwise', /
805 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
806 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
807 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
808 $ ' type ', I2, ', test(', I2, ')=', G10.3 )
809 9992 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=',
810 $ G10.3 )
811 9991 FORMAT( ' SDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
812 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
813 *
814 RETURN
815 *
816 * End of SDRVSX
817 *
818 END