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