LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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,
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
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
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:112
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:275
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
DORGTR
Definition: dorgtr.f:125
subroutine dsytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
DSYTRD
Definition: dsytrd.f:194
subroutine dlatmr(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)
DLATMR
Definition: dlatmr.f:473
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dpteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DPTEQR
Definition: dpteqr.f:147
subroutine dchkst(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, IWORK, LIWORK, RESULT, INFO)
DCHKST
Definition: dchkst.f:593
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:133
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dstt22(N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, LDWORK, RESULT)
DSTT22
Definition: dstt22.f:141
subroutine dspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
DSPT21
Definition: dspt21.f:221
subroutine dstech(N, A, B, EIG, TOL, WORK, INFO)
DSTECH
Definition: dstech.f:103
subroutine dsptrd(UPLO, N, AP, D, E, TAU, INFO)
DSPTRD
Definition: dsptrd.f:152
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:45
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:191
subroutine dopgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
DOPGTR
Definition: dopgtr.f:116
subroutine dsyt21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RESULT)
DSYT21
Definition: dsyt21.f:207
subroutine dstt21(N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RESULT)
DSTT21
Definition: dstt21.f:129
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:323
subroutine dstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEMR
Definition: dstemr.f:323
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:176