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