1 SUBROUTINE ZCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, W1,
3 $ W3, EVECTL, EVECTR, EVECTY, EVECTX, UU, TAU,
4 $ WORK, NWORK, RWORK, IWORK, SELECT, RESULT,
5 $ INFO )
6 *
7 * -- LAPACK test routine (version 3.1.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * February 2007
10 *
11 * .. Scalar Arguments ..
12 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
13 DOUBLE PRECISION THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * ), SELECT( * )
17 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
18 DOUBLE PRECISION RESULT( 14 ), RWORK( * )
19 COMPLEX*16 A( LDA, * ), EVECTL( LDU, * ),
20 $ EVECTR( LDU, * ), EVECTX( LDU, * ),
21 $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
22 $ T2( LDA, * ), TAU( * ), U( LDU, * ),
23 $ UU( LDU, * ), UZ( LDU, * ), W1( * ), W3( * ),
24 $ WORK( * ), Z( LDU, * )
25 * ..
26 *
27 * Purpose
28 * =======
29 *
30 * ZCHKHS checks the nonsymmetric eigenvalue problem routines.
31 *
32 * ZGEHRD factors A as U H U' , where ' means conjugate
33 * transpose, H is hessenberg, and U is unitary.
34 *
35 * ZUNGHR generates the unitary matrix U.
36 *
37 * ZUNMHR multiplies a matrix by the unitary matrix U.
38 *
39 * ZHSEQR factors H as Z T Z' , where Z is unitary and T
40 * is upper triangular. It also computes the eigenvalues,
41 * w(1), ..., w(n); we define a diagonal matrix W whose
42 * (diagonal) entries are the eigenvalues.
43 *
44 * ZTREVC computes the left eigenvector matrix L and the
45 * right eigenvector matrix R for the matrix T. The
46 * columns of L are the complex conjugates of the left
47 * eigenvectors of T. The columns of R are the right
48 * eigenvectors of T. L is lower triangular, and R is
49 * upper triangular.
50 *
51 * ZHSEIN computes the left eigenvector matrix Y and the
52 * right eigenvector matrix X for the matrix H. The
53 * columns of Y are the complex conjugates of the left
54 * eigenvectors of H. The columns of X are the right
55 * eigenvectors of H. Y is lower triangular, and X is
56 * upper triangular.
57 *
58 * When ZCHKHS is called, a number of matrix "sizes" ("n's") and a
59 * number of matrix "types" are specified. For each size ("n")
60 * and each type of matrix, one matrix will be generated and used
61 * to test the nonsymmetric eigenroutines. For each matrix, 14
62 * tests will be performed:
63 *
64 * (1) | A - U H U**H | / ( |A| n ulp )
65 *
66 * (2) | I - UU**H | / ( n ulp )
67 *
68 * (3) | H - Z T Z**H | / ( |H| n ulp )
69 *
70 * (4) | I - ZZ**H | / ( n ulp )
71 *
72 * (5) | A - UZ H (UZ)**H | / ( |A| n ulp )
73 *
74 * (6) | I - UZ (UZ)**H | / ( n ulp )
75 *
76 * (7) | T(Z computed) - T(Z not computed) | / ( |T| ulp )
77 *
78 * (8) | W(Z computed) - W(Z not computed) | / ( |W| ulp )
79 *
80 * (9) | TR - RW | / ( |T| |R| ulp )
81 *
82 * (10) | L**H T - W**H L | / ( |T| |L| ulp )
83 *
84 * (11) | HX - XW | / ( |H| |X| ulp )
85 *
86 * (12) | Y**H H - W**H Y | / ( |H| |Y| ulp )
87 *
88 * (13) | AX - XW | / ( |A| |X| ulp )
89 *
90 * (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
91 *
92 * The "sizes" are specified by an array NN(1:NSIZES); the value of
93 * each element NN(j) specifies one size.
94 * The "types" are specified by a logical array DOTYPE( 1:NTYPES );
95 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
96 * Currently, the list of possible types is:
97 *
98 * (1) The zero matrix.
99 * (2) The identity matrix.
100 * (3) A (transposed) Jordan block, with 1's on the diagonal.
101 *
102 * (4) A diagonal matrix with evenly spaced entries
103 * 1, ..., ULP and random complex angles.
104 * (ULP = (first number larger than 1) - 1 )
105 * (5) A diagonal matrix with geometrically spaced entries
106 * 1, ..., ULP and random complex angles.
107 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
108 * and random complex angles.
109 *
110 * (7) Same as (4), but multiplied by SQRT( overflow threshold )
111 * (8) Same as (4), but multiplied by SQRT( underflow threshold )
112 *
113 * (9) A matrix of the form U' T U, where U is unitary and
114 * T has evenly spaced entries 1, ..., ULP with random complex
115 * angles on the diagonal and random O(1) entries in the upper
116 * triangle.
117 *
118 * (10) A matrix of the form U' T U, where U is unitary and
119 * T has geometrically spaced entries 1, ..., ULP with random
120 * complex angles on the diagonal and random O(1) entries in
121 * the upper triangle.
122 *
123 * (11) A matrix of the form U' T U, where U is unitary and
124 * T has "clustered" entries 1, ULP,..., ULP with random
125 * complex angles on the diagonal and random O(1) entries in
126 * the upper triangle.
127 *
128 * (12) A matrix of the form U' T U, where U is unitary and
129 * T has complex eigenvalues randomly chosen from
130 * ULP < |z| < 1 and random O(1) entries in the upper
131 * triangle.
132 *
133 * (13) A matrix of the form X' T X, where X has condition
134 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
135 * with random complex angles on the diagonal and random O(1)
136 * entries in the upper triangle.
137 *
138 * (14) A matrix of the form X' T X, where X has condition
139 * SQRT( ULP ) and T has geometrically spaced entries
140 * 1, ..., ULP with random complex angles on the diagonal
141 * and random O(1) entries in the upper triangle.
142 *
143 * (15) A matrix of the form X' T X, where X has condition
144 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
145 * with random complex angles on the diagonal and random O(1)
146 * entries in the upper triangle.
147 *
148 * (16) A matrix of the form X' T X, where X has condition
149 * SQRT( ULP ) and T has complex eigenvalues randomly chosen
150 * from ULP < |z| < 1 and random O(1) entries in the upper
151 * triangle.
152 *
153 * (17) Same as (16), but multiplied by SQRT( overflow threshold )
154 * (18) Same as (16), but multiplied by SQRT( underflow threshold )
155 *
156 * (19) Nonsymmetric matrix with random entries chosen from |z| < 1
157 * (20) Same as (19), but multiplied by SQRT( overflow threshold )
158 * (21) Same as (19), but multiplied by SQRT( underflow threshold )
159 *
160 * Arguments
161 * ==========
162 *
163 * NSIZES - INTEGER
164 * The number of sizes of matrices to use. If it is zero,
165 * ZCHKHS does nothing. It must be at least zero.
166 * Not modified.
167 *
168 * NN - INTEGER array, dimension (NSIZES)
169 * An array containing the sizes to be used for the matrices.
170 * Zero values will be skipped. The values must be at least
171 * zero.
172 * Not modified.
173 *
174 * NTYPES - INTEGER
175 * The number of elements in DOTYPE. If it is zero, ZCHKHS
176 * does nothing. It must be at least zero. If it is MAXTYP+1
177 * and NSIZES is 1, then an additional type, MAXTYP+1 is
178 * defined, which is to use whatever matrix is in A. This
179 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
180 * DOTYPE(MAXTYP+1) is .TRUE. .
181 * Not modified.
182 *
183 * DOTYPE - LOGICAL array, dimension (NTYPES)
184 * If DOTYPE(j) is .TRUE., then for each size in NN a
185 * matrix of that size and of type j will be generated.
186 * If NTYPES is smaller than the maximum number of types
187 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
188 * MAXTYP will not be generated. If NTYPES is larger
189 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
190 * will be ignored.
191 * Not modified.
192 *
193 * ISEED - INTEGER array, dimension (4)
194 * On entry ISEED specifies the seed of the random number
195 * generator. The array elements should be between 0 and 4095;
196 * if not they will be reduced mod 4096. Also, ISEED(4) must
197 * be odd. The random number generator uses a linear
198 * congruential sequence limited to small integers, and so
199 * should produce machine independent random numbers. The
200 * values of ISEED are changed on exit, and can be used in the
201 * next call to ZCHKHS to continue the same random number
202 * sequence.
203 * Modified.
204 *
205 * THRESH - DOUBLE PRECISION
206 * A test will count as "failed" if the "error", computed as
207 * described above, exceeds THRESH. Note that the error
208 * is scaled to be O(1), so THRESH should be a reasonably
209 * small multiple of 1, e.g., 10 or 100. In particular,
210 * it should not depend on the precision (single vs. double)
211 * or the size of the matrix. It must be at least zero.
212 * Not modified.
213 *
214 * NOUNIT - INTEGER
215 * The FORTRAN unit number for printing out error messages
216 * (e.g., if a routine returns IINFO not equal to 0.)
217 * Not modified.
218 *
219 * A - COMPLEX*16 array, dimension (LDA,max(NN))
220 * Used to hold the matrix whose eigenvalues are to be
221 * computed. On exit, A contains the last matrix actually
222 * used.
223 * Modified.
224 *
225 * LDA - INTEGER
226 * The leading dimension of A, H, T1 and T2. It must be at
227 * least 1 and at least max( NN ).
228 * Not modified.
229 *
230 * H - COMPLEX*16 array, dimension (LDA,max(NN))
231 * The upper hessenberg matrix computed by ZGEHRD. On exit,
232 * H contains the Hessenberg form of the matrix in A.
233 * Modified.
234 *
235 * T1 - COMPLEX*16 array, dimension (LDA,max(NN))
236 * The Schur (="quasi-triangular") matrix computed by ZHSEQR
237 * if Z is computed. On exit, T1 contains the Schur form of
238 * the matrix in A.
239 * Modified.
240 *
241 * T2 - COMPLEX*16 array, dimension (LDA,max(NN))
242 * The Schur matrix computed by ZHSEQR when Z is not computed.
243 * This should be identical to T1.
244 * Modified.
245 *
246 * LDU - INTEGER
247 * The leading dimension of U, Z, UZ and UU. It must be at
248 * least 1 and at least max( NN ).
249 * Not modified.
250 *
251 * U - COMPLEX*16 array, dimension (LDU,max(NN))
252 * The unitary matrix computed by ZGEHRD.
253 * Modified.
254 *
255 * Z - COMPLEX*16 array, dimension (LDU,max(NN))
256 * The unitary matrix computed by ZHSEQR.
257 * Modified.
258 *
259 * UZ - COMPLEX*16 array, dimension (LDU,max(NN))
260 * The product of U times Z.
261 * Modified.
262 *
263 * W1 - COMPLEX*16 array, dimension (max(NN))
264 * The eigenvalues of A, as computed by a full Schur
265 * decomposition H = Z T Z'. On exit, W1 contains the
266 * eigenvalues of the matrix in A.
267 * Modified.
268 *
269 * W3 - COMPLEX*16 array, dimension (max(NN))
270 * The eigenvalues of A, as computed by a partial Schur
271 * decomposition (Z not computed, T only computed as much
272 * as is necessary for determining eigenvalues). On exit,
273 * W3 contains the eigenvalues of the matrix in A, possibly
274 * perturbed by ZHSEIN.
275 * Modified.
276 *
277 * EVECTL - COMPLEX*16 array, dimension (LDU,max(NN))
278 * The conjugate transpose of the (upper triangular) left
279 * eigenvector matrix for the matrix in T1.
280 * Modified.
281 *
282 * EVEZTR - COMPLEX*16 array, dimension (LDU,max(NN))
283 * The (upper triangular) right eigenvector matrix for the
284 * matrix in T1.
285 * Modified.
286 *
287 * EVECTY - COMPLEX*16 array, dimension (LDU,max(NN))
288 * The conjugate transpose of the left eigenvector matrix
289 * for the matrix in H.
290 * Modified.
291 *
292 * EVECTX - COMPLEX*16 array, dimension (LDU,max(NN))
293 * The right eigenvector matrix for the matrix in H.
294 * Modified.
295 *
296 * UU - COMPLEX*16 array, dimension (LDU,max(NN))
297 * Details of the unitary matrix computed by ZGEHRD.
298 * Modified.
299 *
300 * TAU - COMPLEX*16 array, dimension (max(NN))
301 * Further details of the unitary matrix computed by ZGEHRD.
302 * Modified.
303 *
304 * WORK - COMPLEX*16 array, dimension (NWORK)
305 * Workspace.
306 * Modified.
307 *
308 * NWORK - INTEGER
309 * The number of entries in WORK. NWORK >= 4*NN(j)*NN(j) + 2.
310 *
311 * RWORK - DOUBLE PRECISION array, dimension (max(NN))
312 * Workspace. Could be equivalenced to IWORK, but not SELECT.
313 * Modified.
314 *
315 * IWORK - INTEGER array, dimension (max(NN))
316 * Workspace.
317 * Modified.
318 *
319 * SELECT - LOGICAL array, dimension (max(NN))
320 * Workspace. Could be equivalenced to IWORK, but not RWORK.
321 * Modified.
322 *
323 * RESULT - DOUBLE PRECISION array, dimension (14)
324 * The values computed by the fourteen tests described above.
325 * The values are currently limited to 1/ulp, to avoid
326 * overflow.
327 * Modified.
328 *
329 * INFO - INTEGER
330 * If 0, then everything ran OK.
331 * -1: NSIZES < 0
332 * -2: Some NN(j) < 0
333 * -3: NTYPES < 0
334 * -6: THRESH < 0
335 * -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
336 * -14: LDU < 1 or LDU < NMAX.
337 * -26: NWORK too small.
338 * If ZLATMR, CLATMS, or CLATME returns an error code, the
339 * absolute value of it is returned.
340 * If 1, then ZHSEQR could not find all the shifts.
341 * If 2, then the EISPACK code (for small blocks) failed.
342 * If >2, then 30*N iterations were not enough to find an
343 * eigenvalue or to decompose the problem.
344 * Modified.
345 *
346 *-----------------------------------------------------------------------
347 *
348 * Some Local Variables and Parameters:
349 * ---- ----- --------- --- ----------
350 *
351 * ZERO, ONE Real 0 and 1.
352 * MAXTYP The number of types defined.
353 * MTEST The number of tests defined: care must be taken
354 * that (1) the size of RESULT, (2) the number of
355 * tests actually performed, and (3) MTEST agree.
356 * NTEST The number of tests performed on this matrix
357 * so far. This should be less than MTEST, and
358 * equal to it by the last test. It will be less
359 * if any of the routines being tested indicates
360 * that it could not compute the matrices that
361 * would be tested.
362 * NMAX Largest value in NN.
363 * NMATS The number of matrices generated so far.
364 * NERRS The number of tests which have exceeded THRESH
365 * so far (computed by DLAFTS).
366 * COND, CONDS,
367 * IMODE Values to be passed to the matrix generators.
368 * ANORM Norm of A; passed to matrix generators.
369 *
370 * OVFL, UNFL Overflow and underflow thresholds.
371 * ULP, ULPINV Finest relative precision and its inverse.
372 * RTOVFL, RTUNFL,
373 * RTULP, RTULPI Square roots of the previous 4 values.
374 *
375 * The following four arrays decode JTYPE:
376 * KTYPE(j) The general type (1-10) for type "j".
377 * KMODE(j) The MODE value to be passed to the matrix
378 * generator for type "j".
379 * KMAGN(j) The order of magnitude ( O(1),
380 * O(overflow^(1/2) ), O(underflow^(1/2) )
381 * KCONDS(j) Selects whether CONDS is to be 1 or
382 * 1/sqrt(ulp). (0 means irrelevant.)
383 *
384 * =====================================================================
385 *
386 * .. Parameters ..
387 DOUBLE PRECISION ZERO, ONE
388 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
389 COMPLEX*16 CZERO, CONE
390 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
391 $ CONE = ( 1.0D+0, 0.0D+0 ) )
392 INTEGER MAXTYP
393 PARAMETER ( MAXTYP = 21 )
394 * ..
395 * .. Local Scalars ..
396 LOGICAL BADNN, MATCH
397 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
398 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
399 $ NMATS, NMAX, NTEST, NTESTT
400 DOUBLE PRECISION ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
401 $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL
402 * ..
403 * .. Local Arrays ..
404 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
405 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
406 $ KTYPE( MAXTYP )
407 DOUBLE PRECISION DUMMA( 4 )
408 COMPLEX*16 CDUMMA( 4 )
409 * ..
410 * .. External Functions ..
411 DOUBLE PRECISION DLAMCH
412 EXTERNAL DLAMCH
413 * ..
414 * .. External Subroutines ..
415 EXTERNAL DLABAD, DLAFTS, DLASUM, XERBLA, ZCOPY, ZGEHRD,
416 $ ZGEMM, ZGET10, ZGET22, ZHSEIN, ZHSEQR, ZHST01,
417 $ ZLACPY, ZLASET, ZLATME, ZLATMR, ZLATMS, ZTREVC,
418 $ ZUNGHR, ZUNMHR
419 * ..
420 * .. Intrinsic Functions ..
421 INTRINSIC ABS, DBLE, MAX, MIN, SQRT
422 * ..
423 * .. Data statements ..
424 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
425 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
426 $ 3, 1, 2, 3 /
427 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
428 $ 1, 5, 5, 5, 4, 3, 1 /
429 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
430 * ..
431 * .. Executable Statements ..
432 *
433 * Check for errors
434 *
435 NTESTT = 0
436 INFO = 0
437 *
438 BADNN = .FALSE.
439 NMAX = 0
440 DO 10 J = 1, NSIZES
441 NMAX = MAX( NMAX, NN( J ) )
442 IF( NN( J ).LT.0 )
443 $ BADNN = .TRUE.
444 10 CONTINUE
445 *
446 * Check for errors
447 *
448 IF( NSIZES.LT.0 ) THEN
449 INFO = -1
450 ELSE IF( BADNN ) THEN
451 INFO = -2
452 ELSE IF( NTYPES.LT.0 ) THEN
453 INFO = -3
454 ELSE IF( THRESH.LT.ZERO ) THEN
455 INFO = -6
456 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
457 INFO = -9
458 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
459 INFO = -14
460 ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN
461 INFO = -26
462 END IF
463 *
464 IF( INFO.NE.0 ) THEN
465 CALL XERBLA( 'ZCHKHS', -INFO )
466 RETURN
467 END IF
468 *
469 * Quick return if possible
470 *
471 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
472 $ RETURN
473 *
474 * More important constants
475 *
476 UNFL = DLAMCH( 'Safe minimum' )
477 OVFL = DLAMCH( 'Overflow' )
478 CALL DLABAD( UNFL, OVFL )
479 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
480 ULPINV = ONE / ULP
481 RTUNFL = SQRT( UNFL )
482 RTOVFL = SQRT( OVFL )
483 RTULP = SQRT( ULP )
484 RTULPI = ONE / RTULP
485 *
486 * Loop over sizes, types
487 *
488 NERRS = 0
489 NMATS = 0
490 *
491 DO 260 JSIZE = 1, NSIZES
492 N = NN( JSIZE )
493 N1 = MAX( 1, N )
494 ANINV = ONE / DBLE( N1 )
495 *
496 IF( NSIZES.NE.1 ) THEN
497 MTYPES = MIN( MAXTYP, NTYPES )
498 ELSE
499 MTYPES = MIN( MAXTYP+1, NTYPES )
500 END IF
501 *
502 DO 250 JTYPE = 1, MTYPES
503 IF( .NOT.DOTYPE( JTYPE ) )
504 $ GO TO 250
505 NMATS = NMATS + 1
506 NTEST = 0
507 *
508 * Save ISEED in case of an error.
509 *
510 DO 20 J = 1, 4
511 IOLDSD( J ) = ISEED( J )
512 20 CONTINUE
513 *
514 * Initialize RESULT
515 *
516 DO 30 J = 1, 14
517 RESULT( J ) = ZERO
518 30 CONTINUE
519 *
520 * Compute "A"
521 *
522 * Control parameters:
523 *
524 * KMAGN KCONDS KMODE KTYPE
525 * =1 O(1) 1 clustered 1 zero
526 * =2 large large clustered 2 identity
527 * =3 small exponential Jordan
528 * =4 arithmetic diagonal, (w/ eigenvalues)
529 * =5 random log hermitian, w/ eigenvalues
530 * =6 random general, w/ eigenvalues
531 * =7 random diagonal
532 * =8 random hermitian
533 * =9 random general
534 * =10 random triangular
535 *
536 IF( MTYPES.GT.MAXTYP )
537 $ GO TO 100
538 *
539 ITYPE = KTYPE( JTYPE )
540 IMODE = KMODE( JTYPE )
541 *
542 * Compute norm
543 *
544 GO TO ( 40, 50, 60 )KMAGN( JTYPE )
545 *
546 40 CONTINUE
547 ANORM = ONE
548 GO TO 70
549 *
550 50 CONTINUE
551 ANORM = ( RTOVFL*ULP )*ANINV
552 GO TO 70
553 *
554 60 CONTINUE
555 ANORM = RTUNFL*N*ULPINV
556 GO TO 70
557 *
558 70 CONTINUE
559 *
560 CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
561 IINFO = 0
562 COND = ULPINV
563 *
564 * Special Matrices
565 *
566 IF( ITYPE.EQ.1 ) THEN
567 *
568 * Zero
569 *
570 IINFO = 0
571 ELSE IF( ITYPE.EQ.2 ) THEN
572 *
573 * Identity
574 *
575 DO 80 JCOL = 1, N
576 A( JCOL, JCOL ) = ANORM
577 80 CONTINUE
578 *
579 ELSE IF( ITYPE.EQ.3 ) THEN
580 *
581 * Jordan Block
582 *
583 DO 90 JCOL = 1, N
584 A( JCOL, JCOL ) = ANORM
585 IF( JCOL.GT.1 )
586 $ A( JCOL, JCOL-1 ) = ONE
587 90 CONTINUE
588 *
589 ELSE IF( ITYPE.EQ.4 ) THEN
590 *
591 * Diagonal Matrix, [Eigen]values Specified
592 *
593 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, IMODE, COND,
594 $ CONE, 'T', 'N', WORK( N+1 ), 1, ONE,
595 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
596 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
597 *
598 ELSE IF( ITYPE.EQ.5 ) THEN
599 *
600 * Hermitian, eigenvalues specified
601 *
602 CALL ZLATMS( N, N, 'D', ISEED, 'H', RWORK, IMODE, COND,
603 $ ANORM, N, N, 'N', A, LDA, WORK, IINFO )
604 *
605 ELSE IF( ITYPE.EQ.6 ) THEN
606 *
607 * General, eigenvalues specified
608 *
609 IF( KCONDS( JTYPE ).EQ.1 ) THEN
610 CONDS = ONE
611 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
612 CONDS = RTULPI
613 ELSE
614 CONDS = ZERO
615 END IF
616 *
617 CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
618 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
619 $ A, LDA, WORK( N+1 ), IINFO )
620 *
621 ELSE IF( ITYPE.EQ.7 ) THEN
622 *
623 * Diagonal, random eigenvalues
624 *
625 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
626 $ 'T', 'N', WORK( N+1 ), 1, ONE,
627 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
628 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
629 *
630 ELSE IF( ITYPE.EQ.8 ) THEN
631 *
632 * Hermitian, random eigenvalues
633 *
634 CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
635 $ 'T', 'N', WORK( N+1 ), 1, ONE,
636 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
637 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
638 *
639 ELSE IF( ITYPE.EQ.9 ) THEN
640 *
641 * General, random eigenvalues
642 *
643 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
644 $ 'T', 'N', WORK( N+1 ), 1, ONE,
645 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
646 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
647 *
648 ELSE IF( ITYPE.EQ.10 ) THEN
649 *
650 * Triangular, random eigenvalues
651 *
652 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
653 $ 'T', 'N', WORK( N+1 ), 1, ONE,
654 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
655 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
656 *
657 ELSE
658 *
659 IINFO = 1
660 END IF
661 *
662 IF( IINFO.NE.0 ) THEN
663 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
664 $ IOLDSD
665 INFO = ABS( IINFO )
666 RETURN
667 END IF
668 *
669 100 CONTINUE
670 *
671 * Call ZGEHRD to compute H and U, do tests.
672 *
673 CALL ZLACPY( ' ', N, N, A, LDA, H, LDA )
674 NTEST = 1
675 *
676 ILO = 1
677 IHI = N
678 *
679 CALL ZGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ),
680 $ NWORK-N, IINFO )
681 *
682 IF( IINFO.NE.0 ) THEN
683 RESULT( 1 ) = ULPINV
684 WRITE( NOUNIT, FMT = 9999 )'ZGEHRD', IINFO, N, JTYPE,
685 $ IOLDSD
686 INFO = ABS( IINFO )
687 GO TO 240
688 END IF
689 *
690 DO 120 J = 1, N - 1
691 UU( J+1, J ) = CZERO
692 DO 110 I = J + 2, N
693 U( I, J ) = H( I, J )
694 UU( I, J ) = H( I, J )
695 H( I, J ) = CZERO
696 110 CONTINUE
697 120 CONTINUE
698 CALL ZCOPY( N-1, WORK, 1, TAU, 1 )
699 CALL ZUNGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ),
700 $ NWORK-N, IINFO )
701 NTEST = 2
702 *
703 CALL ZHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK,
704 $ NWORK, RWORK, RESULT( 1 ) )
705 *
706 * Call ZHSEQR to compute T1, T2 and Z, do tests.
707 *
708 * Eigenvalues only (W3)
709 *
710 CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA )
711 NTEST = 3
712 RESULT( 3 ) = ULPINV
713 *
714 CALL ZHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, W3, UZ, LDU,
715 $ WORK, NWORK, IINFO )
716 IF( IINFO.NE.0 ) THEN
717 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(E)', IINFO, N, JTYPE,
718 $ IOLDSD
719 IF( IINFO.LE.N+2 ) THEN
720 INFO = ABS( IINFO )
721 GO TO 240
722 END IF
723 END IF
724 *
725 * Eigenvalues (W1) and Full Schur Form (T2)
726 *
727 CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA )
728 *
729 CALL ZHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, W1, UZ, LDU,
730 $ WORK, NWORK, IINFO )
731 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
732 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(S)', IINFO, N, JTYPE,
733 $ IOLDSD
734 INFO = ABS( IINFO )
735 GO TO 240
736 END IF
737 *
738 * Eigenvalues (W1), Schur Form (T1), and Schur Vectors (UZ)
739 *
740 CALL ZLACPY( ' ', N, N, H, LDA, T1, LDA )
741 CALL ZLACPY( ' ', N, N, U, LDU, UZ, LDU )
742 *
743 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, W1, UZ, LDU,
744 $ WORK, NWORK, IINFO )
745 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
746 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(V)', IINFO, N, JTYPE,
747 $ IOLDSD
748 INFO = ABS( IINFO )
749 GO TO 240
750 END IF
751 *
752 * Compute Z = U' UZ
753 *
754 CALL ZGEMM( 'C', 'N', N, N, N, CONE, U, LDU, UZ, LDU, CZERO,
755 $ Z, LDU )
756 NTEST = 8
757 *
758 * Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
759 * and 4: | I - Z Z' | / ( n ulp )
760 *
761 CALL ZHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK,
762 $ NWORK, RWORK, RESULT( 3 ) )
763 *
764 * Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
765 * and 6: | I - UZ (UZ)' | / ( n ulp )
766 *
767 CALL ZHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK,
768 $ NWORK, RWORK, RESULT( 5 ) )
769 *
770 * Do Test 7: | T2 - T1 | / ( |T| n ulp )
771 *
772 CALL ZGET10( N, N, T2, LDA, T1, LDA, WORK, RWORK,
773 $ RESULT( 7 ) )
774 *
775 * Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp )
776 *
777 TEMP1 = ZERO
778 TEMP2 = ZERO
779 DO 130 J = 1, N
780 TEMP1 = MAX( TEMP1, ABS( W1( J ) ), ABS( W3( J ) ) )
781 TEMP2 = MAX( TEMP2, ABS( W1( J )-W3( J ) ) )
782 130 CONTINUE
783 *
784 RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
785 *
786 * Compute the Left and Right Eigenvectors of T
787 *
788 * Compute the Right eigenvector Matrix:
789 *
790 NTEST = 9
791 RESULT( 9 ) = ULPINV
792 *
793 * Select every other eigenvector
794 *
795 DO 140 J = 1, N
796 SELECT( J ) = .FALSE.
797 140 CONTINUE
798 DO 150 J = 1, N, 2
799 SELECT( J ) = .TRUE.
800 150 CONTINUE
801 CALL ZTREVC( 'Right', 'All', SELECT, N, T1, LDA, CDUMMA,
802 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
803 IF( IINFO.NE.0 ) THEN
804 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,A)', IINFO, N,
805 $ JTYPE, IOLDSD
806 INFO = ABS( IINFO )
807 GO TO 240
808 END IF
809 *
810 * Test 9: | TR - RW | / ( |T| |R| ulp )
811 *
812 CALL ZGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, W1,
813 $ WORK, RWORK, DUMMA( 1 ) )
814 RESULT( 9 ) = DUMMA( 1 )
815 IF( DUMMA( 2 ).GT.THRESH ) THEN
816 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTREVC',
817 $ DUMMA( 2 ), N, JTYPE, IOLDSD
818 END IF
819 *
820 * Compute selected right eigenvectors and confirm that
821 * they agree with previous right eigenvectors
822 *
823 CALL ZTREVC( 'Right', 'Some', SELECT, N, T1, LDA, CDUMMA,
824 $ LDU, EVECTL, LDU, N, IN, WORK, RWORK, IINFO )
825 IF( IINFO.NE.0 ) THEN
826 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,S)', IINFO, N,
827 $ JTYPE, IOLDSD
828 INFO = ABS( IINFO )
829 GO TO 240
830 END IF
831 *
832 K = 1
833 MATCH = .TRUE.
834 DO 170 J = 1, N
835 IF( SELECT( J ) ) THEN
836 DO 160 JJ = 1, N
837 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN
838 MATCH = .FALSE.
839 GO TO 180
840 END IF
841 160 CONTINUE
842 K = K + 1
843 END IF
844 170 CONTINUE
845 180 CONTINUE
846 IF( .NOT.MATCH )
847 $ WRITE( NOUNIT, FMT = 9997 )'Right', 'ZTREVC', N, JTYPE,
848 $ IOLDSD
849 *
850 * Compute the Left eigenvector Matrix:
851 *
852 NTEST = 10
853 RESULT( 10 ) = ULPINV
854 CALL ZTREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU,
855 $ CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
856 IF( IINFO.NE.0 ) THEN
857 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,A)', IINFO, N,
858 $ JTYPE, IOLDSD
859 INFO = ABS( IINFO )
860 GO TO 240
861 END IF
862 *
863 * Test 10: | LT - WL | / ( |T| |L| ulp )
864 *
865 CALL ZGET22( 'C', 'N', 'C', N, T1, LDA, EVECTL, LDU, W1,
866 $ WORK, RWORK, DUMMA( 3 ) )
867 RESULT( 10 ) = DUMMA( 3 )
868 IF( DUMMA( 4 ).GT.THRESH ) THEN
869 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTREVC', DUMMA( 4 ),
870 $ N, JTYPE, IOLDSD
871 END IF
872 *
873 * Compute selected left eigenvectors and confirm that
874 * they agree with previous left eigenvectors
875 *
876 CALL ZTREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR,
877 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
878 IF( IINFO.NE.0 ) THEN
879 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,S)', IINFO, N,
880 $ JTYPE, IOLDSD
881 INFO = ABS( IINFO )
882 GO TO 240
883 END IF
884 *
885 K = 1
886 MATCH = .TRUE.
887 DO 200 J = 1, N
888 IF( SELECT( J ) ) THEN
889 DO 190 JJ = 1, N
890 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN
891 MATCH = .FALSE.
892 GO TO 210
893 END IF
894 190 CONTINUE
895 K = K + 1
896 END IF
897 200 CONTINUE
898 210 CONTINUE
899 IF( .NOT.MATCH )
900 $ WRITE( NOUNIT, FMT = 9997 )'Left', 'ZTREVC', N, JTYPE,
901 $ IOLDSD
902 *
903 * Call ZHSEIN for Right eigenvectors of H, do test 11
904 *
905 NTEST = 11
906 RESULT( 11 ) = ULPINV
907 DO 220 J = 1, N
908 SELECT( J ) = .TRUE.
909 220 CONTINUE
910 *
911 CALL ZHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
912 $ CDUMMA, LDU, EVECTX, LDU, N1, IN, WORK, RWORK,
913 $ IWORK, IWORK, IINFO )
914 IF( IINFO.NE.0 ) THEN
915 WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(R)', IINFO, N, JTYPE,
916 $ IOLDSD
917 INFO = ABS( IINFO )
918 IF( IINFO.LT.0 )
919 $ GO TO 240
920 ELSE
921 *
922 * Test 11: | HX - XW | / ( |H| |X| ulp )
923 *
924 * (from inverse iteration)
925 *
926 CALL ZGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, W3,
927 $ WORK, RWORK, DUMMA( 1 ) )
928 IF( DUMMA( 1 ).LT.ULPINV )
929 $ RESULT( 11 ) = DUMMA( 1 )*ANINV
930 IF( DUMMA( 2 ).GT.THRESH ) THEN
931 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZHSEIN',
932 $ DUMMA( 2 ), N, JTYPE, IOLDSD
933 END IF
934 END IF
935 *
936 * Call ZHSEIN for Left eigenvectors of H, do test 12
937 *
938 NTEST = 12
939 RESULT( 12 ) = ULPINV
940 DO 230 J = 1, N
941 SELECT( J ) = .TRUE.
942 230 CONTINUE
943 *
944 CALL ZHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
945 $ EVECTY, LDU, CDUMMA, LDU, N1, IN, WORK, RWORK,
946 $ IWORK, IWORK, IINFO )
947 IF( IINFO.NE.0 ) THEN
948 WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(L)', IINFO, N, JTYPE,
949 $ IOLDSD
950 INFO = ABS( IINFO )
951 IF( IINFO.LT.0 )
952 $ GO TO 240
953 ELSE
954 *
955 * Test 12: | YH - WY | / ( |H| |Y| ulp )
956 *
957 * (from inverse iteration)
958 *
959 CALL ZGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, W3,
960 $ WORK, RWORK, DUMMA( 3 ) )
961 IF( DUMMA( 3 ).LT.ULPINV )
962 $ RESULT( 12 ) = DUMMA( 3 )*ANINV
963 IF( DUMMA( 4 ).GT.THRESH ) THEN
964 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZHSEIN',
965 $ DUMMA( 4 ), N, JTYPE, IOLDSD
966 END IF
967 END IF
968 *
969 * Call ZUNMHR for Right eigenvectors of A, do test 13
970 *
971 NTEST = 13
972 RESULT( 13 ) = ULPINV
973 *
974 CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
975 $ LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO )
976 IF( IINFO.NE.0 ) THEN
977 WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE,
978 $ IOLDSD
979 INFO = ABS( IINFO )
980 IF( IINFO.LT.0 )
981 $ GO TO 240
982 ELSE
983 *
984 * Test 13: | AX - XW | / ( |A| |X| ulp )
985 *
986 * (from inverse iteration)
987 *
988 CALL ZGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, W3,
989 $ WORK, RWORK, DUMMA( 1 ) )
990 IF( DUMMA( 1 ).LT.ULPINV )
991 $ RESULT( 13 ) = DUMMA( 1 )*ANINV
992 END IF
993 *
994 * Call ZUNMHR for Left eigenvectors of A, do test 14
995 *
996 NTEST = 14
997 RESULT( 14 ) = ULPINV
998 *
999 CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
1000 $ LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO )
1001 IF( IINFO.NE.0 ) THEN
1002 WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE,
1003 $ IOLDSD
1004 INFO = ABS( IINFO )
1005 IF( IINFO.LT.0 )
1006 $ GO TO 240
1007 ELSE
1008 *
1009 * Test 14: | YA - WY | / ( |A| |Y| ulp )
1010 *
1011 * (from inverse iteration)
1012 *
1013 CALL ZGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, W3,
1014 $ WORK, RWORK, DUMMA( 3 ) )
1015 IF( DUMMA( 3 ).LT.ULPINV )
1016 $ RESULT( 14 ) = DUMMA( 3 )*ANINV
1017 END IF
1018 *
1019 * End of Loop -- Check for RESULT(j) > THRESH
1020 *
1021 240 CONTINUE
1022 *
1023 NTESTT = NTESTT + NTEST
1024 CALL DLAFTS( 'ZHS', N, N, JTYPE, NTEST, RESULT, IOLDSD,
1025 $ THRESH, NOUNIT, NERRS )
1026 *
1027 250 CONTINUE
1028 260 CONTINUE
1029 *
1030 * Summary
1031 *
1032 CALL DLASUM( 'ZHS', NOUNIT, NERRS, NTESTT )
1033 *
1034 RETURN
1035 *
1036 9999 FORMAT( ' ZCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
1037 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
1038 9998 FORMAT( ' ZCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ',
1039 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
1040 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
1041 $ ')' )
1042 9997 FORMAT( ' ZCHKHS: Selected ', A, ' Eigenvectors from ', A,
1043 $ ' do not match other eigenvectors ', 9X, 'N=', I6,
1044 $ ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
1045 *
1046 * End of ZCHKHS
1047 *
1048 END
2 $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, W1,
3 $ W3, EVECTL, EVECTR, EVECTY, EVECTX, UU, TAU,
4 $ WORK, NWORK, RWORK, IWORK, SELECT, RESULT,
5 $ INFO )
6 *
7 * -- LAPACK test routine (version 3.1.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * February 2007
10 *
11 * .. Scalar Arguments ..
12 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
13 DOUBLE PRECISION THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * ), SELECT( * )
17 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
18 DOUBLE PRECISION RESULT( 14 ), RWORK( * )
19 COMPLEX*16 A( LDA, * ), EVECTL( LDU, * ),
20 $ EVECTR( LDU, * ), EVECTX( LDU, * ),
21 $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
22 $ T2( LDA, * ), TAU( * ), U( LDU, * ),
23 $ UU( LDU, * ), UZ( LDU, * ), W1( * ), W3( * ),
24 $ WORK( * ), Z( LDU, * )
25 * ..
26 *
27 * Purpose
28 * =======
29 *
30 * ZCHKHS checks the nonsymmetric eigenvalue problem routines.
31 *
32 * ZGEHRD factors A as U H U' , where ' means conjugate
33 * transpose, H is hessenberg, and U is unitary.
34 *
35 * ZUNGHR generates the unitary matrix U.
36 *
37 * ZUNMHR multiplies a matrix by the unitary matrix U.
38 *
39 * ZHSEQR factors H as Z T Z' , where Z is unitary and T
40 * is upper triangular. It also computes the eigenvalues,
41 * w(1), ..., w(n); we define a diagonal matrix W whose
42 * (diagonal) entries are the eigenvalues.
43 *
44 * ZTREVC computes the left eigenvector matrix L and the
45 * right eigenvector matrix R for the matrix T. The
46 * columns of L are the complex conjugates of the left
47 * eigenvectors of T. The columns of R are the right
48 * eigenvectors of T. L is lower triangular, and R is
49 * upper triangular.
50 *
51 * ZHSEIN computes the left eigenvector matrix Y and the
52 * right eigenvector matrix X for the matrix H. The
53 * columns of Y are the complex conjugates of the left
54 * eigenvectors of H. The columns of X are the right
55 * eigenvectors of H. Y is lower triangular, and X is
56 * upper triangular.
57 *
58 * When ZCHKHS is called, a number of matrix "sizes" ("n's") and a
59 * number of matrix "types" are specified. For each size ("n")
60 * and each type of matrix, one matrix will be generated and used
61 * to test the nonsymmetric eigenroutines. For each matrix, 14
62 * tests will be performed:
63 *
64 * (1) | A - U H U**H | / ( |A| n ulp )
65 *
66 * (2) | I - UU**H | / ( n ulp )
67 *
68 * (3) | H - Z T Z**H | / ( |H| n ulp )
69 *
70 * (4) | I - ZZ**H | / ( n ulp )
71 *
72 * (5) | A - UZ H (UZ)**H | / ( |A| n ulp )
73 *
74 * (6) | I - UZ (UZ)**H | / ( n ulp )
75 *
76 * (7) | T(Z computed) - T(Z not computed) | / ( |T| ulp )
77 *
78 * (8) | W(Z computed) - W(Z not computed) | / ( |W| ulp )
79 *
80 * (9) | TR - RW | / ( |T| |R| ulp )
81 *
82 * (10) | L**H T - W**H L | / ( |T| |L| ulp )
83 *
84 * (11) | HX - XW | / ( |H| |X| ulp )
85 *
86 * (12) | Y**H H - W**H Y | / ( |H| |Y| ulp )
87 *
88 * (13) | AX - XW | / ( |A| |X| ulp )
89 *
90 * (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
91 *
92 * The "sizes" are specified by an array NN(1:NSIZES); the value of
93 * each element NN(j) specifies one size.
94 * The "types" are specified by a logical array DOTYPE( 1:NTYPES );
95 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
96 * Currently, the list of possible types is:
97 *
98 * (1) The zero matrix.
99 * (2) The identity matrix.
100 * (3) A (transposed) Jordan block, with 1's on the diagonal.
101 *
102 * (4) A diagonal matrix with evenly spaced entries
103 * 1, ..., ULP and random complex angles.
104 * (ULP = (first number larger than 1) - 1 )
105 * (5) A diagonal matrix with geometrically spaced entries
106 * 1, ..., ULP and random complex angles.
107 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
108 * and random complex angles.
109 *
110 * (7) Same as (4), but multiplied by SQRT( overflow threshold )
111 * (8) Same as (4), but multiplied by SQRT( underflow threshold )
112 *
113 * (9) A matrix of the form U' T U, where U is unitary and
114 * T has evenly spaced entries 1, ..., ULP with random complex
115 * angles on the diagonal and random O(1) entries in the upper
116 * triangle.
117 *
118 * (10) A matrix of the form U' T U, where U is unitary and
119 * T has geometrically spaced entries 1, ..., ULP with random
120 * complex angles on the diagonal and random O(1) entries in
121 * the upper triangle.
122 *
123 * (11) A matrix of the form U' T U, where U is unitary and
124 * T has "clustered" entries 1, ULP,..., ULP with random
125 * complex angles on the diagonal and random O(1) entries in
126 * the upper triangle.
127 *
128 * (12) A matrix of the form U' T U, where U is unitary and
129 * T has complex eigenvalues randomly chosen from
130 * ULP < |z| < 1 and random O(1) entries in the upper
131 * triangle.
132 *
133 * (13) A matrix of the form X' T X, where X has condition
134 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
135 * with random complex angles on the diagonal and random O(1)
136 * entries in the upper triangle.
137 *
138 * (14) A matrix of the form X' T X, where X has condition
139 * SQRT( ULP ) and T has geometrically spaced entries
140 * 1, ..., ULP with random complex angles on the diagonal
141 * and random O(1) entries in the upper triangle.
142 *
143 * (15) A matrix of the form X' T X, where X has condition
144 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
145 * with random complex angles on the diagonal and random O(1)
146 * entries in the upper triangle.
147 *
148 * (16) A matrix of the form X' T X, where X has condition
149 * SQRT( ULP ) and T has complex eigenvalues randomly chosen
150 * from ULP < |z| < 1 and random O(1) entries in the upper
151 * triangle.
152 *
153 * (17) Same as (16), but multiplied by SQRT( overflow threshold )
154 * (18) Same as (16), but multiplied by SQRT( underflow threshold )
155 *
156 * (19) Nonsymmetric matrix with random entries chosen from |z| < 1
157 * (20) Same as (19), but multiplied by SQRT( overflow threshold )
158 * (21) Same as (19), but multiplied by SQRT( underflow threshold )
159 *
160 * Arguments
161 * ==========
162 *
163 * NSIZES - INTEGER
164 * The number of sizes of matrices to use. If it is zero,
165 * ZCHKHS does nothing. It must be at least zero.
166 * Not modified.
167 *
168 * NN - INTEGER array, dimension (NSIZES)
169 * An array containing the sizes to be used for the matrices.
170 * Zero values will be skipped. The values must be at least
171 * zero.
172 * Not modified.
173 *
174 * NTYPES - INTEGER
175 * The number of elements in DOTYPE. If it is zero, ZCHKHS
176 * does nothing. It must be at least zero. If it is MAXTYP+1
177 * and NSIZES is 1, then an additional type, MAXTYP+1 is
178 * defined, which is to use whatever matrix is in A. This
179 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
180 * DOTYPE(MAXTYP+1) is .TRUE. .
181 * Not modified.
182 *
183 * DOTYPE - LOGICAL array, dimension (NTYPES)
184 * If DOTYPE(j) is .TRUE., then for each size in NN a
185 * matrix of that size and of type j will be generated.
186 * If NTYPES is smaller than the maximum number of types
187 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
188 * MAXTYP will not be generated. If NTYPES is larger
189 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
190 * will be ignored.
191 * Not modified.
192 *
193 * ISEED - INTEGER array, dimension (4)
194 * On entry ISEED specifies the seed of the random number
195 * generator. The array elements should be between 0 and 4095;
196 * if not they will be reduced mod 4096. Also, ISEED(4) must
197 * be odd. The random number generator uses a linear
198 * congruential sequence limited to small integers, and so
199 * should produce machine independent random numbers. The
200 * values of ISEED are changed on exit, and can be used in the
201 * next call to ZCHKHS to continue the same random number
202 * sequence.
203 * Modified.
204 *
205 * THRESH - DOUBLE PRECISION
206 * A test will count as "failed" if the "error", computed as
207 * described above, exceeds THRESH. Note that the error
208 * is scaled to be O(1), so THRESH should be a reasonably
209 * small multiple of 1, e.g., 10 or 100. In particular,
210 * it should not depend on the precision (single vs. double)
211 * or the size of the matrix. It must be at least zero.
212 * Not modified.
213 *
214 * NOUNIT - INTEGER
215 * The FORTRAN unit number for printing out error messages
216 * (e.g., if a routine returns IINFO not equal to 0.)
217 * Not modified.
218 *
219 * A - COMPLEX*16 array, dimension (LDA,max(NN))
220 * Used to hold the matrix whose eigenvalues are to be
221 * computed. On exit, A contains the last matrix actually
222 * used.
223 * Modified.
224 *
225 * LDA - INTEGER
226 * The leading dimension of A, H, T1 and T2. It must be at
227 * least 1 and at least max( NN ).
228 * Not modified.
229 *
230 * H - COMPLEX*16 array, dimension (LDA,max(NN))
231 * The upper hessenberg matrix computed by ZGEHRD. On exit,
232 * H contains the Hessenberg form of the matrix in A.
233 * Modified.
234 *
235 * T1 - COMPLEX*16 array, dimension (LDA,max(NN))
236 * The Schur (="quasi-triangular") matrix computed by ZHSEQR
237 * if Z is computed. On exit, T1 contains the Schur form of
238 * the matrix in A.
239 * Modified.
240 *
241 * T2 - COMPLEX*16 array, dimension (LDA,max(NN))
242 * The Schur matrix computed by ZHSEQR when Z is not computed.
243 * This should be identical to T1.
244 * Modified.
245 *
246 * LDU - INTEGER
247 * The leading dimension of U, Z, UZ and UU. It must be at
248 * least 1 and at least max( NN ).
249 * Not modified.
250 *
251 * U - COMPLEX*16 array, dimension (LDU,max(NN))
252 * The unitary matrix computed by ZGEHRD.
253 * Modified.
254 *
255 * Z - COMPLEX*16 array, dimension (LDU,max(NN))
256 * The unitary matrix computed by ZHSEQR.
257 * Modified.
258 *
259 * UZ - COMPLEX*16 array, dimension (LDU,max(NN))
260 * The product of U times Z.
261 * Modified.
262 *
263 * W1 - COMPLEX*16 array, dimension (max(NN))
264 * The eigenvalues of A, as computed by a full Schur
265 * decomposition H = Z T Z'. On exit, W1 contains the
266 * eigenvalues of the matrix in A.
267 * Modified.
268 *
269 * W3 - COMPLEX*16 array, dimension (max(NN))
270 * The eigenvalues of A, as computed by a partial Schur
271 * decomposition (Z not computed, T only computed as much
272 * as is necessary for determining eigenvalues). On exit,
273 * W3 contains the eigenvalues of the matrix in A, possibly
274 * perturbed by ZHSEIN.
275 * Modified.
276 *
277 * EVECTL - COMPLEX*16 array, dimension (LDU,max(NN))
278 * The conjugate transpose of the (upper triangular) left
279 * eigenvector matrix for the matrix in T1.
280 * Modified.
281 *
282 * EVEZTR - COMPLEX*16 array, dimension (LDU,max(NN))
283 * The (upper triangular) right eigenvector matrix for the
284 * matrix in T1.
285 * Modified.
286 *
287 * EVECTY - COMPLEX*16 array, dimension (LDU,max(NN))
288 * The conjugate transpose of the left eigenvector matrix
289 * for the matrix in H.
290 * Modified.
291 *
292 * EVECTX - COMPLEX*16 array, dimension (LDU,max(NN))
293 * The right eigenvector matrix for the matrix in H.
294 * Modified.
295 *
296 * UU - COMPLEX*16 array, dimension (LDU,max(NN))
297 * Details of the unitary matrix computed by ZGEHRD.
298 * Modified.
299 *
300 * TAU - COMPLEX*16 array, dimension (max(NN))
301 * Further details of the unitary matrix computed by ZGEHRD.
302 * Modified.
303 *
304 * WORK - COMPLEX*16 array, dimension (NWORK)
305 * Workspace.
306 * Modified.
307 *
308 * NWORK - INTEGER
309 * The number of entries in WORK. NWORK >= 4*NN(j)*NN(j) + 2.
310 *
311 * RWORK - DOUBLE PRECISION array, dimension (max(NN))
312 * Workspace. Could be equivalenced to IWORK, but not SELECT.
313 * Modified.
314 *
315 * IWORK - INTEGER array, dimension (max(NN))
316 * Workspace.
317 * Modified.
318 *
319 * SELECT - LOGICAL array, dimension (max(NN))
320 * Workspace. Could be equivalenced to IWORK, but not RWORK.
321 * Modified.
322 *
323 * RESULT - DOUBLE PRECISION array, dimension (14)
324 * The values computed by the fourteen tests described above.
325 * The values are currently limited to 1/ulp, to avoid
326 * overflow.
327 * Modified.
328 *
329 * INFO - INTEGER
330 * If 0, then everything ran OK.
331 * -1: NSIZES < 0
332 * -2: Some NN(j) < 0
333 * -3: NTYPES < 0
334 * -6: THRESH < 0
335 * -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
336 * -14: LDU < 1 or LDU < NMAX.
337 * -26: NWORK too small.
338 * If ZLATMR, CLATMS, or CLATME returns an error code, the
339 * absolute value of it is returned.
340 * If 1, then ZHSEQR could not find all the shifts.
341 * If 2, then the EISPACK code (for small blocks) failed.
342 * If >2, then 30*N iterations were not enough to find an
343 * eigenvalue or to decompose the problem.
344 * Modified.
345 *
346 *-----------------------------------------------------------------------
347 *
348 * Some Local Variables and Parameters:
349 * ---- ----- --------- --- ----------
350 *
351 * ZERO, ONE Real 0 and 1.
352 * MAXTYP The number of types defined.
353 * MTEST The number of tests defined: care must be taken
354 * that (1) the size of RESULT, (2) the number of
355 * tests actually performed, and (3) MTEST agree.
356 * NTEST The number of tests performed on this matrix
357 * so far. This should be less than MTEST, and
358 * equal to it by the last test. It will be less
359 * if any of the routines being tested indicates
360 * that it could not compute the matrices that
361 * would be tested.
362 * NMAX Largest value in NN.
363 * NMATS The number of matrices generated so far.
364 * NERRS The number of tests which have exceeded THRESH
365 * so far (computed by DLAFTS).
366 * COND, CONDS,
367 * IMODE Values to be passed to the matrix generators.
368 * ANORM Norm of A; passed to matrix generators.
369 *
370 * OVFL, UNFL Overflow and underflow thresholds.
371 * ULP, ULPINV Finest relative precision and its inverse.
372 * RTOVFL, RTUNFL,
373 * RTULP, RTULPI Square roots of the previous 4 values.
374 *
375 * The following four arrays decode JTYPE:
376 * KTYPE(j) The general type (1-10) for type "j".
377 * KMODE(j) The MODE value to be passed to the matrix
378 * generator for type "j".
379 * KMAGN(j) The order of magnitude ( O(1),
380 * O(overflow^(1/2) ), O(underflow^(1/2) )
381 * KCONDS(j) Selects whether CONDS is to be 1 or
382 * 1/sqrt(ulp). (0 means irrelevant.)
383 *
384 * =====================================================================
385 *
386 * .. Parameters ..
387 DOUBLE PRECISION ZERO, ONE
388 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
389 COMPLEX*16 CZERO, CONE
390 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
391 $ CONE = ( 1.0D+0, 0.0D+0 ) )
392 INTEGER MAXTYP
393 PARAMETER ( MAXTYP = 21 )
394 * ..
395 * .. Local Scalars ..
396 LOGICAL BADNN, MATCH
397 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
398 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
399 $ NMATS, NMAX, NTEST, NTESTT
400 DOUBLE PRECISION ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
401 $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL
402 * ..
403 * .. Local Arrays ..
404 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
405 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
406 $ KTYPE( MAXTYP )
407 DOUBLE PRECISION DUMMA( 4 )
408 COMPLEX*16 CDUMMA( 4 )
409 * ..
410 * .. External Functions ..
411 DOUBLE PRECISION DLAMCH
412 EXTERNAL DLAMCH
413 * ..
414 * .. External Subroutines ..
415 EXTERNAL DLABAD, DLAFTS, DLASUM, XERBLA, ZCOPY, ZGEHRD,
416 $ ZGEMM, ZGET10, ZGET22, ZHSEIN, ZHSEQR, ZHST01,
417 $ ZLACPY, ZLASET, ZLATME, ZLATMR, ZLATMS, ZTREVC,
418 $ ZUNGHR, ZUNMHR
419 * ..
420 * .. Intrinsic Functions ..
421 INTRINSIC ABS, DBLE, MAX, MIN, SQRT
422 * ..
423 * .. Data statements ..
424 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
425 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
426 $ 3, 1, 2, 3 /
427 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
428 $ 1, 5, 5, 5, 4, 3, 1 /
429 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
430 * ..
431 * .. Executable Statements ..
432 *
433 * Check for errors
434 *
435 NTESTT = 0
436 INFO = 0
437 *
438 BADNN = .FALSE.
439 NMAX = 0
440 DO 10 J = 1, NSIZES
441 NMAX = MAX( NMAX, NN( J ) )
442 IF( NN( J ).LT.0 )
443 $ BADNN = .TRUE.
444 10 CONTINUE
445 *
446 * Check for errors
447 *
448 IF( NSIZES.LT.0 ) THEN
449 INFO = -1
450 ELSE IF( BADNN ) THEN
451 INFO = -2
452 ELSE IF( NTYPES.LT.0 ) THEN
453 INFO = -3
454 ELSE IF( THRESH.LT.ZERO ) THEN
455 INFO = -6
456 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
457 INFO = -9
458 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
459 INFO = -14
460 ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN
461 INFO = -26
462 END IF
463 *
464 IF( INFO.NE.0 ) THEN
465 CALL XERBLA( 'ZCHKHS', -INFO )
466 RETURN
467 END IF
468 *
469 * Quick return if possible
470 *
471 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
472 $ RETURN
473 *
474 * More important constants
475 *
476 UNFL = DLAMCH( 'Safe minimum' )
477 OVFL = DLAMCH( 'Overflow' )
478 CALL DLABAD( UNFL, OVFL )
479 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
480 ULPINV = ONE / ULP
481 RTUNFL = SQRT( UNFL )
482 RTOVFL = SQRT( OVFL )
483 RTULP = SQRT( ULP )
484 RTULPI = ONE / RTULP
485 *
486 * Loop over sizes, types
487 *
488 NERRS = 0
489 NMATS = 0
490 *
491 DO 260 JSIZE = 1, NSIZES
492 N = NN( JSIZE )
493 N1 = MAX( 1, N )
494 ANINV = ONE / DBLE( N1 )
495 *
496 IF( NSIZES.NE.1 ) THEN
497 MTYPES = MIN( MAXTYP, NTYPES )
498 ELSE
499 MTYPES = MIN( MAXTYP+1, NTYPES )
500 END IF
501 *
502 DO 250 JTYPE = 1, MTYPES
503 IF( .NOT.DOTYPE( JTYPE ) )
504 $ GO TO 250
505 NMATS = NMATS + 1
506 NTEST = 0
507 *
508 * Save ISEED in case of an error.
509 *
510 DO 20 J = 1, 4
511 IOLDSD( J ) = ISEED( J )
512 20 CONTINUE
513 *
514 * Initialize RESULT
515 *
516 DO 30 J = 1, 14
517 RESULT( J ) = ZERO
518 30 CONTINUE
519 *
520 * Compute "A"
521 *
522 * Control parameters:
523 *
524 * KMAGN KCONDS KMODE KTYPE
525 * =1 O(1) 1 clustered 1 zero
526 * =2 large large clustered 2 identity
527 * =3 small exponential Jordan
528 * =4 arithmetic diagonal, (w/ eigenvalues)
529 * =5 random log hermitian, w/ eigenvalues
530 * =6 random general, w/ eigenvalues
531 * =7 random diagonal
532 * =8 random hermitian
533 * =9 random general
534 * =10 random triangular
535 *
536 IF( MTYPES.GT.MAXTYP )
537 $ GO TO 100
538 *
539 ITYPE = KTYPE( JTYPE )
540 IMODE = KMODE( JTYPE )
541 *
542 * Compute norm
543 *
544 GO TO ( 40, 50, 60 )KMAGN( JTYPE )
545 *
546 40 CONTINUE
547 ANORM = ONE
548 GO TO 70
549 *
550 50 CONTINUE
551 ANORM = ( RTOVFL*ULP )*ANINV
552 GO TO 70
553 *
554 60 CONTINUE
555 ANORM = RTUNFL*N*ULPINV
556 GO TO 70
557 *
558 70 CONTINUE
559 *
560 CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
561 IINFO = 0
562 COND = ULPINV
563 *
564 * Special Matrices
565 *
566 IF( ITYPE.EQ.1 ) THEN
567 *
568 * Zero
569 *
570 IINFO = 0
571 ELSE IF( ITYPE.EQ.2 ) THEN
572 *
573 * Identity
574 *
575 DO 80 JCOL = 1, N
576 A( JCOL, JCOL ) = ANORM
577 80 CONTINUE
578 *
579 ELSE IF( ITYPE.EQ.3 ) THEN
580 *
581 * Jordan Block
582 *
583 DO 90 JCOL = 1, N
584 A( JCOL, JCOL ) = ANORM
585 IF( JCOL.GT.1 )
586 $ A( JCOL, JCOL-1 ) = ONE
587 90 CONTINUE
588 *
589 ELSE IF( ITYPE.EQ.4 ) THEN
590 *
591 * Diagonal Matrix, [Eigen]values Specified
592 *
593 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, IMODE, COND,
594 $ CONE, 'T', 'N', WORK( N+1 ), 1, ONE,
595 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
596 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
597 *
598 ELSE IF( ITYPE.EQ.5 ) THEN
599 *
600 * Hermitian, eigenvalues specified
601 *
602 CALL ZLATMS( N, N, 'D', ISEED, 'H', RWORK, IMODE, COND,
603 $ ANORM, N, N, 'N', A, LDA, WORK, IINFO )
604 *
605 ELSE IF( ITYPE.EQ.6 ) THEN
606 *
607 * General, eigenvalues specified
608 *
609 IF( KCONDS( JTYPE ).EQ.1 ) THEN
610 CONDS = ONE
611 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
612 CONDS = RTULPI
613 ELSE
614 CONDS = ZERO
615 END IF
616 *
617 CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
618 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
619 $ A, LDA, WORK( N+1 ), IINFO )
620 *
621 ELSE IF( ITYPE.EQ.7 ) THEN
622 *
623 * Diagonal, random eigenvalues
624 *
625 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
626 $ 'T', 'N', WORK( N+1 ), 1, ONE,
627 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
628 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
629 *
630 ELSE IF( ITYPE.EQ.8 ) THEN
631 *
632 * Hermitian, random eigenvalues
633 *
634 CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
635 $ 'T', 'N', WORK( N+1 ), 1, ONE,
636 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
637 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
638 *
639 ELSE IF( ITYPE.EQ.9 ) THEN
640 *
641 * General, random eigenvalues
642 *
643 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
644 $ 'T', 'N', WORK( N+1 ), 1, ONE,
645 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
646 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
647 *
648 ELSE IF( ITYPE.EQ.10 ) THEN
649 *
650 * Triangular, random eigenvalues
651 *
652 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
653 $ 'T', 'N', WORK( N+1 ), 1, ONE,
654 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
655 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
656 *
657 ELSE
658 *
659 IINFO = 1
660 END IF
661 *
662 IF( IINFO.NE.0 ) THEN
663 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
664 $ IOLDSD
665 INFO = ABS( IINFO )
666 RETURN
667 END IF
668 *
669 100 CONTINUE
670 *
671 * Call ZGEHRD to compute H and U, do tests.
672 *
673 CALL ZLACPY( ' ', N, N, A, LDA, H, LDA )
674 NTEST = 1
675 *
676 ILO = 1
677 IHI = N
678 *
679 CALL ZGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ),
680 $ NWORK-N, IINFO )
681 *
682 IF( IINFO.NE.0 ) THEN
683 RESULT( 1 ) = ULPINV
684 WRITE( NOUNIT, FMT = 9999 )'ZGEHRD', IINFO, N, JTYPE,
685 $ IOLDSD
686 INFO = ABS( IINFO )
687 GO TO 240
688 END IF
689 *
690 DO 120 J = 1, N - 1
691 UU( J+1, J ) = CZERO
692 DO 110 I = J + 2, N
693 U( I, J ) = H( I, J )
694 UU( I, J ) = H( I, J )
695 H( I, J ) = CZERO
696 110 CONTINUE
697 120 CONTINUE
698 CALL ZCOPY( N-1, WORK, 1, TAU, 1 )
699 CALL ZUNGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ),
700 $ NWORK-N, IINFO )
701 NTEST = 2
702 *
703 CALL ZHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK,
704 $ NWORK, RWORK, RESULT( 1 ) )
705 *
706 * Call ZHSEQR to compute T1, T2 and Z, do tests.
707 *
708 * Eigenvalues only (W3)
709 *
710 CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA )
711 NTEST = 3
712 RESULT( 3 ) = ULPINV
713 *
714 CALL ZHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, W3, UZ, LDU,
715 $ WORK, NWORK, IINFO )
716 IF( IINFO.NE.0 ) THEN
717 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(E)', IINFO, N, JTYPE,
718 $ IOLDSD
719 IF( IINFO.LE.N+2 ) THEN
720 INFO = ABS( IINFO )
721 GO TO 240
722 END IF
723 END IF
724 *
725 * Eigenvalues (W1) and Full Schur Form (T2)
726 *
727 CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA )
728 *
729 CALL ZHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, W1, UZ, LDU,
730 $ WORK, NWORK, IINFO )
731 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
732 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(S)', IINFO, N, JTYPE,
733 $ IOLDSD
734 INFO = ABS( IINFO )
735 GO TO 240
736 END IF
737 *
738 * Eigenvalues (W1), Schur Form (T1), and Schur Vectors (UZ)
739 *
740 CALL ZLACPY( ' ', N, N, H, LDA, T1, LDA )
741 CALL ZLACPY( ' ', N, N, U, LDU, UZ, LDU )
742 *
743 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, W1, UZ, LDU,
744 $ WORK, NWORK, IINFO )
745 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
746 WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(V)', IINFO, N, JTYPE,
747 $ IOLDSD
748 INFO = ABS( IINFO )
749 GO TO 240
750 END IF
751 *
752 * Compute Z = U' UZ
753 *
754 CALL ZGEMM( 'C', 'N', N, N, N, CONE, U, LDU, UZ, LDU, CZERO,
755 $ Z, LDU )
756 NTEST = 8
757 *
758 * Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
759 * and 4: | I - Z Z' | / ( n ulp )
760 *
761 CALL ZHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK,
762 $ NWORK, RWORK, RESULT( 3 ) )
763 *
764 * Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
765 * and 6: | I - UZ (UZ)' | / ( n ulp )
766 *
767 CALL ZHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK,
768 $ NWORK, RWORK, RESULT( 5 ) )
769 *
770 * Do Test 7: | T2 - T1 | / ( |T| n ulp )
771 *
772 CALL ZGET10( N, N, T2, LDA, T1, LDA, WORK, RWORK,
773 $ RESULT( 7 ) )
774 *
775 * Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp )
776 *
777 TEMP1 = ZERO
778 TEMP2 = ZERO
779 DO 130 J = 1, N
780 TEMP1 = MAX( TEMP1, ABS( W1( J ) ), ABS( W3( J ) ) )
781 TEMP2 = MAX( TEMP2, ABS( W1( J )-W3( J ) ) )
782 130 CONTINUE
783 *
784 RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
785 *
786 * Compute the Left and Right Eigenvectors of T
787 *
788 * Compute the Right eigenvector Matrix:
789 *
790 NTEST = 9
791 RESULT( 9 ) = ULPINV
792 *
793 * Select every other eigenvector
794 *
795 DO 140 J = 1, N
796 SELECT( J ) = .FALSE.
797 140 CONTINUE
798 DO 150 J = 1, N, 2
799 SELECT( J ) = .TRUE.
800 150 CONTINUE
801 CALL ZTREVC( 'Right', 'All', SELECT, N, T1, LDA, CDUMMA,
802 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
803 IF( IINFO.NE.0 ) THEN
804 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,A)', IINFO, N,
805 $ JTYPE, IOLDSD
806 INFO = ABS( IINFO )
807 GO TO 240
808 END IF
809 *
810 * Test 9: | TR - RW | / ( |T| |R| ulp )
811 *
812 CALL ZGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, W1,
813 $ WORK, RWORK, DUMMA( 1 ) )
814 RESULT( 9 ) = DUMMA( 1 )
815 IF( DUMMA( 2 ).GT.THRESH ) THEN
816 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTREVC',
817 $ DUMMA( 2 ), N, JTYPE, IOLDSD
818 END IF
819 *
820 * Compute selected right eigenvectors and confirm that
821 * they agree with previous right eigenvectors
822 *
823 CALL ZTREVC( 'Right', 'Some', SELECT, N, T1, LDA, CDUMMA,
824 $ LDU, EVECTL, LDU, N, IN, WORK, RWORK, IINFO )
825 IF( IINFO.NE.0 ) THEN
826 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,S)', IINFO, N,
827 $ JTYPE, IOLDSD
828 INFO = ABS( IINFO )
829 GO TO 240
830 END IF
831 *
832 K = 1
833 MATCH = .TRUE.
834 DO 170 J = 1, N
835 IF( SELECT( J ) ) THEN
836 DO 160 JJ = 1, N
837 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN
838 MATCH = .FALSE.
839 GO TO 180
840 END IF
841 160 CONTINUE
842 K = K + 1
843 END IF
844 170 CONTINUE
845 180 CONTINUE
846 IF( .NOT.MATCH )
847 $ WRITE( NOUNIT, FMT = 9997 )'Right', 'ZTREVC', N, JTYPE,
848 $ IOLDSD
849 *
850 * Compute the Left eigenvector Matrix:
851 *
852 NTEST = 10
853 RESULT( 10 ) = ULPINV
854 CALL ZTREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU,
855 $ CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
856 IF( IINFO.NE.0 ) THEN
857 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,A)', IINFO, N,
858 $ JTYPE, IOLDSD
859 INFO = ABS( IINFO )
860 GO TO 240
861 END IF
862 *
863 * Test 10: | LT - WL | / ( |T| |L| ulp )
864 *
865 CALL ZGET22( 'C', 'N', 'C', N, T1, LDA, EVECTL, LDU, W1,
866 $ WORK, RWORK, DUMMA( 3 ) )
867 RESULT( 10 ) = DUMMA( 3 )
868 IF( DUMMA( 4 ).GT.THRESH ) THEN
869 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTREVC', DUMMA( 4 ),
870 $ N, JTYPE, IOLDSD
871 END IF
872 *
873 * Compute selected left eigenvectors and confirm that
874 * they agree with previous left eigenvectors
875 *
876 CALL ZTREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR,
877 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
878 IF( IINFO.NE.0 ) THEN
879 WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,S)', IINFO, N,
880 $ JTYPE, IOLDSD
881 INFO = ABS( IINFO )
882 GO TO 240
883 END IF
884 *
885 K = 1
886 MATCH = .TRUE.
887 DO 200 J = 1, N
888 IF( SELECT( J ) ) THEN
889 DO 190 JJ = 1, N
890 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN
891 MATCH = .FALSE.
892 GO TO 210
893 END IF
894 190 CONTINUE
895 K = K + 1
896 END IF
897 200 CONTINUE
898 210 CONTINUE
899 IF( .NOT.MATCH )
900 $ WRITE( NOUNIT, FMT = 9997 )'Left', 'ZTREVC', N, JTYPE,
901 $ IOLDSD
902 *
903 * Call ZHSEIN for Right eigenvectors of H, do test 11
904 *
905 NTEST = 11
906 RESULT( 11 ) = ULPINV
907 DO 220 J = 1, N
908 SELECT( J ) = .TRUE.
909 220 CONTINUE
910 *
911 CALL ZHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
912 $ CDUMMA, LDU, EVECTX, LDU, N1, IN, WORK, RWORK,
913 $ IWORK, IWORK, IINFO )
914 IF( IINFO.NE.0 ) THEN
915 WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(R)', IINFO, N, JTYPE,
916 $ IOLDSD
917 INFO = ABS( IINFO )
918 IF( IINFO.LT.0 )
919 $ GO TO 240
920 ELSE
921 *
922 * Test 11: | HX - XW | / ( |H| |X| ulp )
923 *
924 * (from inverse iteration)
925 *
926 CALL ZGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, W3,
927 $ WORK, RWORK, DUMMA( 1 ) )
928 IF( DUMMA( 1 ).LT.ULPINV )
929 $ RESULT( 11 ) = DUMMA( 1 )*ANINV
930 IF( DUMMA( 2 ).GT.THRESH ) THEN
931 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZHSEIN',
932 $ DUMMA( 2 ), N, JTYPE, IOLDSD
933 END IF
934 END IF
935 *
936 * Call ZHSEIN for Left eigenvectors of H, do test 12
937 *
938 NTEST = 12
939 RESULT( 12 ) = ULPINV
940 DO 230 J = 1, N
941 SELECT( J ) = .TRUE.
942 230 CONTINUE
943 *
944 CALL ZHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
945 $ EVECTY, LDU, CDUMMA, LDU, N1, IN, WORK, RWORK,
946 $ IWORK, IWORK, IINFO )
947 IF( IINFO.NE.0 ) THEN
948 WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(L)', IINFO, N, JTYPE,
949 $ IOLDSD
950 INFO = ABS( IINFO )
951 IF( IINFO.LT.0 )
952 $ GO TO 240
953 ELSE
954 *
955 * Test 12: | YH - WY | / ( |H| |Y| ulp )
956 *
957 * (from inverse iteration)
958 *
959 CALL ZGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, W3,
960 $ WORK, RWORK, DUMMA( 3 ) )
961 IF( DUMMA( 3 ).LT.ULPINV )
962 $ RESULT( 12 ) = DUMMA( 3 )*ANINV
963 IF( DUMMA( 4 ).GT.THRESH ) THEN
964 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZHSEIN',
965 $ DUMMA( 4 ), N, JTYPE, IOLDSD
966 END IF
967 END IF
968 *
969 * Call ZUNMHR for Right eigenvectors of A, do test 13
970 *
971 NTEST = 13
972 RESULT( 13 ) = ULPINV
973 *
974 CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
975 $ LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO )
976 IF( IINFO.NE.0 ) THEN
977 WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE,
978 $ IOLDSD
979 INFO = ABS( IINFO )
980 IF( IINFO.LT.0 )
981 $ GO TO 240
982 ELSE
983 *
984 * Test 13: | AX - XW | / ( |A| |X| ulp )
985 *
986 * (from inverse iteration)
987 *
988 CALL ZGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, W3,
989 $ WORK, RWORK, DUMMA( 1 ) )
990 IF( DUMMA( 1 ).LT.ULPINV )
991 $ RESULT( 13 ) = DUMMA( 1 )*ANINV
992 END IF
993 *
994 * Call ZUNMHR for Left eigenvectors of A, do test 14
995 *
996 NTEST = 14
997 RESULT( 14 ) = ULPINV
998 *
999 CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
1000 $ LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO )
1001 IF( IINFO.NE.0 ) THEN
1002 WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE,
1003 $ IOLDSD
1004 INFO = ABS( IINFO )
1005 IF( IINFO.LT.0 )
1006 $ GO TO 240
1007 ELSE
1008 *
1009 * Test 14: | YA - WY | / ( |A| |Y| ulp )
1010 *
1011 * (from inverse iteration)
1012 *
1013 CALL ZGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, W3,
1014 $ WORK, RWORK, DUMMA( 3 ) )
1015 IF( DUMMA( 3 ).LT.ULPINV )
1016 $ RESULT( 14 ) = DUMMA( 3 )*ANINV
1017 END IF
1018 *
1019 * End of Loop -- Check for RESULT(j) > THRESH
1020 *
1021 240 CONTINUE
1022 *
1023 NTESTT = NTESTT + NTEST
1024 CALL DLAFTS( 'ZHS', N, N, JTYPE, NTEST, RESULT, IOLDSD,
1025 $ THRESH, NOUNIT, NERRS )
1026 *
1027 250 CONTINUE
1028 260 CONTINUE
1029 *
1030 * Summary
1031 *
1032 CALL DLASUM( 'ZHS', NOUNIT, NERRS, NTESTT )
1033 *
1034 RETURN
1035 *
1036 9999 FORMAT( ' ZCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
1037 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
1038 9998 FORMAT( ' ZCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ',
1039 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
1040 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
1041 $ ')' )
1042 9997 FORMAT( ' ZCHKHS: Selected ', A, ' Eigenvectors from ', A,
1043 $ ' do not match other eigenvectors ', 9X, 'N=', I6,
1044 $ ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
1045 *
1046 * End of ZCHKHS
1047 *
1048 END