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