LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zchkst2stg.f
Go to the documentation of this file.
1 *> \brief \b ZCHKST2STG
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 ZCHKST2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
13 * WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
14 * LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT,
15 * INFO )
16 *
17 * .. Scalar Arguments ..
18 * INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT,
19 * $ NSIZES, NTYPES
20 * DOUBLE PRECISION THRESH
21 * ..
22 * .. Array Arguments ..
23 * LOGICAL DOTYPE( * )
24 * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
25 * DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ),
26 * $ RESULT( * ), RWORK( * ), SD( * ), SE( * ),
27 * $ WA1( * ), WA2( * ), WA3( * ), WR( * )
28 * COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ),
29 * $ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZCHKST2STG checks the Hermitian eigenvalue problem routines
39 *> using the 2-stage reduction techniques. Since the generation
40 *> of Q or the vectors is not available in this release, we only
41 *> compare the eigenvalue resulting when using the 2-stage to the
42 *> one considered as reference using the standard 1-stage reduction
43 *> ZHETRD. For that, we call the standard ZHETRD and compute D1 using
44 *> DSTEQR, then we call the 2-stage ZHETRD_2STAGE with Upper and Lower
45 *> and we compute D2 and D3 using DSTEQR and then we replaced tests
46 *> 3 and 4 by tests 11 and 12. test 1 and 2 remain to verify that
47 *> the 1-stage results are OK and can be trusted.
48 *> This testing routine will converge to the ZCHKST in the next
49 *> release when vectors and generation of Q will be implemented.
50 *>
51 *> ZHETRD factors A as U S U* , where * means conjugate transpose,
52 *> S is real symmetric tridiagonal, and U is unitary.
53 *> ZHETRD can use either just the lower or just the upper triangle
54 *> of A; ZCHKST2STG checks both cases.
55 *> U is represented as a product of Householder
56 *> transformations, whose vectors are stored in the first
57 *> n-1 columns of V, and whose scale factors are in TAU.
58 *>
59 *> ZHPTRD does the same as ZHETRD, except that A and V are stored
60 *> in "packed" format.
61 *>
62 *> ZUNGTR constructs the matrix U from the contents of V and TAU.
63 *>
64 *> ZUPGTR constructs the matrix U from the contents of VP and TAU.
65 *>
66 *> ZSTEQR factors S as Z D1 Z* , where Z is the unitary
67 *> matrix of eigenvectors and D1 is a diagonal matrix with
68 *> the eigenvalues on the diagonal. D2 is the matrix of
69 *> eigenvalues computed when Z is not computed.
70 *>
71 *> DSTERF computes D3, the matrix of eigenvalues, by the
72 *> PWK method, which does not yield eigenvectors.
73 *>
74 *> ZPTEQR factors S as Z4 D4 Z4* , for a
75 *> Hermitian positive definite tridiagonal matrix.
76 *> D5 is the matrix of eigenvalues computed when Z is not
77 *> computed.
78 *>
79 *> DSTEBZ computes selected eigenvalues. WA1, WA2, and
80 *> WA3 will denote eigenvalues computed to high
81 *> absolute accuracy, with different range options.
82 *> WR will denote eigenvalues computed to high relative
83 *> accuracy.
84 *>
85 *> ZSTEIN computes Y, the eigenvectors of S, given the
86 *> eigenvalues.
87 *>
88 *> ZSTEDC factors S as Z D1 Z* , where Z is the unitary
89 *> matrix of eigenvectors and D1 is a diagonal matrix with
90 *> the eigenvalues on the diagonal ('I' option). It may also
91 *> update an input unitary matrix, usually the output
92 *> from ZHETRD/ZUNGTR or ZHPTRD/ZUPGTR ('V' option). It may
93 *> also just compute eigenvalues ('N' option).
94 *>
95 *> ZSTEMR factors S as Z D1 Z* , where Z is the unitary
96 *> matrix of eigenvectors and D1 is a diagonal matrix with
97 *> the eigenvalues on the diagonal ('I' option). ZSTEMR
98 *> uses the Relatively Robust Representation whenever possible.
99 *>
100 *> When ZCHKST2STG is called, a number of matrix "sizes" ("n's") and a
101 *> number of matrix "types" are specified. For each size ("n")
102 *> and each type of matrix, one matrix will be generated and used
103 *> to test the Hermitian eigenroutines. For each matrix, a number
104 *> of tests will be performed:
105 *>
106 *> (1) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='U', ... )
107 *>
108 *> (2) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='U', ... )
109 *>
110 *> (3) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='L', ... )
111 *> replaced by | D1 - D2 | / ( |D1| ulp ) where D1 is the
112 *> eigenvalue matrix computed using S and D2 is the
113 *> eigenvalue matrix computed using S_2stage the output of
114 *> ZHETRD_2STAGE("N", "U",....). D1 and D2 are computed
115 *> via DSTEQR('N',...)
116 *>
117 *> (4) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='L', ... )
118 *> replaced by | D1 - D3 | / ( |D1| ulp ) where D1 is the
119 *> eigenvalue matrix computed using S and D3 is the
120 *> eigenvalue matrix computed using S_2stage the output of
121 *> ZHETRD_2STAGE("N", "L",....). D1 and D3 are computed
122 *> via DSTEQR('N',...)
123 *>
124 *> (5-8) Same as 1-4, but for ZHPTRD and ZUPGTR.
125 *>
126 *> (9) | S - Z D Z* | / ( |S| n ulp ) ZSTEQR('V',...)
127 *>
128 *> (10) | I - ZZ* | / ( n ulp ) ZSTEQR('V',...)
129 *>
130 *> (11) | D1 - D2 | / ( |D1| ulp ) ZSTEQR('N',...)
131 *>
132 *> (12) | D1 - D3 | / ( |D1| ulp ) DSTERF
133 *>
134 *> (13) 0 if the true eigenvalues (computed by sturm count)
135 *> of S are within THRESH of
136 *> those in D1. 2*THRESH if they are not. (Tested using
137 *> DSTECH)
138 *>
139 *> For S positive definite,
140 *>
141 *> (14) | S - Z4 D4 Z4* | / ( |S| n ulp ) ZPTEQR('V',...)
142 *>
143 *> (15) | I - Z4 Z4* | / ( n ulp ) ZPTEQR('V',...)
144 *>
145 *> (16) | D4 - D5 | / ( 100 |D4| ulp ) ZPTEQR('N',...)
146 *>
147 *> When S is also diagonally dominant by the factor gamma < 1,
148 *>
149 *> (17) max | D4(i) - WR(i) | / ( |D4(i)| omega ) ,
150 *> i
151 *> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
152 *> DSTEBZ( 'A', 'E', ...)
153 *>
154 *> (18) | WA1 - D3 | / ( |D3| ulp ) DSTEBZ( 'A', 'E', ...)
155 *>
156 *> (19) ( max { min | WA2(i)-WA3(j) | } +
157 *> i j
158 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
159 *> i j
160 *> DSTEBZ( 'I', 'E', ...)
161 *>
162 *> (20) | S - Y WA1 Y* | / ( |S| n ulp ) DSTEBZ, ZSTEIN
163 *>
164 *> (21) | I - Y Y* | / ( n ulp ) DSTEBZ, ZSTEIN
165 *>
166 *> (22) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('I')
167 *>
168 *> (23) | I - ZZ* | / ( n ulp ) ZSTEDC('I')
169 *>
170 *> (24) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('V')
171 *>
172 *> (25) | I - ZZ* | / ( n ulp ) ZSTEDC('V')
173 *>
174 *> (26) | D1 - D2 | / ( |D1| ulp ) ZSTEDC('V') and
175 *> ZSTEDC('N')
176 *>
177 *> Test 27 is disabled at the moment because ZSTEMR does not
178 *> guarantee high relatvie accuracy.
179 *>
180 *> (27) max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
181 *> i
182 *> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
183 *> ZSTEMR('V', 'A')
184 *>
185 *> (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
186 *> i
187 *> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
188 *> ZSTEMR('V', 'I')
189 *>
190 *> Tests 29 through 34 are disable at present because ZSTEMR
191 *> does not handle partial spectrum requests.
192 *>
193 *> (29) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'I')
194 *>
195 *> (30) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'I')
196 *>
197 *> (31) ( max { min | WA2(i)-WA3(j) | } +
198 *> i j
199 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
200 *> i j
201 *> ZSTEMR('N', 'I') vs. CSTEMR('V', 'I')
202 *>
203 *> (32) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'V')
204 *>
205 *> (33) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'V')
206 *>
207 *> (34) ( max { min | WA2(i)-WA3(j) | } +
208 *> i j
209 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
210 *> i j
211 *> ZSTEMR('N', 'V') vs. CSTEMR('V', 'V')
212 *>
213 *> (35) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'A')
214 *>
215 *> (36) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'A')
216 *>
217 *> (37) ( max { min | WA2(i)-WA3(j) | } +
218 *> i j
219 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
220 *> i j
221 *> ZSTEMR('N', 'A') vs. CSTEMR('V', 'A')
222 *>
223 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
224 *> each element NN(j) specifies one size.
225 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
226 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
227 *> Currently, the list of possible types is:
228 *>
229 *> (1) The zero matrix.
230 *> (2) The identity matrix.
231 *>
232 *> (3) A diagonal matrix with evenly spaced entries
233 *> 1, ..., ULP and random signs.
234 *> (ULP = (first number larger than 1) - 1 )
235 *> (4) A diagonal matrix with geometrically spaced entries
236 *> 1, ..., ULP and random signs.
237 *> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
238 *> and random signs.
239 *>
240 *> (6) Same as (4), but multiplied by SQRT( overflow threshold )
241 *> (7) Same as (4), but multiplied by SQRT( underflow threshold )
242 *>
243 *> (8) A matrix of the form U* D U, where U is unitary and
244 *> D has evenly spaced entries 1, ..., ULP with random signs
245 *> on the diagonal.
246 *>
247 *> (9) A matrix of the form U* D U, where U is unitary and
248 *> D has geometrically spaced entries 1, ..., ULP with random
249 *> signs on the diagonal.
250 *>
251 *> (10) A matrix of the form U* D U, where U is unitary and
252 *> D has "clustered" entries 1, ULP,..., ULP with random
253 *> signs on the diagonal.
254 *>
255 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
256 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
257 *>
258 *> (13) Hermitian matrix with random entries chosen from (-1,1).
259 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
260 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
261 *> (16) Same as (8), but diagonal elements are all positive.
262 *> (17) Same as (9), but diagonal elements are all positive.
263 *> (18) Same as (10), but diagonal elements are all positive.
264 *> (19) Same as (16), but multiplied by SQRT( overflow threshold )
265 *> (20) Same as (16), but multiplied by SQRT( underflow threshold )
266 *> (21) A diagonally dominant tridiagonal matrix with geometrically
267 *> spaced diagonal entries 1, ..., ULP.
268 *> \endverbatim
269 *
270 * Arguments:
271 * ==========
272 *
273 *> \param[in] NSIZES
274 *> \verbatim
275 *> NSIZES is INTEGER
276 *> The number of sizes of matrices to use. If it is zero,
277 *> ZCHKST2STG does nothing. It must be at least zero.
278 *> \endverbatim
279 *>
280 *> \param[in] NN
281 *> \verbatim
282 *> NN is INTEGER array, dimension (NSIZES)
283 *> An array containing the sizes to be used for the matrices.
284 *> Zero values will be skipped. The values must be at least
285 *> zero.
286 *> \endverbatim
287 *>
288 *> \param[in] NTYPES
289 *> \verbatim
290 *> NTYPES is INTEGER
291 *> The number of elements in DOTYPE. If it is zero, ZCHKST2STG
292 *> does nothing. It must be at least zero. If it is MAXTYP+1
293 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
294 *> defined, which is to use whatever matrix is in A. This
295 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
296 *> DOTYPE(MAXTYP+1) is .TRUE. .
297 *> \endverbatim
298 *>
299 *> \param[in] DOTYPE
300 *> \verbatim
301 *> DOTYPE is LOGICAL array, dimension (NTYPES)
302 *> If DOTYPE(j) is .TRUE., then for each size in NN a
303 *> matrix of that size and of type j will be generated.
304 *> If NTYPES is smaller than the maximum number of types
305 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
306 *> MAXTYP will not be generated. If NTYPES is larger
307 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
308 *> will be ignored.
309 *> \endverbatim
310 *>
311 *> \param[in,out] ISEED
312 *> \verbatim
313 *> ISEED is INTEGER array, dimension (4)
314 *> On entry ISEED specifies the seed of the random number
315 *> generator. The array elements should be between 0 and 4095;
316 *> if not they will be reduced mod 4096. Also, ISEED(4) must
317 *> be odd. The random number generator uses a linear
318 *> congruential sequence limited to small integers, and so
319 *> should produce machine independent random numbers. The
320 *> values of ISEED are changed on exit, and can be used in the
321 *> next call to ZCHKST2STG to continue the same random number
322 *> sequence.
323 *> \endverbatim
324 *>
325 *> \param[in] THRESH
326 *> \verbatim
327 *> THRESH is DOUBLE PRECISION
328 *> A test will count as "failed" if the "error", computed as
329 *> described above, exceeds THRESH. Note that the error
330 *> is scaled to be O(1), so THRESH should be a reasonably
331 *> small multiple of 1, e.g., 10 or 100. In particular,
332 *> it should not depend on the precision (single vs. double)
333 *> or the size of the matrix. It must be at least zero.
334 *> \endverbatim
335 *>
336 *> \param[in] NOUNIT
337 *> \verbatim
338 *> NOUNIT is INTEGER
339 *> The FORTRAN unit number for printing out error messages
340 *> (e.g., if a routine returns IINFO not equal to 0.)
341 *> \endverbatim
342 *>
343 *> \param[in,out] A
344 *> \verbatim
345 *> A is COMPLEX*16 array of
346 *> dimension ( LDA , max(NN) )
347 *> Used to hold the matrix whose eigenvalues are to be
348 *> computed. On exit, A contains the last matrix actually
349 *> used.
350 *> \endverbatim
351 *>
352 *> \param[in] LDA
353 *> \verbatim
354 *> LDA is INTEGER
355 *> The leading dimension of A. It must be at
356 *> least 1 and at least max( NN ).
357 *> \endverbatim
358 *>
359 *> \param[out] AP
360 *> \verbatim
361 *> AP is COMPLEX*16 array of
362 *> dimension( max(NN)*max(NN+1)/2 )
363 *> The matrix A stored in packed format.
364 *> \endverbatim
365 *>
366 *> \param[out] SD
367 *> \verbatim
368 *> SD is DOUBLE PRECISION array of
369 *> dimension( max(NN) )
370 *> The diagonal of the tridiagonal matrix computed by ZHETRD.
371 *> On exit, SD and SE contain the tridiagonal form of the
372 *> matrix in A.
373 *> \endverbatim
374 *>
375 *> \param[out] SE
376 *> \verbatim
377 *> SE is DOUBLE PRECISION array of
378 *> dimension( max(NN) )
379 *> The off-diagonal of the tridiagonal matrix computed by
380 *> ZHETRD. On exit, SD and SE contain the tridiagonal form of
381 *> the matrix in A.
382 *> \endverbatim
383 *>
384 *> \param[out] D1
385 *> \verbatim
386 *> D1 is DOUBLE PRECISION array of
387 *> dimension( max(NN) )
388 *> The eigenvalues of A, as computed by ZSTEQR simlutaneously
389 *> with Z. On exit, the eigenvalues in D1 correspond with the
390 *> matrix in A.
391 *> \endverbatim
392 *>
393 *> \param[out] D2
394 *> \verbatim
395 *> D2 is DOUBLE PRECISION array of
396 *> dimension( max(NN) )
397 *> The eigenvalues of A, as computed by ZSTEQR if Z is not
398 *> computed. On exit, the eigenvalues in D2 correspond with
399 *> the matrix in A.
400 *> \endverbatim
401 *>
402 *> \param[out] D3
403 *> \verbatim
404 *> D3 is DOUBLE PRECISION array of
405 *> dimension( max(NN) )
406 *> The eigenvalues of A, as computed by DSTERF. On exit, the
407 *> eigenvalues in D3 correspond with the matrix in A.
408 *> \endverbatim
409 *>
410 *> \param[out] D4
411 *> \verbatim
412 *> D4 is DOUBLE PRECISION array of
413 *> dimension( max(NN) )
414 *> The eigenvalues of A, as computed by ZPTEQR(V).
415 *> ZPTEQR factors S as Z4 D4 Z4*
416 *> On exit, the eigenvalues in D4 correspond with the matrix in A.
417 *> \endverbatim
418 *>
419 *> \param[out] D5
420 *> \verbatim
421 *> D5 is DOUBLE PRECISION array of
422 *> dimension( max(NN) )
423 *> The eigenvalues of A, as computed by ZPTEQR(N)
424 *> when Z is not computed. On exit, the
425 *> eigenvalues in D4 correspond with the matrix in A.
426 *> \endverbatim
427 *>
428 *> \param[out] WA1
429 *> \verbatim
430 *> WA1 is DOUBLE PRECISION array of
431 *> dimension( max(NN) )
432 *> All eigenvalues of A, computed to high
433 *> absolute accuracy, with different range options.
434 *> as computed by DSTEBZ.
435 *> \endverbatim
436 *>
437 *> \param[out] WA2
438 *> \verbatim
439 *> WA2 is DOUBLE PRECISION array of
440 *> dimension( max(NN) )
441 *> Selected eigenvalues of A, computed to high
442 *> absolute accuracy, with different range options.
443 *> as computed by DSTEBZ.
444 *> Choose random values for IL and IU, and ask for the
445 *> IL-th through IU-th eigenvalues.
446 *> \endverbatim
447 *>
448 *> \param[out] WA3
449 *> \verbatim
450 *> WA3 is DOUBLE PRECISION array of
451 *> dimension( max(NN) )
452 *> Selected eigenvalues of A, computed to high
453 *> absolute accuracy, with different range options.
454 *> as computed by DSTEBZ.
455 *> Determine the values VL and VU of the IL-th and IU-th
456 *> eigenvalues and ask for all eigenvalues in this range.
457 *> \endverbatim
458 *>
459 *> \param[out] WR
460 *> \verbatim
461 *> WR is DOUBLE PRECISION array of
462 *> dimension( max(NN) )
463 *> All eigenvalues of A, computed to high
464 *> absolute accuracy, with different options.
465 *> as computed by DSTEBZ.
466 *> \endverbatim
467 *>
468 *> \param[out] U
469 *> \verbatim
470 *> U is COMPLEX*16 array of
471 *> dimension( LDU, max(NN) ).
472 *> The unitary matrix computed by ZHETRD + ZUNGTR.
473 *> \endverbatim
474 *>
475 *> \param[in] LDU
476 *> \verbatim
477 *> LDU is INTEGER
478 *> The leading dimension of U, Z, and V. It must be at least 1
479 *> and at least max( NN ).
480 *> \endverbatim
481 *>
482 *> \param[out] V
483 *> \verbatim
484 *> V is COMPLEX*16 array of
485 *> dimension( LDU, max(NN) ).
486 *> The Housholder vectors computed by ZHETRD in reducing A to
487 *> tridiagonal form. The vectors computed with UPLO='U' are
488 *> in the upper triangle, and the vectors computed with UPLO='L'
489 *> are in the lower triangle. (As described in ZHETRD, the
490 *> sub- and superdiagonal are not set to 1, although the
491 *> true Householder vector has a 1 in that position. The
492 *> routines that use V, such as ZUNGTR, set those entries to
493 *> 1 before using them, and then restore them later.)
494 *> \endverbatim
495 *>
496 *> \param[out] VP
497 *> \verbatim
498 *> VP is COMPLEX*16 array of
499 *> dimension( max(NN)*max(NN+1)/2 )
500 *> The matrix V stored in packed format.
501 *> \endverbatim
502 *>
503 *> \param[out] TAU
504 *> \verbatim
505 *> TAU is COMPLEX*16 array of
506 *> dimension( max(NN) )
507 *> The Householder factors computed by ZHETRD in reducing A
508 *> to tridiagonal form.
509 *> \endverbatim
510 *>
511 *> \param[out] Z
512 *> \verbatim
513 *> Z is COMPLEX*16 array of
514 *> dimension( LDU, max(NN) ).
515 *> The unitary matrix of eigenvectors computed by ZSTEQR,
516 *> ZPTEQR, and ZSTEIN.
517 *> \endverbatim
518 *>
519 *> \param[out] WORK
520 *> \verbatim
521 *> WORK is COMPLEX*16 array of
522 *> dimension( LWORK )
523 *> \endverbatim
524 *>
525 *> \param[in] LWORK
526 *> \verbatim
527 *> LWORK is INTEGER
528 *> The number of entries in WORK. This must be at least
529 *> 1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2
530 *> where Nmax = max( NN(j), 2 ) and lg = log base 2.
531 *> \endverbatim
532 *>
533 *> \param[out] IWORK
534 *> \verbatim
535 *> IWORK is INTEGER array,
536 *> Workspace.
537 *> \endverbatim
538 *>
539 *> \param[out] LIWORK
540 *> \verbatim
541 *> LIWORK is INTEGER
542 *> The number of entries in IWORK. This must be at least
543 *> 6 + 6*Nmax + 5 * Nmax * lg Nmax
544 *> where Nmax = max( NN(j), 2 ) and lg = log base 2.
545 *> \endverbatim
546 *>
547 *> \param[out] RWORK
548 *> \verbatim
549 *> RWORK is DOUBLE PRECISION array
550 *> \endverbatim
551 *>
552 *> \param[in] LRWORK
553 *> \verbatim
554 *> LRWORK is INTEGER
555 *> The number of entries in LRWORK (dimension( ??? )
556 *> \endverbatim
557 *>
558 *> \param[out] RESULT
559 *> \verbatim
560 *> RESULT is DOUBLE PRECISION array, dimension (26)
561 *> The values computed by the tests described above.
562 *> The values are currently limited to 1/ulp, to avoid
563 *> overflow.
564 *> \endverbatim
565 *>
566 *> \param[out] INFO
567 *> \verbatim
568 *> INFO is INTEGER
569 *> If 0, then everything ran OK.
570 *> -1: NSIZES < 0
571 *> -2: Some NN(j) < 0
572 *> -3: NTYPES < 0
573 *> -5: THRESH < 0
574 *> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
575 *> -23: LDU < 1 or LDU < NMAX.
576 *> -29: LWORK too small.
577 *> If ZLATMR, CLATMS, ZHETRD, ZUNGTR, ZSTEQR, DSTERF,
578 *> or ZUNMC2 returns an error code, the
579 *> absolute value of it is returned.
580 *>
581 *>-----------------------------------------------------------------------
582 *>
583 *> Some Local Variables and Parameters:
584 *> ---- ----- --------- --- ----------
585 *> ZERO, ONE Real 0 and 1.
586 *> MAXTYP The number of types defined.
587 *> NTEST The number of tests performed, or which can
588 *> be performed so far, for the current matrix.
589 *> NTESTT The total number of tests performed so far.
590 *> NBLOCK Blocksize as returned by ENVIR.
591 *> NMAX Largest value in NN.
592 *> NMATS The number of matrices generated so far.
593 *> NERRS The number of tests which have exceeded THRESH
594 *> so far.
595 *> COND, IMODE Values to be passed to the matrix generators.
596 *> ANORM Norm of A; passed to matrix generators.
597 *>
598 *> OVFL, UNFL Overflow and underflow thresholds.
599 *> ULP, ULPINV Finest relative precision and its inverse.
600 *> RTOVFL, RTUNFL Square roots of the previous 2 values.
601 *> The following four arrays decode JTYPE:
602 *> KTYPE(j) The general type (1-10) for type "j".
603 *> KMODE(j) The MODE value to be passed to the matrix
604 *> generator for type "j".
605 *> KMAGN(j) The order of magnitude ( O(1),
606 *> O(overflow^(1/2) ), O(underflow^(1/2) )
607 *> \endverbatim
608 *
609 * Authors:
610 * ========
611 *
612 *> \author Univ. of Tennessee
613 *> \author Univ. of California Berkeley
614 *> \author Univ. of Colorado Denver
615 *> \author NAG Ltd.
616 *
617 *> \ingroup complex16_eig
618 *
619 * =====================================================================
620  SUBROUTINE zchkst2stg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
621  $ NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
622  $ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
623  $ LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT,
624  $ INFO )
625 *
626 * -- LAPACK test routine --
627 * -- LAPACK is a software package provided by Univ. of Tennessee, --
628 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
629 *
630 * .. Scalar Arguments ..
631  INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT,
632  $ NSIZES, NTYPES
633  DOUBLE PRECISION THRESH
634 * ..
635 * .. Array Arguments ..
636  LOGICAL DOTYPE( * )
637  INTEGER ISEED( 4 ), IWORK( * ), NN( * )
638  DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ),
639  $ RESULT( * ), RWORK( * ), SD( * ), SE( * ),
640  $ wa1( * ), wa2( * ), wa3( * ), wr( * )
641  COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ),
642  $ v( ldu, * ), vp( * ), work( * ), z( ldu, * )
643 * ..
644 *
645 * =====================================================================
646 *
647 * .. Parameters ..
648  DOUBLE PRECISION ZERO, ONE, TWO, EIGHT, TEN, HUN
649  PARAMETER ( ZERO = 0.0d0, one = 1.0d0, two = 2.0d0,
650  $ eight = 8.0d0, ten = 10.0d0, hun = 100.0d0 )
651  COMPLEX*16 CZERO, CONE
652  parameter( czero = ( 0.0d+0, 0.0d+0 ),
653  $ cone = ( 1.0d+0, 0.0d+0 ) )
654  DOUBLE PRECISION HALF
655  parameter( half = one / two )
656  INTEGER MAXTYP
657  PARAMETER ( MAXTYP = 21 )
658  LOGICAL CRANGE
659  parameter( crange = .false. )
660  LOGICAL CREL
661  parameter( crel = .false. )
662 * ..
663 * .. Local Scalars ..
664  LOGICAL BADNN, TRYRAC
665  INTEGER I, IINFO, IL, IMODE, INDE, INDRWK, ITEMP,
666  $ ITYPE, IU, J, JC, JR, JSIZE, JTYPE, LGN,
667  $ LIWEDC, LOG2UI, LRWEDC, LWEDC, M, M2, M3,
668  $ mtypes, n, nap, nblock, nerrs, nmats, nmax,
669  $ nsplit, ntest, ntestt, lh, lw
670  DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
671  $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
672  $ ULPINV, UNFL, VL, VU
673 * ..
674 * .. Local Arrays ..
675  INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
676  $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
677  $ KTYPE( MAXTYP )
678  DOUBLE PRECISION DUMMA( 1 )
679 * ..
680 * .. External Functions ..
681  INTEGER ILAENV
682  DOUBLE PRECISION DLAMCH, DLARND, DSXT1
683  EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1
684 * ..
685 * .. External Subroutines ..
686  EXTERNAL dcopy, dlabad, dlasum, dstebz, dstech, dsterf,
691 * ..
692 * .. Intrinsic Functions ..
693  INTRINSIC abs, dble, dconjg, int, log, max, min, sqrt
694 * ..
695 * .. Data statements ..
696  DATA ktype / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
697  $ 8, 8, 9, 9, 9, 9, 9, 10 /
698  DATA kmagn / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
699  $ 2, 3, 1, 1, 1, 2, 3, 1 /
700  DATA kmode / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
701  $ 0, 0, 4, 3, 1, 4, 4, 3 /
702 * ..
703 * .. Executable Statements ..
704 *
705 * Keep ftnchek happy
706  idumma( 1 ) = 1
707 *
708 * Check for errors
709 *
710  ntestt = 0
711  info = 0
712 *
713 * Important constants
714 *
715  badnn = .false.
716  tryrac = .true.
717  nmax = 1
718  DO 10 j = 1, nsizes
719  nmax = max( nmax, nn( j ) )
720  IF( nn( j ).LT.0 )
721  $ badnn = .true.
722  10 CONTINUE
723 *
724  nblock = ilaenv( 1, 'ZHETRD', 'L', nmax, -1, -1, -1 )
725  nblock = min( nmax, max( 1, nblock ) )
726 *
727 * Check for errors
728 *
729  IF( nsizes.LT.0 ) THEN
730  info = -1
731  ELSE IF( badnn ) THEN
732  info = -2
733  ELSE IF( ntypes.LT.0 ) THEN
734  info = -3
735  ELSE IF( lda.LT.nmax ) THEN
736  info = -9
737  ELSE IF( ldu.LT.nmax ) THEN
738  info = -23
739  ELSE IF( 2*max( 2, nmax )**2.GT.lwork ) THEN
740  info = -29
741  END IF
742 *
743  IF( info.NE.0 ) THEN
744  CALL xerbla( 'ZCHKST2STG', -info )
745  RETURN
746  END IF
747 *
748 * Quick return if possible
749 *
750  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
751  $ RETURN
752 *
753 * More Important constants
754 *
755  unfl = dlamch( 'Safe minimum' )
756  ovfl = one / unfl
757  CALL dlabad( unfl, ovfl )
758  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
759  ulpinv = one / ulp
760  log2ui = int( log( ulpinv ) / log( two ) )
761  rtunfl = sqrt( unfl )
762  rtovfl = sqrt( ovfl )
763 *
764 * Loop over sizes, types
765 *
766  DO 20 i = 1, 4
767  iseed2( i ) = iseed( i )
768  20 CONTINUE
769  nerrs = 0
770  nmats = 0
771 *
772  DO 310 jsize = 1, nsizes
773  n = nn( jsize )
774  IF( n.GT.0 ) THEN
775  lgn = int( log( dble( n ) ) / log( two ) )
776  IF( 2**lgn.LT.n )
777  $ lgn = lgn + 1
778  IF( 2**lgn.LT.n )
779  $ lgn = lgn + 1
780  lwedc = 1 + 4*n + 2*n*lgn + 4*n**2
781  lrwedc = 1 + 3*n + 2*n*lgn + 4*n**2
782  liwedc = 6 + 6*n + 5*n*lgn
783  ELSE
784  lwedc = 8
785  lrwedc = 7
786  liwedc = 12
787  END IF
788  nap = ( n*( n+1 ) ) / 2
789  aninv = one / dble( max( 1, n ) )
790 *
791  IF( nsizes.NE.1 ) THEN
792  mtypes = min( maxtyp, ntypes )
793  ELSE
794  mtypes = min( maxtyp+1, ntypes )
795  END IF
796 *
797  DO 300 jtype = 1, mtypes
798  IF( .NOT.dotype( jtype ) )
799  $ GO TO 300
800  nmats = nmats + 1
801  ntest = 0
802 *
803  DO 30 j = 1, 4
804  ioldsd( j ) = iseed( j )
805  30 CONTINUE
806 *
807 * Compute "A"
808 *
809 * Control parameters:
810 *
811 * KMAGN KMODE KTYPE
812 * =1 O(1) clustered 1 zero
813 * =2 large clustered 2 identity
814 * =3 small exponential (none)
815 * =4 arithmetic diagonal, (w/ eigenvalues)
816 * =5 random log Hermitian, w/ eigenvalues
817 * =6 random (none)
818 * =7 random diagonal
819 * =8 random Hermitian
820 * =9 positive definite
821 * =10 diagonally dominant tridiagonal
822 *
823  IF( mtypes.GT.maxtyp )
824  $ GO TO 100
825 *
826  itype = ktype( jtype )
827  imode = kmode( jtype )
828 *
829 * Compute norm
830 *
831  GO TO ( 40, 50, 60 )kmagn( jtype )
832 *
833  40 CONTINUE
834  anorm = one
835  GO TO 70
836 *
837  50 CONTINUE
838  anorm = ( rtovfl*ulp )*aninv
839  GO TO 70
840 *
841  60 CONTINUE
842  anorm = rtunfl*n*ulpinv
843  GO TO 70
844 *
845  70 CONTINUE
846 *
847  CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
848  iinfo = 0
849  IF( jtype.LE.15 ) THEN
850  cond = ulpinv
851  ELSE
852  cond = ulpinv*aninv / ten
853  END IF
854 *
855 * Special Matrices -- Identity & Jordan block
856 *
857 * Zero
858 *
859  IF( itype.EQ.1 ) THEN
860  iinfo = 0
861 *
862  ELSE IF( itype.EQ.2 ) THEN
863 *
864 * Identity
865 *
866  DO 80 jc = 1, n
867  a( jc, jc ) = anorm
868  80 CONTINUE
869 *
870  ELSE IF( itype.EQ.4 ) THEN
871 *
872 * Diagonal Matrix, [Eigen]values Specified
873 *
874  CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
875  $ anorm, 0, 0, 'N', a, lda, work, iinfo )
876 *
877 *
878  ELSE IF( itype.EQ.5 ) THEN
879 *
880 * Hermitian, eigenvalues specified
881 *
882  CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
883  $ anorm, n, n, 'N', a, lda, work, iinfo )
884 *
885  ELSE IF( itype.EQ.7 ) THEN
886 *
887 * Diagonal, random eigenvalues
888 *
889  CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one, cone,
890  $ 'T', 'N', work( n+1 ), 1, one,
891  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
892  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
893 *
894  ELSE IF( itype.EQ.8 ) THEN
895 *
896 * Hermitian, random eigenvalues
897 *
898  CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one, cone,
899  $ 'T', 'N', work( n+1 ), 1, one,
900  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
901  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
902 *
903  ELSE IF( itype.EQ.9 ) THEN
904 *
905 * Positive definite, eigenvalues specified.
906 *
907  CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode, cond,
908  $ anorm, n, n, 'N', a, lda, work, iinfo )
909 *
910  ELSE IF( itype.EQ.10 ) THEN
911 *
912 * Positive definite tridiagonal, eigenvalues specified.
913 *
914  CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode, cond,
915  $ anorm, 1, 1, 'N', a, lda, work, iinfo )
916  DO 90 i = 2, n
917  temp1 = abs( a( i-1, i ) )
918  temp2 = sqrt( abs( a( i-1, i-1 )*a( i, i ) ) )
919  IF( temp1.GT.half*temp2 ) THEN
920  a( i-1, i ) = a( i-1, i )*
921  $ ( half*temp2 / ( unfl+temp1 ) )
922  a( i, i-1 ) = dconjg( a( i-1, i ) )
923  END IF
924  90 CONTINUE
925 *
926  ELSE
927 *
928  iinfo = 1
929  END IF
930 *
931  IF( iinfo.NE.0 ) THEN
932  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
933  $ ioldsd
934  info = abs( iinfo )
935  RETURN
936  END IF
937 *
938  100 CONTINUE
939 *
940 * Call ZHETRD and ZUNGTR to compute S and U from
941 * upper triangle.
942 *
943  CALL zlacpy( 'U', n, n, a, lda, v, ldu )
944 *
945  ntest = 1
946  CALL zhetrd( 'U', n, v, ldu, sd, se, tau, work, lwork,
947  $ iinfo )
948 *
949  IF( iinfo.NE.0 ) THEN
950  WRITE( nounit, fmt = 9999 )'ZHETRD(U)', iinfo, n, jtype,
951  $ ioldsd
952  info = abs( iinfo )
953  IF( iinfo.LT.0 ) THEN
954  RETURN
955  ELSE
956  result( 1 ) = ulpinv
957  GO TO 280
958  END IF
959  END IF
960 *
961  CALL zlacpy( 'U', n, n, v, ldu, u, ldu )
962 *
963  ntest = 2
964  CALL zungtr( 'U', n, u, ldu, tau, work, lwork, iinfo )
965  IF( iinfo.NE.0 ) THEN
966  WRITE( nounit, fmt = 9999 )'ZUNGTR(U)', iinfo, n, jtype,
967  $ ioldsd
968  info = abs( iinfo )
969  IF( iinfo.LT.0 ) THEN
970  RETURN
971  ELSE
972  result( 2 ) = ulpinv
973  GO TO 280
974  END IF
975  END IF
976 *
977 * Do tests 1 and 2
978 *
979  CALL zhet21( 2, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
980  $ ldu, tau, work, rwork, result( 1 ) )
981  CALL zhet21( 3, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
982  $ ldu, tau, work, rwork, result( 2 ) )
983 *
984 * Compute D1 the eigenvalues resulting from the tridiagonal
985 * form using the standard 1-stage algorithm and use it as a
986 * reference to compare with the 2-stage technique
987 *
988 * Compute D1 from the 1-stage and used as reference for the
989 * 2-stage
990 *
991  CALL dcopy( n, sd, 1, d1, 1 )
992  IF( n.GT.0 )
993  $ CALL dcopy( n-1, se, 1, rwork, 1 )
994 *
995  CALL zsteqr( 'N', n, d1, rwork, work, ldu, rwork( n+1 ),
996  $ iinfo )
997  IF( iinfo.NE.0 ) THEN
998  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n, jtype,
999  $ ioldsd
1000  info = abs( iinfo )
1001  IF( iinfo.LT.0 ) THEN
1002  RETURN
1003  ELSE
1004  result( 3 ) = ulpinv
1005  GO TO 280
1006  END IF
1007  END IF
1008 *
1009 * 2-STAGE TRD Upper case is used to compute D2.
1010 * Note to set SD and SE to zero to be sure not reusing
1011 * the one from above. Compare it with D1 computed
1012 * using the 1-stage.
1013 *
1014  CALL dlaset( 'Full', n, 1, zero, zero, sd, n )
1015  CALL dlaset( 'Full', n, 1, zero, zero, se, n )
1016  CALL zlacpy( 'U', n, n, a, lda, v, ldu )
1017  lh = max(1, 4*n)
1018  lw = lwork - lh
1019  CALL zhetrd_2stage( 'N', "U", n, v, ldu, sd, se, tau,
1020  $ work, lh, work( lh+1 ), lw, iinfo )
1021 *
1022 * Compute D2 from the 2-stage Upper case
1023 *
1024  CALL dcopy( n, sd, 1, d2, 1 )
1025  IF( n.GT.0 )
1026  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1027 *
1028  ntest = 3
1029  CALL zsteqr( 'N', n, d2, rwork, work, ldu, rwork( n+1 ),
1030  $ iinfo )
1031  IF( iinfo.NE.0 ) THEN
1032  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n, jtype,
1033  $ ioldsd
1034  info = abs( iinfo )
1035  IF( iinfo.LT.0 ) THEN
1036  RETURN
1037  ELSE
1038  result( 3 ) = ulpinv
1039  GO TO 280
1040  END IF
1041  END IF
1042 *
1043 * 2-STAGE TRD Lower case is used to compute D3.
1044 * Note to set SD and SE to zero to be sure not reusing
1045 * the one from above. Compare it with D1 computed
1046 * using the 1-stage.
1047 *
1048  CALL dlaset( 'Full', n, 1, zero, zero, sd, n )
1049  CALL dlaset( 'Full', n, 1, zero, zero, se, n )
1050  CALL zlacpy( 'L', n, n, a, lda, v, ldu )
1051  CALL zhetrd_2stage( 'N', "L", n, v, ldu, sd, se, tau,
1052  $ work, lh, work( lh+1 ), lw, iinfo )
1053 *
1054 * Compute D3 from the 2-stage Upper case
1055 *
1056  CALL dcopy( n, sd, 1, d3, 1 )
1057  IF( n.GT.0 )
1058  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1059 *
1060  ntest = 4
1061  CALL zsteqr( 'N', n, d3, rwork, work, ldu, rwork( n+1 ),
1062  $ iinfo )
1063  IF( iinfo.NE.0 ) THEN
1064  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n, jtype,
1065  $ ioldsd
1066  info = abs( iinfo )
1067  IF( iinfo.LT.0 ) THEN
1068  RETURN
1069  ELSE
1070  result( 4 ) = ulpinv
1071  GO TO 280
1072  END IF
1073  END IF
1074 *
1075 *
1076 * Do Tests 3 and 4 which are similar to 11 and 12 but with the
1077 * D1 computed using the standard 1-stage reduction as reference
1078 *
1079  ntest = 4
1080  temp1 = zero
1081  temp2 = zero
1082  temp3 = zero
1083  temp4 = zero
1084 *
1085  DO 151 j = 1, n
1086  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1087  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1088  temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
1089  temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
1090  151 CONTINUE
1091 *
1092  result( 3 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1093  result( 4 ) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
1094 *
1095 * Store the upper triangle of A in AP
1096 *
1097  i = 0
1098  DO 120 jc = 1, n
1099  DO 110 jr = 1, jc
1100  i = i + 1
1101  ap( i ) = a( jr, jc )
1102  110 CONTINUE
1103  120 CONTINUE
1104 *
1105 * Call ZHPTRD and ZUPGTR to compute S and U from AP
1106 *
1107  CALL zcopy( nap, ap, 1, vp, 1 )
1108 *
1109  ntest = 5
1110  CALL zhptrd( 'U', n, vp, sd, se, tau, iinfo )
1111 *
1112  IF( iinfo.NE.0 ) THEN
1113  WRITE( nounit, fmt = 9999 )'ZHPTRD(U)', iinfo, n, jtype,
1114  $ ioldsd
1115  info = abs( iinfo )
1116  IF( iinfo.LT.0 ) THEN
1117  RETURN
1118  ELSE
1119  result( 5 ) = ulpinv
1120  GO TO 280
1121  END IF
1122  END IF
1123 *
1124  ntest = 6
1125  CALL zupgtr( 'U', n, vp, tau, u, ldu, work, iinfo )
1126  IF( iinfo.NE.0 ) THEN
1127  WRITE( nounit, fmt = 9999 )'ZUPGTR(U)', iinfo, n, jtype,
1128  $ ioldsd
1129  info = abs( iinfo )
1130  IF( iinfo.LT.0 ) THEN
1131  RETURN
1132  ELSE
1133  result( 6 ) = ulpinv
1134  GO TO 280
1135  END IF
1136  END IF
1137 *
1138 * Do tests 5 and 6
1139 *
1140  CALL zhpt21( 2, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1141  $ work, rwork, result( 5 ) )
1142  CALL zhpt21( 3, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1143  $ work, rwork, result( 6 ) )
1144 *
1145 * Store the lower triangle of A in AP
1146 *
1147  i = 0
1148  DO 140 jc = 1, n
1149  DO 130 jr = jc, n
1150  i = i + 1
1151  ap( i ) = a( jr, jc )
1152  130 CONTINUE
1153  140 CONTINUE
1154 *
1155 * Call ZHPTRD and ZUPGTR to compute S and U from AP
1156 *
1157  CALL zcopy( nap, ap, 1, vp, 1 )
1158 *
1159  ntest = 7
1160  CALL zhptrd( 'L', n, vp, sd, se, tau, iinfo )
1161 *
1162  IF( iinfo.NE.0 ) THEN
1163  WRITE( nounit, fmt = 9999 )'ZHPTRD(L)', iinfo, n, jtype,
1164  $ ioldsd
1165  info = abs( iinfo )
1166  IF( iinfo.LT.0 ) THEN
1167  RETURN
1168  ELSE
1169  result( 7 ) = ulpinv
1170  GO TO 280
1171  END IF
1172  END IF
1173 *
1174  ntest = 8
1175  CALL zupgtr( 'L', n, vp, tau, u, ldu, work, iinfo )
1176  IF( iinfo.NE.0 ) THEN
1177  WRITE( nounit, fmt = 9999 )'ZUPGTR(L)', iinfo, n, jtype,
1178  $ ioldsd
1179  info = abs( iinfo )
1180  IF( iinfo.LT.0 ) THEN
1181  RETURN
1182  ELSE
1183  result( 8 ) = ulpinv
1184  GO TO 280
1185  END IF
1186  END IF
1187 *
1188  CALL zhpt21( 2, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1189  $ work, rwork, result( 7 ) )
1190  CALL zhpt21( 3, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1191  $ work, rwork, result( 8 ) )
1192 *
1193 * Call ZSTEQR to compute D1, D2, and Z, do tests.
1194 *
1195 * Compute D1 and Z
1196 *
1197  CALL dcopy( n, sd, 1, d1, 1 )
1198  IF( n.GT.0 )
1199  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1200  CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1201 *
1202  ntest = 9
1203  CALL zsteqr( 'V', n, d1, rwork, z, ldu, rwork( n+1 ),
1204  $ iinfo )
1205  IF( iinfo.NE.0 ) THEN
1206  WRITE( nounit, fmt = 9999 )'ZSTEQR(V)', iinfo, n, jtype,
1207  $ ioldsd
1208  info = abs( iinfo )
1209  IF( iinfo.LT.0 ) THEN
1210  RETURN
1211  ELSE
1212  result( 9 ) = ulpinv
1213  GO TO 280
1214  END IF
1215  END IF
1216 *
1217 * Compute D2
1218 *
1219  CALL dcopy( n, sd, 1, d2, 1 )
1220  IF( n.GT.0 )
1221  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1222 *
1223  ntest = 11
1224  CALL zsteqr( 'N', n, d2, rwork, work, ldu, rwork( n+1 ),
1225  $ iinfo )
1226  IF( iinfo.NE.0 ) THEN
1227  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n, jtype,
1228  $ ioldsd
1229  info = abs( iinfo )
1230  IF( iinfo.LT.0 ) THEN
1231  RETURN
1232  ELSE
1233  result( 11 ) = ulpinv
1234  GO TO 280
1235  END IF
1236  END IF
1237 *
1238 * Compute D3 (using PWK method)
1239 *
1240  CALL dcopy( n, sd, 1, d3, 1 )
1241  IF( n.GT.0 )
1242  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1243 *
1244  ntest = 12
1245  CALL dsterf( n, d3, rwork, iinfo )
1246  IF( iinfo.NE.0 ) THEN
1247  WRITE( nounit, fmt = 9999 )'DSTERF', iinfo, n, jtype,
1248  $ ioldsd
1249  info = abs( iinfo )
1250  IF( iinfo.LT.0 ) THEN
1251  RETURN
1252  ELSE
1253  result( 12 ) = ulpinv
1254  GO TO 280
1255  END IF
1256  END IF
1257 *
1258 * Do Tests 9 and 10
1259 *
1260  CALL zstt21( n, 0, sd, se, d1, dumma, z, ldu, work, rwork,
1261  $ result( 9 ) )
1262 *
1263 * Do Tests 11 and 12
1264 *
1265  temp1 = zero
1266  temp2 = zero
1267  temp3 = zero
1268  temp4 = zero
1269 *
1270  DO 150 j = 1, n
1271  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1272  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1273  temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
1274  temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
1275  150 CONTINUE
1276 *
1277  result( 11 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1278  result( 12 ) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
1279 *
1280 * Do Test 13 -- Sturm Sequence Test of Eigenvalues
1281 * Go up by factors of two until it succeeds
1282 *
1283  ntest = 13
1284  temp1 = thresh*( half-ulp )
1285 *
1286  DO 160 j = 0, log2ui
1287  CALL dstech( n, sd, se, d1, temp1, rwork, iinfo )
1288  IF( iinfo.EQ.0 )
1289  $ GO TO 170
1290  temp1 = temp1*two
1291  160 CONTINUE
1292 *
1293  170 CONTINUE
1294  result( 13 ) = temp1
1295 *
1296 * For positive definite matrices ( JTYPE.GT.15 ) call ZPTEQR
1297 * and do tests 14, 15, and 16 .
1298 *
1299  IF( jtype.GT.15 ) THEN
1300 *
1301 * Compute D4 and Z4
1302 *
1303  CALL dcopy( n, sd, 1, d4, 1 )
1304  IF( n.GT.0 )
1305  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1306  CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1307 *
1308  ntest = 14
1309  CALL zpteqr( 'V', n, d4, rwork, z, ldu, rwork( n+1 ),
1310  $ iinfo )
1311  IF( iinfo.NE.0 ) THEN
1312  WRITE( nounit, fmt = 9999 )'ZPTEQR(V)', iinfo, n,
1313  $ jtype, ioldsd
1314  info = abs( iinfo )
1315  IF( iinfo.LT.0 ) THEN
1316  RETURN
1317  ELSE
1318  result( 14 ) = ulpinv
1319  GO TO 280
1320  END IF
1321  END IF
1322 *
1323 * Do Tests 14 and 15
1324 *
1325  CALL zstt21( n, 0, sd, se, d4, dumma, z, ldu, work,
1326  $ rwork, result( 14 ) )
1327 *
1328 * Compute D5
1329 *
1330  CALL dcopy( n, sd, 1, d5, 1 )
1331  IF( n.GT.0 )
1332  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1333 *
1334  ntest = 16
1335  CALL zpteqr( 'N', n, d5, rwork, z, ldu, rwork( n+1 ),
1336  $ iinfo )
1337  IF( iinfo.NE.0 ) THEN
1338  WRITE( nounit, fmt = 9999 )'ZPTEQR(N)', iinfo, n,
1339  $ jtype, ioldsd
1340  info = abs( iinfo )
1341  IF( iinfo.LT.0 ) THEN
1342  RETURN
1343  ELSE
1344  result( 16 ) = ulpinv
1345  GO TO 280
1346  END IF
1347  END IF
1348 *
1349 * Do Test 16
1350 *
1351  temp1 = zero
1352  temp2 = zero
1353  DO 180 j = 1, n
1354  temp1 = max( temp1, abs( d4( j ) ), abs( d5( j ) ) )
1355  temp2 = max( temp2, abs( d4( j )-d5( j ) ) )
1356  180 CONTINUE
1357 *
1358  result( 16 ) = temp2 / max( unfl,
1359  $ hun*ulp*max( temp1, temp2 ) )
1360  ELSE
1361  result( 14 ) = zero
1362  result( 15 ) = zero
1363  result( 16 ) = zero
1364  END IF
1365 *
1366 * Call DSTEBZ with different options and do tests 17-18.
1367 *
1368 * If S is positive definite and diagonally dominant,
1369 * ask for all eigenvalues with high relative accuracy.
1370 *
1371  vl = zero
1372  vu = zero
1373  il = 0
1374  iu = 0
1375  IF( jtype.EQ.21 ) THEN
1376  ntest = 17
1377  abstol = unfl + unfl
1378  CALL dstebz( 'A', 'E', n, vl, vu, il, iu, abstol, sd, se,
1379  $ m, nsplit, wr, iwork( 1 ), iwork( n+1 ),
1380  $ rwork, iwork( 2*n+1 ), iinfo )
1381  IF( iinfo.NE.0 ) THEN
1382  WRITE( nounit, fmt = 9999 )'DSTEBZ(A,rel)', iinfo, n,
1383  $ jtype, ioldsd
1384  info = abs( iinfo )
1385  IF( iinfo.LT.0 ) THEN
1386  RETURN
1387  ELSE
1388  result( 17 ) = ulpinv
1389  GO TO 280
1390  END IF
1391  END IF
1392 *
1393 * Do test 17
1394 *
1395  temp2 = two*( two*n-one )*ulp*( one+eight*half**2 ) /
1396  $ ( one-half )**4
1397 *
1398  temp1 = zero
1399  DO 190 j = 1, n
1400  temp1 = max( temp1, abs( d4( j )-wr( n-j+1 ) ) /
1401  $ ( abstol+abs( d4( j ) ) ) )
1402  190 CONTINUE
1403 *
1404  result( 17 ) = temp1 / temp2
1405  ELSE
1406  result( 17 ) = zero
1407  END IF
1408 *
1409 * Now ask for all eigenvalues with high absolute accuracy.
1410 *
1411  ntest = 18
1412  abstol = unfl + unfl
1413  CALL dstebz( 'A', 'E', n, vl, vu, il, iu, abstol, sd, se, m,
1414  $ nsplit, wa1, iwork( 1 ), iwork( n+1 ), rwork,
1415  $ iwork( 2*n+1 ), iinfo )
1416  IF( iinfo.NE.0 ) THEN
1417  WRITE( nounit, fmt = 9999 )'DSTEBZ(A)', iinfo, n, jtype,
1418  $ ioldsd
1419  info = abs( iinfo )
1420  IF( iinfo.LT.0 ) THEN
1421  RETURN
1422  ELSE
1423  result( 18 ) = ulpinv
1424  GO TO 280
1425  END IF
1426  END IF
1427 *
1428 * Do test 18
1429 *
1430  temp1 = zero
1431  temp2 = zero
1432  DO 200 j = 1, n
1433  temp1 = max( temp1, abs( d3( j ) ), abs( wa1( j ) ) )
1434  temp2 = max( temp2, abs( d3( j )-wa1( j ) ) )
1435  200 CONTINUE
1436 *
1437  result( 18 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1438 *
1439 * Choose random values for IL and IU, and ask for the
1440 * IL-th through IU-th eigenvalues.
1441 *
1442  ntest = 19
1443  IF( n.LE.1 ) THEN
1444  il = 1
1445  iu = n
1446  ELSE
1447  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1448  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1449  IF( iu.LT.il ) THEN
1450  itemp = iu
1451  iu = il
1452  il = itemp
1453  END IF
1454  END IF
1455 *
1456  CALL dstebz( 'I', 'E', n, vl, vu, il, iu, abstol, sd, se,
1457  $ m2, nsplit, wa2, iwork( 1 ), iwork( n+1 ),
1458  $ rwork, iwork( 2*n+1 ), iinfo )
1459  IF( iinfo.NE.0 ) THEN
1460  WRITE( nounit, fmt = 9999 )'DSTEBZ(I)', iinfo, n, jtype,
1461  $ ioldsd
1462  info = abs( iinfo )
1463  IF( iinfo.LT.0 ) THEN
1464  RETURN
1465  ELSE
1466  result( 19 ) = ulpinv
1467  GO TO 280
1468  END IF
1469  END IF
1470 *
1471 * Determine the values VL and VU of the IL-th and IU-th
1472 * eigenvalues and ask for all eigenvalues in this range.
1473 *
1474  IF( n.GT.0 ) THEN
1475  IF( il.NE.1 ) THEN
1476  vl = wa1( il ) - max( half*( wa1( il )-wa1( il-1 ) ),
1477  $ ulp*anorm, two*rtunfl )
1478  ELSE
1479  vl = wa1( 1 ) - max( half*( wa1( n )-wa1( 1 ) ),
1480  $ ulp*anorm, two*rtunfl )
1481  END IF
1482  IF( iu.NE.n ) THEN
1483  vu = wa1( iu ) + max( half*( wa1( iu+1 )-wa1( iu ) ),
1484  $ ulp*anorm, two*rtunfl )
1485  ELSE
1486  vu = wa1( n ) + max( half*( wa1( n )-wa1( 1 ) ),
1487  $ ulp*anorm, two*rtunfl )
1488  END IF
1489  ELSE
1490  vl = zero
1491  vu = one
1492  END IF
1493 *
1494  CALL dstebz( 'V', 'E', n, vl, vu, il, iu, abstol, sd, se,
1495  $ m3, nsplit, wa3, iwork( 1 ), iwork( n+1 ),
1496  $ rwork, iwork( 2*n+1 ), iinfo )
1497  IF( iinfo.NE.0 ) THEN
1498  WRITE( nounit, fmt = 9999 )'DSTEBZ(V)', iinfo, n, jtype,
1499  $ ioldsd
1500  info = abs( iinfo )
1501  IF( iinfo.LT.0 ) THEN
1502  RETURN
1503  ELSE
1504  result( 19 ) = ulpinv
1505  GO TO 280
1506  END IF
1507  END IF
1508 *
1509  IF( m3.EQ.0 .AND. n.NE.0 ) THEN
1510  result( 19 ) = ulpinv
1511  GO TO 280
1512  END IF
1513 *
1514 * Do test 19
1515 *
1516  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1517  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1518  IF( n.GT.0 ) THEN
1519  temp3 = max( abs( wa1( n ) ), abs( wa1( 1 ) ) )
1520  ELSE
1521  temp3 = zero
1522  END IF
1523 *
1524  result( 19 ) = ( temp1+temp2 ) / max( unfl, temp3*ulp )
1525 *
1526 * Call ZSTEIN to compute eigenvectors corresponding to
1527 * eigenvalues in WA1. (First call DSTEBZ again, to make sure
1528 * it returns these eigenvalues in the correct order.)
1529 *
1530  ntest = 21
1531  CALL dstebz( 'A', 'B', n, vl, vu, il, iu, abstol, sd, se, m,
1532  $ nsplit, wa1, iwork( 1 ), iwork( n+1 ), rwork,
1533  $ iwork( 2*n+1 ), iinfo )
1534  IF( iinfo.NE.0 ) THEN
1535  WRITE( nounit, fmt = 9999 )'DSTEBZ(A,B)', iinfo, n,
1536  $ jtype, ioldsd
1537  info = abs( iinfo )
1538  IF( iinfo.LT.0 ) THEN
1539  RETURN
1540  ELSE
1541  result( 20 ) = ulpinv
1542  result( 21 ) = ulpinv
1543  GO TO 280
1544  END IF
1545  END IF
1546 *
1547  CALL zstein( n, sd, se, m, wa1, iwork( 1 ), iwork( n+1 ), z,
1548  $ ldu, rwork, iwork( 2*n+1 ), iwork( 3*n+1 ),
1549  $ iinfo )
1550  IF( iinfo.NE.0 ) THEN
1551  WRITE( nounit, fmt = 9999 )'ZSTEIN', iinfo, n, jtype,
1552  $ ioldsd
1553  info = abs( iinfo )
1554  IF( iinfo.LT.0 ) THEN
1555  RETURN
1556  ELSE
1557  result( 20 ) = ulpinv
1558  result( 21 ) = ulpinv
1559  GO TO 280
1560  END IF
1561  END IF
1562 *
1563 * Do tests 20 and 21
1564 *
1565  CALL zstt21( n, 0, sd, se, wa1, dumma, z, ldu, work, rwork,
1566  $ result( 20 ) )
1567 *
1568 * Call ZSTEDC(I) to compute D1 and Z, do tests.
1569 *
1570 * Compute D1 and Z
1571 *
1572  inde = 1
1573  indrwk = inde + n
1574  CALL dcopy( n, sd, 1, d1, 1 )
1575  IF( n.GT.0 )
1576  $ CALL dcopy( n-1, se, 1, rwork( inde ), 1 )
1577  CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1578 *
1579  ntest = 22
1580  CALL zstedc( 'I', n, d1, rwork( inde ), z, ldu, work, lwedc,
1581  $ rwork( indrwk ), lrwedc, iwork, liwedc, iinfo )
1582  IF( iinfo.NE.0 ) THEN
1583  WRITE( nounit, fmt = 9999 )'ZSTEDC(I)', iinfo, n, jtype,
1584  $ ioldsd
1585  info = abs( iinfo )
1586  IF( iinfo.LT.0 ) THEN
1587  RETURN
1588  ELSE
1589  result( 22 ) = ulpinv
1590  GO TO 280
1591  END IF
1592  END IF
1593 *
1594 * Do Tests 22 and 23
1595 *
1596  CALL zstt21( n, 0, sd, se, d1, dumma, z, ldu, work, rwork,
1597  $ result( 22 ) )
1598 *
1599 * Call ZSTEDC(V) to compute D1 and Z, do tests.
1600 *
1601 * Compute D1 and Z
1602 *
1603  CALL dcopy( n, sd, 1, d1, 1 )
1604  IF( n.GT.0 )
1605  $ CALL dcopy( n-1, se, 1, rwork( inde ), 1 )
1606  CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1607 *
1608  ntest = 24
1609  CALL zstedc( 'V', n, d1, rwork( inde ), z, ldu, work, lwedc,
1610  $ rwork( indrwk ), lrwedc, iwork, liwedc, iinfo )
1611  IF( iinfo.NE.0 ) THEN
1612  WRITE( nounit, fmt = 9999 )'ZSTEDC(V)', iinfo, n, jtype,
1613  $ ioldsd
1614  info = abs( iinfo )
1615  IF( iinfo.LT.0 ) THEN
1616  RETURN
1617  ELSE
1618  result( 24 ) = ulpinv
1619  GO TO 280
1620  END IF
1621  END IF
1622 *
1623 * Do Tests 24 and 25
1624 *
1625  CALL zstt21( n, 0, sd, se, d1, dumma, z, ldu, work, rwork,
1626  $ result( 24 ) )
1627 *
1628 * Call ZSTEDC(N) to compute D2, do tests.
1629 *
1630 * Compute D2
1631 *
1632  CALL dcopy( n, sd, 1, d2, 1 )
1633  IF( n.GT.0 )
1634  $ CALL dcopy( n-1, se, 1, rwork( inde ), 1 )
1635  CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1636 *
1637  ntest = 26
1638  CALL zstedc( 'N', n, d2, rwork( inde ), z, ldu, work, lwedc,
1639  $ rwork( indrwk ), lrwedc, iwork, liwedc, iinfo )
1640  IF( iinfo.NE.0 ) THEN
1641  WRITE( nounit, fmt = 9999 )'ZSTEDC(N)', iinfo, n, jtype,
1642  $ ioldsd
1643  info = abs( iinfo )
1644  IF( iinfo.LT.0 ) THEN
1645  RETURN
1646  ELSE
1647  result( 26 ) = ulpinv
1648  GO TO 280
1649  END IF
1650  END IF
1651 *
1652 * Do Test 26
1653 *
1654  temp1 = zero
1655  temp2 = zero
1656 *
1657  DO 210 j = 1, n
1658  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1659  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1660  210 CONTINUE
1661 *
1662  result( 26 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1663 *
1664 * Only test ZSTEMR if IEEE compliant
1665 *
1666  IF( ilaenv( 10, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
1667  $ ilaenv( 11, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
1668 *
1669 * Call ZSTEMR, do test 27 (relative eigenvalue accuracy)
1670 *
1671 * If S is positive definite and diagonally dominant,
1672 * ask for all eigenvalues with high relative accuracy.
1673 *
1674  vl = zero
1675  vu = zero
1676  il = 0
1677  iu = 0
1678  IF( jtype.EQ.21 .AND. crel ) THEN
1679  ntest = 27
1680  abstol = unfl + unfl
1681  CALL zstemr( 'V', 'A', n, sd, se, vl, vu, il, iu,
1682  $ m, wr, z, ldu, n, iwork( 1 ), tryrac,
1683  $ rwork, lrwork, iwork( 2*n+1 ), lwork-2*n,
1684  $ iinfo )
1685  IF( iinfo.NE.0 ) THEN
1686  WRITE( nounit, fmt = 9999 )'ZSTEMR(V,A,rel)',
1687  $ iinfo, n, jtype, ioldsd
1688  info = abs( iinfo )
1689  IF( iinfo.LT.0 ) THEN
1690  RETURN
1691  ELSE
1692  result( 27 ) = ulpinv
1693  GO TO 270
1694  END IF
1695  END IF
1696 *
1697 * Do test 27
1698 *
1699  temp2 = two*( two*n-one )*ulp*( one+eight*half**2 ) /
1700  $ ( one-half )**4
1701 *
1702  temp1 = zero
1703  DO 220 j = 1, n
1704  temp1 = max( temp1, abs( d4( j )-wr( n-j+1 ) ) /
1705  $ ( abstol+abs( d4( j ) ) ) )
1706  220 CONTINUE
1707 *
1708  result( 27 ) = temp1 / temp2
1709 *
1710  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1711  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1712  IF( iu.LT.il ) THEN
1713  itemp = iu
1714  iu = il
1715  il = itemp
1716  END IF
1717 *
1718  IF( crange ) THEN
1719  ntest = 28
1720  abstol = unfl + unfl
1721  CALL zstemr( 'V', 'I', n, sd, se, vl, vu, il, iu,
1722  $ m, wr, z, ldu, n, iwork( 1 ), tryrac,
1723  $ rwork, lrwork, iwork( 2*n+1 ),
1724  $ lwork-2*n, iinfo )
1725 *
1726  IF( iinfo.NE.0 ) THEN
1727  WRITE( nounit, fmt = 9999 )'ZSTEMR(V,I,rel)',
1728  $ iinfo, n, jtype, ioldsd
1729  info = abs( iinfo )
1730  IF( iinfo.LT.0 ) THEN
1731  RETURN
1732  ELSE
1733  result( 28 ) = ulpinv
1734  GO TO 270
1735  END IF
1736  END IF
1737 *
1738 *
1739 * Do test 28
1740 *
1741  temp2 = two*( two*n-one )*ulp*
1742  $ ( one+eight*half**2 ) / ( one-half )**4
1743 *
1744  temp1 = zero
1745  DO 230 j = il, iu
1746  temp1 = max( temp1, abs( wr( j-il+1 )-d4( n-j+
1747  $ 1 ) ) / ( abstol+abs( wr( j-il+1 ) ) ) )
1748  230 CONTINUE
1749 *
1750  result( 28 ) = temp1 / temp2
1751  ELSE
1752  result( 28 ) = zero
1753  END IF
1754  ELSE
1755  result( 27 ) = zero
1756  result( 28 ) = zero
1757  END IF
1758 *
1759 * Call ZSTEMR(V,I) to compute D1 and Z, do tests.
1760 *
1761 * Compute D1 and Z
1762 *
1763  CALL dcopy( n, sd, 1, d5, 1 )
1764  IF( n.GT.0 )
1765  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1766  CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1767 *
1768  IF( crange ) THEN
1769  ntest = 29
1770  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1771  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1772  IF( iu.LT.il ) THEN
1773  itemp = iu
1774  iu = il
1775  il = itemp
1776  END IF
1777  CALL zstemr( 'V', 'I', n, d5, rwork, vl, vu, il, iu,
1778  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1779  $ rwork( n+1 ), lrwork-n, iwork( 2*n+1 ),
1780  $ liwork-2*n, iinfo )
1781  IF( iinfo.NE.0 ) THEN
1782  WRITE( nounit, fmt = 9999 )'ZSTEMR(V,I)', iinfo,
1783  $ n, jtype, ioldsd
1784  info = abs( iinfo )
1785  IF( iinfo.LT.0 ) THEN
1786  RETURN
1787  ELSE
1788  result( 29 ) = ulpinv
1789  GO TO 280
1790  END IF
1791  END IF
1792 *
1793 * Do Tests 29 and 30
1794 *
1795 *
1796 * Call ZSTEMR to compute D2, do tests.
1797 *
1798 * Compute D2
1799 *
1800  CALL dcopy( n, sd, 1, d5, 1 )
1801  IF( n.GT.0 )
1802  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1803 *
1804  ntest = 31
1805  CALL zstemr( 'N', 'I', n, d5, rwork, vl, vu, il, iu,
1806  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1807  $ rwork( n+1 ), lrwork-n, iwork( 2*n+1 ),
1808  $ liwork-2*n, iinfo )
1809  IF( iinfo.NE.0 ) THEN
1810  WRITE( nounit, fmt = 9999 )'ZSTEMR(N,I)', iinfo,
1811  $ n, jtype, ioldsd
1812  info = abs( iinfo )
1813  IF( iinfo.LT.0 ) THEN
1814  RETURN
1815  ELSE
1816  result( 31 ) = ulpinv
1817  GO TO 280
1818  END IF
1819  END IF
1820 *
1821 * Do Test 31
1822 *
1823  temp1 = zero
1824  temp2 = zero
1825 *
1826  DO 240 j = 1, iu - il + 1
1827  temp1 = max( temp1, abs( d1( j ) ),
1828  $ abs( d2( j ) ) )
1829  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1830  240 CONTINUE
1831 *
1832  result( 31 ) = temp2 / max( unfl,
1833  $ ulp*max( temp1, temp2 ) )
1834 *
1835 *
1836 * Call ZSTEMR(V,V) to compute D1 and Z, do tests.
1837 *
1838 * Compute D1 and Z
1839 *
1840  CALL dcopy( n, sd, 1, d5, 1 )
1841  IF( n.GT.0 )
1842  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1843  CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1844 *
1845  ntest = 32
1846 *
1847  IF( n.GT.0 ) THEN
1848  IF( il.NE.1 ) THEN
1849  vl = d2( il ) - max( half*
1850  $ ( d2( il )-d2( il-1 ) ), ulp*anorm,
1851  $ two*rtunfl )
1852  ELSE
1853  vl = d2( 1 ) - max( half*( d2( n )-d2( 1 ) ),
1854  $ ulp*anorm, two*rtunfl )
1855  END IF
1856  IF( iu.NE.n ) THEN
1857  vu = d2( iu ) + max( half*
1858  $ ( d2( iu+1 )-d2( iu ) ), ulp*anorm,
1859  $ two*rtunfl )
1860  ELSE
1861  vu = d2( n ) + max( half*( d2( n )-d2( 1 ) ),
1862  $ ulp*anorm, two*rtunfl )
1863  END IF
1864  ELSE
1865  vl = zero
1866  vu = one
1867  END IF
1868 *
1869  CALL zstemr( 'V', 'V', n, d5, rwork, vl, vu, il, iu,
1870  $ m, d1, z, ldu, m, iwork( 1 ), tryrac,
1871  $ rwork( n+1 ), lrwork-n, iwork( 2*n+1 ),
1872  $ liwork-2*n, iinfo )
1873  IF( iinfo.NE.0 ) THEN
1874  WRITE( nounit, fmt = 9999 )'ZSTEMR(V,V)', iinfo,
1875  $ n, jtype, ioldsd
1876  info = abs( iinfo )
1877  IF( iinfo.LT.0 ) THEN
1878  RETURN
1879  ELSE
1880  result( 32 ) = ulpinv
1881  GO TO 280
1882  END IF
1883  END IF
1884 *
1885 * Do Tests 32 and 33
1886 *
1887  CALL zstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1888  $ m, rwork, result( 32 ) )
1889 *
1890 * Call ZSTEMR to compute D2, do tests.
1891 *
1892 * Compute D2
1893 *
1894  CALL dcopy( n, sd, 1, d5, 1 )
1895  IF( n.GT.0 )
1896  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1897 *
1898  ntest = 34
1899  CALL zstemr( 'N', 'V', n, d5, rwork, vl, vu, il, iu,
1900  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1901  $ rwork( n+1 ), lrwork-n, iwork( 2*n+1 ),
1902  $ liwork-2*n, iinfo )
1903  IF( iinfo.NE.0 ) THEN
1904  WRITE( nounit, fmt = 9999 )'ZSTEMR(N,V)', iinfo,
1905  $ n, jtype, ioldsd
1906  info = abs( iinfo )
1907  IF( iinfo.LT.0 ) THEN
1908  RETURN
1909  ELSE
1910  result( 34 ) = ulpinv
1911  GO TO 280
1912  END IF
1913  END IF
1914 *
1915 * Do Test 34
1916 *
1917  temp1 = zero
1918  temp2 = zero
1919 *
1920  DO 250 j = 1, iu - il + 1
1921  temp1 = max( temp1, abs( d1( j ) ),
1922  $ abs( d2( j ) ) )
1923  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1924  250 CONTINUE
1925 *
1926  result( 34 ) = temp2 / max( unfl,
1927  $ ulp*max( temp1, temp2 ) )
1928  ELSE
1929  result( 29 ) = zero
1930  result( 30 ) = zero
1931  result( 31 ) = zero
1932  result( 32 ) = zero
1933  result( 33 ) = zero
1934  result( 34 ) = zero
1935  END IF
1936 *
1937 *
1938 * Call ZSTEMR(V,A) to compute D1 and Z, do tests.
1939 *
1940 * Compute D1 and Z
1941 *
1942  CALL dcopy( n, sd, 1, d5, 1 )
1943  IF( n.GT.0 )
1944  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1945 *
1946  ntest = 35
1947 *
1948  CALL zstemr( 'V', 'A', n, d5, rwork, vl, vu, il, iu,
1949  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1950  $ rwork( n+1 ), lrwork-n, iwork( 2*n+1 ),
1951  $ liwork-2*n, iinfo )
1952  IF( iinfo.NE.0 ) THEN
1953  WRITE( nounit, fmt = 9999 )'ZSTEMR(V,A)', iinfo, n,
1954  $ jtype, ioldsd
1955  info = abs( iinfo )
1956  IF( iinfo.LT.0 ) THEN
1957  RETURN
1958  ELSE
1959  result( 35 ) = ulpinv
1960  GO TO 280
1961  END IF
1962  END IF
1963 *
1964 * Do Tests 35 and 36
1965 *
1966  CALL zstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work, m,
1967  $ rwork, result( 35 ) )
1968 *
1969 * Call ZSTEMR to compute D2, do tests.
1970 *
1971 * Compute D2
1972 *
1973  CALL dcopy( n, sd, 1, d5, 1 )
1974  IF( n.GT.0 )
1975  $ CALL dcopy( n-1, se, 1, rwork, 1 )
1976 *
1977  ntest = 37
1978  CALL zstemr( 'N', 'A', n, d5, rwork, vl, vu, il, iu,
1979  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1980  $ rwork( n+1 ), lrwork-n, iwork( 2*n+1 ),
1981  $ liwork-2*n, iinfo )
1982  IF( iinfo.NE.0 ) THEN
1983  WRITE( nounit, fmt = 9999 )'ZSTEMR(N,A)', iinfo, n,
1984  $ jtype, ioldsd
1985  info = abs( iinfo )
1986  IF( iinfo.LT.0 ) THEN
1987  RETURN
1988  ELSE
1989  result( 37 ) = ulpinv
1990  GO TO 280
1991  END IF
1992  END IF
1993 *
1994 * Do Test 34
1995 *
1996  temp1 = zero
1997  temp2 = zero
1998 *
1999  DO 260 j = 1, n
2000  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
2001  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
2002  260 CONTINUE
2003 *
2004  result( 37 ) = temp2 / max( unfl,
2005  $ ulp*max( temp1, temp2 ) )
2006  END IF
2007  270 CONTINUE
2008  280 CONTINUE
2009  ntestt = ntestt + ntest
2010 *
2011 * End of Loop -- Check for RESULT(j) > THRESH
2012 *
2013 *
2014 * Print out tests which fail.
2015 *
2016  DO 290 jr = 1, ntest
2017  IF( result( jr ).GE.thresh ) THEN
2018 *
2019 * If this is the first test to fail,
2020 * print a header to the data file.
2021 *
2022  IF( nerrs.EQ.0 ) THEN
2023  WRITE( nounit, fmt = 9998 )'ZST'
2024  WRITE( nounit, fmt = 9997 )
2025  WRITE( nounit, fmt = 9996 )
2026  WRITE( nounit, fmt = 9995 )'Hermitian'
2027  WRITE( nounit, fmt = 9994 )
2028 *
2029 * Tests performed
2030 *
2031  WRITE( nounit, fmt = 9987 )
2032  END IF
2033  nerrs = nerrs + 1
2034  IF( result( jr ).LT.10000.0d0 ) THEN
2035  WRITE( nounit, fmt = 9989 )n, jtype, ioldsd, jr,
2036  $ result( jr )
2037  ELSE
2038  WRITE( nounit, fmt = 9988 )n, jtype, ioldsd, jr,
2039  $ result( jr )
2040  END IF
2041  END IF
2042  290 CONTINUE
2043  300 CONTINUE
2044  310 CONTINUE
2045 *
2046 * Summary
2047 *
2048  CALL dlasum( 'ZST', nounit, nerrs, ntestt )
2049  RETURN
2050 *
2051  9999 FORMAT( ' ZCHKST2STG: ', a, ' returned INFO=', i6, '.', / 9x,
2052  $ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
2053 *
2054  9998 FORMAT( / 1x, a3, ' -- Complex Hermitian eigenvalue problem' )
2055  9997 FORMAT( ' Matrix types (see ZCHKST2STG for details): ' )
2056 *
2057  9996 FORMAT( / ' Special Matrices:',
2058  $ / ' 1=Zero matrix. ',
2059  $ ' 5=Diagonal: clustered entries.',
2060  $ / ' 2=Identity matrix. ',
2061  $ ' 6=Diagonal: large, evenly spaced.',
2062  $ / ' 3=Diagonal: evenly spaced entries. ',
2063  $ ' 7=Diagonal: small, evenly spaced.',
2064  $ / ' 4=Diagonal: geometr. spaced entries.' )
2065  9995 FORMAT( ' Dense ', a, ' Matrices:',
2066  $ / ' 8=Evenly spaced eigenvals. ',
2067  $ ' 12=Small, evenly spaced eigenvals.',
2068  $ / ' 9=Geometrically spaced eigenvals. ',
2069  $ ' 13=Matrix with random O(1) entries.',
2070  $ / ' 10=Clustered eigenvalues. ',
2071  $ ' 14=Matrix with large random entries.',
2072  $ / ' 11=Large, evenly spaced eigenvals. ',
2073  $ ' 15=Matrix with small random entries.' )
2074  9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues',
2075  $ / ' 17=Positive definite, geometrically spaced eigenvlaues',
2076  $ / ' 18=Positive definite, clustered eigenvalues',
2077  $ / ' 19=Positive definite, small evenly spaced eigenvalues',
2078  $ / ' 20=Positive definite, large evenly spaced eigenvalues',
2079  $ / ' 21=Diagonally dominant tridiagonal, geometrically',
2080  $ ' spaced eigenvalues' )
2081 *
2082  9989 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
2083  $ 4( i4, ',' ), ' result ', i3, ' is', 0p, f8.2 )
2084  9988 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
2085  $ 4( i4, ',' ), ' result ', i3, ' is', 1p, d10.3 )
2086 *
2087  9987 FORMAT( / 'Test performed: see ZCHKST2STG for details.', / )
2088 *
2089 * End of ZCHKST2STG
2090 *
2091  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zstt22(N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, LDWORK, RWORK, RESULT)
ZSTT22
Definition: zstt22.f:145
subroutine zstt21(N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK, RESULT)
ZSTT21
Definition: zstt21.f:133
subroutine zhet21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RWORK, RESULT)
ZHET21
Definition: zhet21.f:214
subroutine zchkst2stg(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5, WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT, INFO)
ZCHKST2STG
Definition: zchkst2stg.f:625
subroutine zhpt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RWORK, RESULT)
ZHPT21
Definition: zhpt21.f:228
subroutine zlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
ZLATMS
Definition: zlatms.f:332
subroutine zlatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
ZLATMR
Definition: zlatmr.f:490
subroutine zhetrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
ZHETRD_2STAGE
subroutine zhetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
ZHETRD
Definition: zhetrd.f:192
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:132
subroutine zstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
ZSTEMR
Definition: zstemr.f:338
subroutine zungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGTR
Definition: zungtr.f:123
subroutine zhptrd(UPLO, N, AP, D, E, TAU, INFO)
ZHPTRD
Definition: zhptrd.f:151
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:182
subroutine zstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
ZSTEDC
Definition: zstedc.f:212
subroutine zupgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
ZUPGTR
Definition: zupgtr.f:114
subroutine zpteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZPTEQR
Definition: zpteqr.f:145
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:43
subroutine dstech(N, A, B, EIG, TOL, WORK, INFO)
DSTECH
Definition: dstech.f:101