1 SUBROUTINE CDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS,
3 $ LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK,
4 $ 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 ), NN( * )
18 REAL RESULT( 17 ), RWORK( * )
19 COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
20 $ VS( LDVS, * ), VS1( LDVS, * ), W( * ),
21 $ WORK( * ), WT( * ), WTMP( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * CDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
28 * expert driver CGEESX.
29 *
30 * CDRVSX 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 CDRVSX 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 W 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 W are eigenvalues of T
74 * 1/ulp otherwise
75 * If workspace sufficient, also compare W 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 complex angles.
112 * (ULP = (first number larger than 1) - 1 )
113 * (5) A diagonal matrix with geometrically spaced entries
114 * 1, ..., ULP and random complex angles.
115 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
116 * and random complex angles.
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 unitary and
124 * T has evenly spaced entries 1, ..., ULP with random
125 * complex angles on the diagonal and random O(1) entries in
126 * the upper triangle.
127 *
128 * (10) A matrix of the form U' T U, where U is unitary and
129 * T has geometrically spaced entries 1, ..., ULP with random
130 * complex angles on the diagonal and random O(1) entries in
131 * the upper 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 * complex angles on the diagonal and random O(1) entries in
136 * the upper triangle.
137 *
138 * (12) A matrix of the form U' T U, where U is unitary and
139 * T has complex eigenvalues randomly chosen from
140 * ULP < |z| < 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 complex angles on the diagonal and random O(1)
146 * entries 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 complex angles on the diagonal
151 * and random 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 complex angles on the diagonal and random O(1)
156 * entries 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 complex eigenvalues randomly chosen
160 * from ULP < |z| < 1 and random O(1) entries in the upper
161 * 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 CGEESX 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 CGEESX 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 CDRVSX 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) COMPLEX 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) COMPLEX array, dimension (LDA, max(NN))
269 * Another copy of the test matrix A, modified by CGEESX.
270 *
271 * HT (workspace) COMPLEX array, dimension (LDA, max(NN))
272 * Yet another copy of the test matrix A, modified by CGEESX.
273 *
274 * W (workspace) COMPLEX array, dimension (max(NN))
275 * The computed eigenvalues of A.
276 *
277 * WT (workspace) COMPLEX array, dimension (max(NN))
278 * Like W, this array contains the eigenvalues of A,
279 * but those computed when CGEESX only computes a partial
280 * eigendecomposition, i.e. not Schur vectors
281 *
282 * WTMP (workspace) COMPLEX array, dimension (max(NN))
283 * More temporary storage for eigenvalues.
284 *
285 * VS (workspace) COMPLEX array, dimension (LDVS, max(NN))
286 * VS holds the computed Schur vectors.
287 *
288 * LDVS (input) INTEGER
289 * Leading dimension of VS. Must be at least max(1,max(NN)).
290 *
291 * VS1 (workspace) COMPLEX array, dimension (LDVS, max(NN))
292 * VS1 holds another copy of the computed Schur vectors.
293 *
294 * RESULT (output) REAL array, dimension (17)
295 * The values computed by the 17 tests described above.
296 * The values are currently limited to 1/ulp, to avoid overflow.
297 *
298 * WORK (workspace) COMPLEX array, dimension (LWORK)
299 *
300 * LWORK (input) INTEGER
301 * The number of entries in WORK. This must be at least
302 * max(1,2*NN(j)**2) for all j.
303 *
304 * RWORK (workspace) REAL array, dimension (max(NN))
305 *
306 * BWORK (workspace) LOGICAL array, dimension (max(NN))
307 *
308 * INFO (output) INTEGER
309 * If 0, successful exit.
310 * <0, input parameter -INFO is incorrect
311 * >0, CLATMR, CLATMS, CLATME or CGET24 returned an error
312 * code and INFO is its absolute value
313 *
314 *-----------------------------------------------------------------------
315 *
316 * Some Local Variables and Parameters:
317 * ---- ----- --------- --- ----------
318 * ZERO, ONE Real 0 and 1.
319 * MAXTYP The number of types defined.
320 * NMAX Largest value in NN.
321 * NERRS The number of tests which have exceeded THRESH
322 * COND, CONDS,
323 * IMODE Values to be passed to the matrix generators.
324 * ANORM Norm of A; passed to matrix generators.
325 *
326 * OVFL, UNFL Overflow and underflow thresholds.
327 * ULP, ULPINV Finest relative precision and its inverse.
328 * RTULP, RTULPI Square roots of the previous 4 values.
329 * The following four arrays decode JTYPE:
330 * KTYPE(j) The general type (1-10) for type "j".
331 * KMODE(j) The MODE value to be passed to the matrix
332 * generator for type "j".
333 * KMAGN(j) The order of magnitude ( O(1),
334 * O(overflow^(1/2) ), O(underflow^(1/2) )
335 * KCONDS(j) Selectw whether CONDS is to be 1 or
336 * 1/sqrt(ulp). (0 means irrelevant.)
337 *
338 * =====================================================================
339 *
340 * .. Parameters ..
341 COMPLEX CZERO
342 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
343 COMPLEX CONE
344 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
345 REAL ZERO, ONE
346 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
347 INTEGER MAXTYP
348 PARAMETER ( MAXTYP = 21 )
349 * ..
350 * .. Local Scalars ..
351 LOGICAL BADNN
352 CHARACTER*3 PATH
353 INTEGER I, IINFO, IMODE, ISRT, ITYPE, IWK, J, JCOL,
354 $ JSIZE, JTYPE, MTYPES, N, NERRS, NFAIL,
355 $ NMAX, NNWORK, NSLCT, NTEST, NTESTF, NTESTT
356 REAL ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
357 $ RTULP, RTULPI, ULP, ULPINV, UNFL
358 * ..
359 * .. Local Arrays ..
360 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
361 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ),
362 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
363 * ..
364 * .. Arrays in Common ..
365 LOGICAL SELVAL( 20 )
366 REAL SELWI( 20 ), SELWR( 20 )
367 * ..
368 * .. Scalars in Common ..
369 INTEGER SELDIM, SELOPT
370 * ..
371 * .. Common blocks ..
372 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
373 * ..
374 * .. External Functions ..
375 REAL SLAMCH
376 EXTERNAL SLAMCH
377 * ..
378 * .. External Subroutines ..
379 EXTERNAL CGET24, CLATME, CLATMR, CLATMS, CLASET, SLABAD,
380 $ SLASUM, XERBLA
381 * ..
382 * .. Intrinsic Functions ..
383 INTRINSIC ABS, MAX, MIN, SQRT
384 * ..
385 * .. Data statements ..
386 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
387 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
388 $ 3, 1, 2, 3 /
389 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
390 $ 1, 5, 5, 5, 4, 3, 1 /
391 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
392 * ..
393 * .. Executable Statements ..
394 *
395 PATH( 1: 1 ) = 'Complex precision'
396 PATH( 2: 3 ) = 'SX'
397 *
398 * Check for errors
399 *
400 NTESTT = 0
401 NTESTF = 0
402 INFO = 0
403 *
404 * Important constants
405 *
406 BADNN = .FALSE.
407 *
408 * 8 is the largest dimension in the input file of precomputed
409 * problems
410 *
411 NMAX = 8
412 DO 10 J = 1, NSIZES
413 NMAX = MAX( NMAX, NN( J ) )
414 IF( NN( J ).LT.0 )
415 $ BADNN = .TRUE.
416 10 CONTINUE
417 *
418 * Check for errors
419 *
420 IF( NSIZES.LT.0 ) THEN
421 INFO = -1
422 ELSE IF( BADNN ) THEN
423 INFO = -2
424 ELSE IF( NTYPES.LT.0 ) THEN
425 INFO = -3
426 ELSE IF( THRESH.LT.ZERO ) THEN
427 INFO = -6
428 ELSE IF( NIUNIT.LE.0 ) THEN
429 INFO = -7
430 ELSE IF( NOUNIT.LE.0 ) THEN
431 INFO = -8
432 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
433 INFO = -10
434 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN
435 INFO = -20
436 ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN
437 INFO = -24
438 END IF
439 *
440 IF( INFO.NE.0 ) THEN
441 CALL XERBLA( 'CDRVSX', -INFO )
442 RETURN
443 END IF
444 *
445 * If nothing to do check on NIUNIT
446 *
447 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
448 $ GO TO 150
449 *
450 * More Important constants
451 *
452 UNFL = SLAMCH( 'Safe minimum' )
453 OVFL = ONE / UNFL
454 CALL SLABAD( UNFL, OVFL )
455 ULP = SLAMCH( 'Precision' )
456 ULPINV = ONE / ULP
457 RTULP = SQRT( ULP )
458 RTULPI = ONE / RTULP
459 *
460 * Loop over sizes, types
461 *
462 NERRS = 0
463 *
464 DO 140 JSIZE = 1, NSIZES
465 N = NN( JSIZE )
466 IF( NSIZES.NE.1 ) THEN
467 MTYPES = MIN( MAXTYP, NTYPES )
468 ELSE
469 MTYPES = MIN( MAXTYP+1, NTYPES )
470 END IF
471 *
472 DO 130 JTYPE = 1, MTYPES
473 IF( .NOT.DOTYPE( JTYPE ) )
474 $ GO TO 130
475 *
476 * Save ISEED in case of an error.
477 *
478 DO 20 J = 1, 4
479 IOLDSD( J ) = ISEED( J )
480 20 CONTINUE
481 *
482 * Compute "A"
483 *
484 * Control parameters:
485 *
486 * KMAGN KCONDS KMODE KTYPE
487 * =1 O(1) 1 clustered 1 zero
488 * =2 large large clustered 2 identity
489 * =3 small exponential Jordan
490 * =4 arithmetic diagonal, (w/ eigenvalues)
491 * =5 random log symmetric, w/ eigenvalues
492 * =6 random general, w/ eigenvalues
493 * =7 random diagonal
494 * =8 random symmetric
495 * =9 random general
496 * =10 random triangular
497 *
498 IF( MTYPES.GT.MAXTYP )
499 $ GO TO 90
500 *
501 ITYPE = KTYPE( JTYPE )
502 IMODE = KMODE( JTYPE )
503 *
504 * Compute norm
505 *
506 GO TO ( 30, 40, 50 )KMAGN( JTYPE )
507 *
508 30 CONTINUE
509 ANORM = ONE
510 GO TO 60
511 *
512 40 CONTINUE
513 ANORM = OVFL*ULP
514 GO TO 60
515 *
516 50 CONTINUE
517 ANORM = UNFL*ULPINV
518 GO TO 60
519 *
520 60 CONTINUE
521 *
522 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
523 IINFO = 0
524 COND = ULPINV
525 *
526 * Special Matrices -- Identity & Jordan block
527 *
528 IF( ITYPE.EQ.1 ) THEN
529 *
530 * Zero
531 *
532 IINFO = 0
533 *
534 ELSE IF( ITYPE.EQ.2 ) THEN
535 *
536 * Identity
537 *
538 DO 70 JCOL = 1, N
539 A( JCOL, JCOL ) = ANORM
540 70 CONTINUE
541 *
542 ELSE IF( ITYPE.EQ.3 ) THEN
543 *
544 * Jordan Block
545 *
546 DO 80 JCOL = 1, N
547 A( JCOL, JCOL ) = ANORM
548 IF( JCOL.GT.1 )
549 $ A( JCOL, JCOL-1 ) = CONE
550 80 CONTINUE
551 *
552 ELSE IF( ITYPE.EQ.4 ) THEN
553 *
554 * Diagonal Matrix, [Eigen]values Specified
555 *
556 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
557 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
558 $ IINFO )
559 *
560 ELSE IF( ITYPE.EQ.5 ) THEN
561 *
562 * Symmetric, eigenvalues specified
563 *
564 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
565 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
566 $ IINFO )
567 *
568 ELSE IF( ITYPE.EQ.6 ) THEN
569 *
570 * General, eigenvalues specified
571 *
572 IF( KCONDS( JTYPE ).EQ.1 ) THEN
573 CONDS = ONE
574 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
575 CONDS = RTULPI
576 ELSE
577 CONDS = ZERO
578 END IF
579 *
580 CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
581 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
582 $ A, LDA, WORK( 2*N+1 ), IINFO )
583 *
584 ELSE IF( ITYPE.EQ.7 ) THEN
585 *
586 * Diagonal, random eigenvalues
587 *
588 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
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, IDUMMA, IINFO )
592 *
593 ELSE IF( ITYPE.EQ.8 ) THEN
594 *
595 * Symmetric, random eigenvalues
596 *
597 CALL CLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
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, IDUMMA, IINFO )
601 *
602 ELSE IF( ITYPE.EQ.9 ) THEN
603 *
604 * General, random eigenvalues
605 *
606 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
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, IDUMMA, IINFO )
610 IF( N.GE.4 ) THEN
611 CALL CLASET( 'Full', 2, N, CZERO, CZERO, A, LDA )
612 CALL CLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ),
613 $ LDA )
614 CALL CLASET( 'Full', N-3, 2, CZERO, CZERO,
615 $ A( 3, N-1 ), LDA )
616 CALL CLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ),
617 $ LDA )
618 END IF
619 *
620 ELSE IF( ITYPE.EQ.10 ) THEN
621 *
622 * Triangular, random eigenvalues
623 *
624 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
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, IDUMMA, 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 = 2*N
648 ELSE
649 NNWORK = MAX( 2*N, N*( N+1 ) / 2 )
650 END IF
651 NNWORK = MAX( NNWORK, 1 )
652 *
653 CALL CGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N,
654 $ A, LDA, H, HT, W, WT, WTMP, VS, LDVS, VS1,
655 $ RCDEIN, RCDVIN, NSLCT, ISLCT, 0, RESULT,
656 $ WORK, NNWORK, RWORK, BWORK, INFO )
657 *
658 * Check for RESULT(j) > THRESH
659 *
660 NTEST = 0
661 NFAIL = 0
662 DO 100 J = 1, 15
663 IF( RESULT( J ).GE.ZERO )
664 $ NTEST = NTEST + 1
665 IF( RESULT( J ).GE.THRESH )
666 $ NFAIL = NFAIL + 1
667 100 CONTINUE
668 *
669 IF( NFAIL.GT.0 )
670 $ NTESTF = NTESTF + 1
671 IF( NTESTF.EQ.1 ) THEN
672 WRITE( NOUNIT, FMT = 9999 )PATH
673 WRITE( NOUNIT, FMT = 9998 )
674 WRITE( NOUNIT, FMT = 9997 )
675 WRITE( NOUNIT, FMT = 9996 )
676 WRITE( NOUNIT, FMT = 9995 )THRESH
677 WRITE( NOUNIT, FMT = 9994 )
678 NTESTF = 2
679 END IF
680 *
681 DO 110 J = 1, 15
682 IF( RESULT( J ).GE.THRESH ) THEN
683 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
684 $ J, RESULT( J )
685 END IF
686 110 CONTINUE
687 *
688 NERRS = NERRS + NFAIL
689 NTESTT = NTESTT + NTEST
690 *
691 120 CONTINUE
692 130 CONTINUE
693 140 CONTINUE
694 *
695 150 CONTINUE
696 *
697 * Read in data from file to check accuracy of condition estimation
698 * Read input data until N=0
699 *
700 JTYPE = 0
701 160 CONTINUE
702 READ( NIUNIT, FMT = *, END = 200 )N, NSLCT, ISRT
703 IF( N.EQ.0 )
704 $ GO TO 200
705 JTYPE = JTYPE + 1
706 ISEED( 1 ) = JTYPE
707 READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT )
708 DO 170 I = 1, N
709 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
710 170 CONTINUE
711 READ( NIUNIT, FMT = * )RCDEIN, RCDVIN
712 *
713 CALL CGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT,
714 $ W, WT, WTMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT,
715 $ ISLCT, ISRT, RESULT, WORK, LWORK, RWORK, BWORK,
716 $ INFO )
717 *
718 * Check for RESULT(j) > THRESH
719 *
720 NTEST = 0
721 NFAIL = 0
722 DO 180 J = 1, 17
723 IF( RESULT( J ).GE.ZERO )
724 $ NTEST = NTEST + 1
725 IF( RESULT( J ).GE.THRESH )
726 $ NFAIL = NFAIL + 1
727 180 CONTINUE
728 *
729 IF( NFAIL.GT.0 )
730 $ NTESTF = NTESTF + 1
731 IF( NTESTF.EQ.1 ) THEN
732 WRITE( NOUNIT, FMT = 9999 )PATH
733 WRITE( NOUNIT, FMT = 9998 )
734 WRITE( NOUNIT, FMT = 9997 )
735 WRITE( NOUNIT, FMT = 9996 )
736 WRITE( NOUNIT, FMT = 9995 )THRESH
737 WRITE( NOUNIT, FMT = 9994 )
738 NTESTF = 2
739 END IF
740 DO 190 J = 1, 17
741 IF( RESULT( J ).GE.THRESH ) THEN
742 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J )
743 END IF
744 190 CONTINUE
745 *
746 NERRS = NERRS + NFAIL
747 NTESTT = NTESTT + NTEST
748 GO TO 160
749 200 CONTINUE
750 *
751 * Summary
752 *
753 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
754 *
755 9999 FORMAT( / 1X, A3, ' -- Complex Schur Form Decomposition Expert ',
756 $ 'Driver', / ' Matrix types (see CDRVSX for details): ' )
757 *
758 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
759 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
760 $ / ' 2=Identity matrix. ', ' 6=Diagona',
761 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
762 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
763 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
764 $ 'mall, evenly spaced.' )
765 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
766 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
767 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
768 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
769 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
770 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
771 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
772 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
773 $ ' complx ' )
774 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
775 $ 'with small random entries.', / ' 20=Matrix with large ran',
776 $ 'dom entries. ', / )
777 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
778 $ / ' ( A denotes A on input and T denotes A on output)',
779 $ / / ' 1 = 0 if T in Schur form (no sort), ',
780 $ ' 1/ulp otherwise', /
781 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
782 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
783 $ / ' 4 = 0 if W are eigenvalues of T (no sort),',
784 $ ' 1/ulp otherwise', /
785 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
786 $ ' 1/ulp otherwise', /
787 $ ' 6 = 0 if W same no matter if VS computed (no sort)',
788 $ ', 1/ulp otherwise' )
789 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
790 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
791 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
792 $ / ' 10 = 0 if W are eigenvalues of T (sort),',
793 $ ' 1/ulp otherwise', /
794 $ ' 11 = 0 if T same no matter what else computed (sort),',
795 $ ' 1/ulp otherwise', /
796 $ ' 12 = 0 if W same no matter what else computed ',
797 $ '(sort), 1/ulp otherwise', /
798 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise',
799 $ / ' 14 = 0 if RCONDE same no matter what else computed,',
800 $ ' 1/ulp otherwise', /
801 $ ' 15 = 0 if RCONDv same no matter what else computed,',
802 $ ' 1/ulp otherwise', /
803 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
804 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
805 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
806 $ ' type ', I2, ', test(', I2, ')=', G10.3 )
807 9992 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=',
808 $ G10.3 )
809 9991 FORMAT( ' CDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
810 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
811 *
812 RETURN
813 *
814 * End of CDRVSX
815 *
816 END
2 $ NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS,
3 $ LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK,
4 $ 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 ), NN( * )
18 REAL RESULT( 17 ), RWORK( * )
19 COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
20 $ VS( LDVS, * ), VS1( LDVS, * ), W( * ),
21 $ WORK( * ), WT( * ), WTMP( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * CDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
28 * expert driver CGEESX.
29 *
30 * CDRVSX 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 CDRVSX 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 W 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 W are eigenvalues of T
74 * 1/ulp otherwise
75 * If workspace sufficient, also compare W 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 complex angles.
112 * (ULP = (first number larger than 1) - 1 )
113 * (5) A diagonal matrix with geometrically spaced entries
114 * 1, ..., ULP and random complex angles.
115 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
116 * and random complex angles.
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 unitary and
124 * T has evenly spaced entries 1, ..., ULP with random
125 * complex angles on the diagonal and random O(1) entries in
126 * the upper triangle.
127 *
128 * (10) A matrix of the form U' T U, where U is unitary and
129 * T has geometrically spaced entries 1, ..., ULP with random
130 * complex angles on the diagonal and random O(1) entries in
131 * the upper 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 * complex angles on the diagonal and random O(1) entries in
136 * the upper triangle.
137 *
138 * (12) A matrix of the form U' T U, where U is unitary and
139 * T has complex eigenvalues randomly chosen from
140 * ULP < |z| < 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 complex angles on the diagonal and random O(1)
146 * entries 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 complex angles on the diagonal
151 * and random 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 complex angles on the diagonal and random O(1)
156 * entries 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 complex eigenvalues randomly chosen
160 * from ULP < |z| < 1 and random O(1) entries in the upper
161 * 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 CGEESX 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 CGEESX 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 CDRVSX 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) COMPLEX 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) COMPLEX array, dimension (LDA, max(NN))
269 * Another copy of the test matrix A, modified by CGEESX.
270 *
271 * HT (workspace) COMPLEX array, dimension (LDA, max(NN))
272 * Yet another copy of the test matrix A, modified by CGEESX.
273 *
274 * W (workspace) COMPLEX array, dimension (max(NN))
275 * The computed eigenvalues of A.
276 *
277 * WT (workspace) COMPLEX array, dimension (max(NN))
278 * Like W, this array contains the eigenvalues of A,
279 * but those computed when CGEESX only computes a partial
280 * eigendecomposition, i.e. not Schur vectors
281 *
282 * WTMP (workspace) COMPLEX array, dimension (max(NN))
283 * More temporary storage for eigenvalues.
284 *
285 * VS (workspace) COMPLEX array, dimension (LDVS, max(NN))
286 * VS holds the computed Schur vectors.
287 *
288 * LDVS (input) INTEGER
289 * Leading dimension of VS. Must be at least max(1,max(NN)).
290 *
291 * VS1 (workspace) COMPLEX array, dimension (LDVS, max(NN))
292 * VS1 holds another copy of the computed Schur vectors.
293 *
294 * RESULT (output) REAL array, dimension (17)
295 * The values computed by the 17 tests described above.
296 * The values are currently limited to 1/ulp, to avoid overflow.
297 *
298 * WORK (workspace) COMPLEX array, dimension (LWORK)
299 *
300 * LWORK (input) INTEGER
301 * The number of entries in WORK. This must be at least
302 * max(1,2*NN(j)**2) for all j.
303 *
304 * RWORK (workspace) REAL array, dimension (max(NN))
305 *
306 * BWORK (workspace) LOGICAL array, dimension (max(NN))
307 *
308 * INFO (output) INTEGER
309 * If 0, successful exit.
310 * <0, input parameter -INFO is incorrect
311 * >0, CLATMR, CLATMS, CLATME or CGET24 returned an error
312 * code and INFO is its absolute value
313 *
314 *-----------------------------------------------------------------------
315 *
316 * Some Local Variables and Parameters:
317 * ---- ----- --------- --- ----------
318 * ZERO, ONE Real 0 and 1.
319 * MAXTYP The number of types defined.
320 * NMAX Largest value in NN.
321 * NERRS The number of tests which have exceeded THRESH
322 * COND, CONDS,
323 * IMODE Values to be passed to the matrix generators.
324 * ANORM Norm of A; passed to matrix generators.
325 *
326 * OVFL, UNFL Overflow and underflow thresholds.
327 * ULP, ULPINV Finest relative precision and its inverse.
328 * RTULP, RTULPI Square roots of the previous 4 values.
329 * The following four arrays decode JTYPE:
330 * KTYPE(j) The general type (1-10) for type "j".
331 * KMODE(j) The MODE value to be passed to the matrix
332 * generator for type "j".
333 * KMAGN(j) The order of magnitude ( O(1),
334 * O(overflow^(1/2) ), O(underflow^(1/2) )
335 * KCONDS(j) Selectw whether CONDS is to be 1 or
336 * 1/sqrt(ulp). (0 means irrelevant.)
337 *
338 * =====================================================================
339 *
340 * .. Parameters ..
341 COMPLEX CZERO
342 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
343 COMPLEX CONE
344 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
345 REAL ZERO, ONE
346 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
347 INTEGER MAXTYP
348 PARAMETER ( MAXTYP = 21 )
349 * ..
350 * .. Local Scalars ..
351 LOGICAL BADNN
352 CHARACTER*3 PATH
353 INTEGER I, IINFO, IMODE, ISRT, ITYPE, IWK, J, JCOL,
354 $ JSIZE, JTYPE, MTYPES, N, NERRS, NFAIL,
355 $ NMAX, NNWORK, NSLCT, NTEST, NTESTF, NTESTT
356 REAL ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
357 $ RTULP, RTULPI, ULP, ULPINV, UNFL
358 * ..
359 * .. Local Arrays ..
360 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
361 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ),
362 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
363 * ..
364 * .. Arrays in Common ..
365 LOGICAL SELVAL( 20 )
366 REAL SELWI( 20 ), SELWR( 20 )
367 * ..
368 * .. Scalars in Common ..
369 INTEGER SELDIM, SELOPT
370 * ..
371 * .. Common blocks ..
372 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
373 * ..
374 * .. External Functions ..
375 REAL SLAMCH
376 EXTERNAL SLAMCH
377 * ..
378 * .. External Subroutines ..
379 EXTERNAL CGET24, CLATME, CLATMR, CLATMS, CLASET, SLABAD,
380 $ SLASUM, XERBLA
381 * ..
382 * .. Intrinsic Functions ..
383 INTRINSIC ABS, MAX, MIN, SQRT
384 * ..
385 * .. Data statements ..
386 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
387 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
388 $ 3, 1, 2, 3 /
389 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
390 $ 1, 5, 5, 5, 4, 3, 1 /
391 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
392 * ..
393 * .. Executable Statements ..
394 *
395 PATH( 1: 1 ) = 'Complex precision'
396 PATH( 2: 3 ) = 'SX'
397 *
398 * Check for errors
399 *
400 NTESTT = 0
401 NTESTF = 0
402 INFO = 0
403 *
404 * Important constants
405 *
406 BADNN = .FALSE.
407 *
408 * 8 is the largest dimension in the input file of precomputed
409 * problems
410 *
411 NMAX = 8
412 DO 10 J = 1, NSIZES
413 NMAX = MAX( NMAX, NN( J ) )
414 IF( NN( J ).LT.0 )
415 $ BADNN = .TRUE.
416 10 CONTINUE
417 *
418 * Check for errors
419 *
420 IF( NSIZES.LT.0 ) THEN
421 INFO = -1
422 ELSE IF( BADNN ) THEN
423 INFO = -2
424 ELSE IF( NTYPES.LT.0 ) THEN
425 INFO = -3
426 ELSE IF( THRESH.LT.ZERO ) THEN
427 INFO = -6
428 ELSE IF( NIUNIT.LE.0 ) THEN
429 INFO = -7
430 ELSE IF( NOUNIT.LE.0 ) THEN
431 INFO = -8
432 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
433 INFO = -10
434 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN
435 INFO = -20
436 ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN
437 INFO = -24
438 END IF
439 *
440 IF( INFO.NE.0 ) THEN
441 CALL XERBLA( 'CDRVSX', -INFO )
442 RETURN
443 END IF
444 *
445 * If nothing to do check on NIUNIT
446 *
447 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
448 $ GO TO 150
449 *
450 * More Important constants
451 *
452 UNFL = SLAMCH( 'Safe minimum' )
453 OVFL = ONE / UNFL
454 CALL SLABAD( UNFL, OVFL )
455 ULP = SLAMCH( 'Precision' )
456 ULPINV = ONE / ULP
457 RTULP = SQRT( ULP )
458 RTULPI = ONE / RTULP
459 *
460 * Loop over sizes, types
461 *
462 NERRS = 0
463 *
464 DO 140 JSIZE = 1, NSIZES
465 N = NN( JSIZE )
466 IF( NSIZES.NE.1 ) THEN
467 MTYPES = MIN( MAXTYP, NTYPES )
468 ELSE
469 MTYPES = MIN( MAXTYP+1, NTYPES )
470 END IF
471 *
472 DO 130 JTYPE = 1, MTYPES
473 IF( .NOT.DOTYPE( JTYPE ) )
474 $ GO TO 130
475 *
476 * Save ISEED in case of an error.
477 *
478 DO 20 J = 1, 4
479 IOLDSD( J ) = ISEED( J )
480 20 CONTINUE
481 *
482 * Compute "A"
483 *
484 * Control parameters:
485 *
486 * KMAGN KCONDS KMODE KTYPE
487 * =1 O(1) 1 clustered 1 zero
488 * =2 large large clustered 2 identity
489 * =3 small exponential Jordan
490 * =4 arithmetic diagonal, (w/ eigenvalues)
491 * =5 random log symmetric, w/ eigenvalues
492 * =6 random general, w/ eigenvalues
493 * =7 random diagonal
494 * =8 random symmetric
495 * =9 random general
496 * =10 random triangular
497 *
498 IF( MTYPES.GT.MAXTYP )
499 $ GO TO 90
500 *
501 ITYPE = KTYPE( JTYPE )
502 IMODE = KMODE( JTYPE )
503 *
504 * Compute norm
505 *
506 GO TO ( 30, 40, 50 )KMAGN( JTYPE )
507 *
508 30 CONTINUE
509 ANORM = ONE
510 GO TO 60
511 *
512 40 CONTINUE
513 ANORM = OVFL*ULP
514 GO TO 60
515 *
516 50 CONTINUE
517 ANORM = UNFL*ULPINV
518 GO TO 60
519 *
520 60 CONTINUE
521 *
522 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
523 IINFO = 0
524 COND = ULPINV
525 *
526 * Special Matrices -- Identity & Jordan block
527 *
528 IF( ITYPE.EQ.1 ) THEN
529 *
530 * Zero
531 *
532 IINFO = 0
533 *
534 ELSE IF( ITYPE.EQ.2 ) THEN
535 *
536 * Identity
537 *
538 DO 70 JCOL = 1, N
539 A( JCOL, JCOL ) = ANORM
540 70 CONTINUE
541 *
542 ELSE IF( ITYPE.EQ.3 ) THEN
543 *
544 * Jordan Block
545 *
546 DO 80 JCOL = 1, N
547 A( JCOL, JCOL ) = ANORM
548 IF( JCOL.GT.1 )
549 $ A( JCOL, JCOL-1 ) = CONE
550 80 CONTINUE
551 *
552 ELSE IF( ITYPE.EQ.4 ) THEN
553 *
554 * Diagonal Matrix, [Eigen]values Specified
555 *
556 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
557 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
558 $ IINFO )
559 *
560 ELSE IF( ITYPE.EQ.5 ) THEN
561 *
562 * Symmetric, eigenvalues specified
563 *
564 CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
565 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
566 $ IINFO )
567 *
568 ELSE IF( ITYPE.EQ.6 ) THEN
569 *
570 * General, eigenvalues specified
571 *
572 IF( KCONDS( JTYPE ).EQ.1 ) THEN
573 CONDS = ONE
574 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
575 CONDS = RTULPI
576 ELSE
577 CONDS = ZERO
578 END IF
579 *
580 CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
581 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
582 $ A, LDA, WORK( 2*N+1 ), IINFO )
583 *
584 ELSE IF( ITYPE.EQ.7 ) THEN
585 *
586 * Diagonal, random eigenvalues
587 *
588 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
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, IDUMMA, IINFO )
592 *
593 ELSE IF( ITYPE.EQ.8 ) THEN
594 *
595 * Symmetric, random eigenvalues
596 *
597 CALL CLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
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, IDUMMA, IINFO )
601 *
602 ELSE IF( ITYPE.EQ.9 ) THEN
603 *
604 * General, random eigenvalues
605 *
606 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
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, IDUMMA, IINFO )
610 IF( N.GE.4 ) THEN
611 CALL CLASET( 'Full', 2, N, CZERO, CZERO, A, LDA )
612 CALL CLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ),
613 $ LDA )
614 CALL CLASET( 'Full', N-3, 2, CZERO, CZERO,
615 $ A( 3, N-1 ), LDA )
616 CALL CLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ),
617 $ LDA )
618 END IF
619 *
620 ELSE IF( ITYPE.EQ.10 ) THEN
621 *
622 * Triangular, random eigenvalues
623 *
624 CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
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, IDUMMA, 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 = 2*N
648 ELSE
649 NNWORK = MAX( 2*N, N*( N+1 ) / 2 )
650 END IF
651 NNWORK = MAX( NNWORK, 1 )
652 *
653 CALL CGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N,
654 $ A, LDA, H, HT, W, WT, WTMP, VS, LDVS, VS1,
655 $ RCDEIN, RCDVIN, NSLCT, ISLCT, 0, RESULT,
656 $ WORK, NNWORK, RWORK, BWORK, INFO )
657 *
658 * Check for RESULT(j) > THRESH
659 *
660 NTEST = 0
661 NFAIL = 0
662 DO 100 J = 1, 15
663 IF( RESULT( J ).GE.ZERO )
664 $ NTEST = NTEST + 1
665 IF( RESULT( J ).GE.THRESH )
666 $ NFAIL = NFAIL + 1
667 100 CONTINUE
668 *
669 IF( NFAIL.GT.0 )
670 $ NTESTF = NTESTF + 1
671 IF( NTESTF.EQ.1 ) THEN
672 WRITE( NOUNIT, FMT = 9999 )PATH
673 WRITE( NOUNIT, FMT = 9998 )
674 WRITE( NOUNIT, FMT = 9997 )
675 WRITE( NOUNIT, FMT = 9996 )
676 WRITE( NOUNIT, FMT = 9995 )THRESH
677 WRITE( NOUNIT, FMT = 9994 )
678 NTESTF = 2
679 END IF
680 *
681 DO 110 J = 1, 15
682 IF( RESULT( J ).GE.THRESH ) THEN
683 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
684 $ J, RESULT( J )
685 END IF
686 110 CONTINUE
687 *
688 NERRS = NERRS + NFAIL
689 NTESTT = NTESTT + NTEST
690 *
691 120 CONTINUE
692 130 CONTINUE
693 140 CONTINUE
694 *
695 150 CONTINUE
696 *
697 * Read in data from file to check accuracy of condition estimation
698 * Read input data until N=0
699 *
700 JTYPE = 0
701 160 CONTINUE
702 READ( NIUNIT, FMT = *, END = 200 )N, NSLCT, ISRT
703 IF( N.EQ.0 )
704 $ GO TO 200
705 JTYPE = JTYPE + 1
706 ISEED( 1 ) = JTYPE
707 READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT )
708 DO 170 I = 1, N
709 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
710 170 CONTINUE
711 READ( NIUNIT, FMT = * )RCDEIN, RCDVIN
712 *
713 CALL CGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT,
714 $ W, WT, WTMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT,
715 $ ISLCT, ISRT, RESULT, WORK, LWORK, RWORK, BWORK,
716 $ INFO )
717 *
718 * Check for RESULT(j) > THRESH
719 *
720 NTEST = 0
721 NFAIL = 0
722 DO 180 J = 1, 17
723 IF( RESULT( J ).GE.ZERO )
724 $ NTEST = NTEST + 1
725 IF( RESULT( J ).GE.THRESH )
726 $ NFAIL = NFAIL + 1
727 180 CONTINUE
728 *
729 IF( NFAIL.GT.0 )
730 $ NTESTF = NTESTF + 1
731 IF( NTESTF.EQ.1 ) THEN
732 WRITE( NOUNIT, FMT = 9999 )PATH
733 WRITE( NOUNIT, FMT = 9998 )
734 WRITE( NOUNIT, FMT = 9997 )
735 WRITE( NOUNIT, FMT = 9996 )
736 WRITE( NOUNIT, FMT = 9995 )THRESH
737 WRITE( NOUNIT, FMT = 9994 )
738 NTESTF = 2
739 END IF
740 DO 190 J = 1, 17
741 IF( RESULT( J ).GE.THRESH ) THEN
742 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J )
743 END IF
744 190 CONTINUE
745 *
746 NERRS = NERRS + NFAIL
747 NTESTT = NTESTT + NTEST
748 GO TO 160
749 200 CONTINUE
750 *
751 * Summary
752 *
753 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
754 *
755 9999 FORMAT( / 1X, A3, ' -- Complex Schur Form Decomposition Expert ',
756 $ 'Driver', / ' Matrix types (see CDRVSX for details): ' )
757 *
758 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
759 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
760 $ / ' 2=Identity matrix. ', ' 6=Diagona',
761 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
762 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
763 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
764 $ 'mall, evenly spaced.' )
765 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
766 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
767 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
768 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
769 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
770 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
771 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
772 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
773 $ ' complx ' )
774 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
775 $ 'with small random entries.', / ' 20=Matrix with large ran',
776 $ 'dom entries. ', / )
777 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
778 $ / ' ( A denotes A on input and T denotes A on output)',
779 $ / / ' 1 = 0 if T in Schur form (no sort), ',
780 $ ' 1/ulp otherwise', /
781 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
782 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
783 $ / ' 4 = 0 if W are eigenvalues of T (no sort),',
784 $ ' 1/ulp otherwise', /
785 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
786 $ ' 1/ulp otherwise', /
787 $ ' 6 = 0 if W same no matter if VS computed (no sort)',
788 $ ', 1/ulp otherwise' )
789 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
790 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
791 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
792 $ / ' 10 = 0 if W are eigenvalues of T (sort),',
793 $ ' 1/ulp otherwise', /
794 $ ' 11 = 0 if T same no matter what else computed (sort),',
795 $ ' 1/ulp otherwise', /
796 $ ' 12 = 0 if W same no matter what else computed ',
797 $ '(sort), 1/ulp otherwise', /
798 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise',
799 $ / ' 14 = 0 if RCONDE same no matter what else computed,',
800 $ ' 1/ulp otherwise', /
801 $ ' 15 = 0 if RCONDv same no matter what else computed,',
802 $ ' 1/ulp otherwise', /
803 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
804 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
805 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
806 $ ' type ', I2, ', test(', I2, ')=', G10.3 )
807 9992 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=',
808 $ G10.3 )
809 9991 FORMAT( ' CDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
810 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
811 *
812 RETURN
813 *
814 * End of CDRVSX
815 *
816 END