LAPACK  3.10.0
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 spectrum 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 *> \ingroup double_eig
585 *
586 * =====================================================================
587  SUBROUTINE dchkst( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
588  $ NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
589  $ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
590  $ LWORK, IWORK, LIWORK, RESULT, INFO )
591 *
592 * -- LAPACK test routine --
593 * -- LAPACK is a software package provided by Univ. of Tennessee, --
594 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
595 *
596 * .. Scalar Arguments ..
597  INTEGER INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
598  $ NTYPES
599  DOUBLE PRECISION THRESH
600 * ..
601 * .. Array Arguments ..
602  LOGICAL DOTYPE( * )
603  INTEGER ISEED( 4 ), IWORK( * ), NN( * )
604  DOUBLE PRECISION A( LDA, * ), AP( * ), D1( * ), D2( * ),
605  $ d3( * ), d4( * ), d5( * ), result( * ),
606  $ sd( * ), se( * ), tau( * ), u( ldu, * ),
607  $ v( ldu, * ), vp( * ), wa1( * ), wa2( * ),
608  $ wa3( * ), work( * ), wr( * ), z( ldu, * )
609 * ..
610 *
611 * =====================================================================
612 *
613 * .. Parameters ..
614  DOUBLE PRECISION ZERO, ONE, TWO, EIGHT, TEN, HUN
615  PARAMETER ( ZERO = 0.0d0, one = 1.0d0, two = 2.0d0,
616  $ eight = 8.0d0, ten = 10.0d0, hun = 100.0d0 )
617  DOUBLE PRECISION HALF
618  parameter( half = one / two )
619  INTEGER MAXTYP
620  parameter( maxtyp = 21 )
621  LOGICAL SRANGE
622  parameter( srange = .false. )
623  LOGICAL SREL
624  parameter( srel = .false. )
625 * ..
626 * .. Local Scalars ..
627  LOGICAL BADNN, TRYRAC
628  INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC,
629  $ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC,
630  $ m, m2, m3, mtypes, n, nap, nblock, nerrs,
631  $ nmats, nmax, nsplit, ntest, ntestt
632  DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
633  $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
634  $ ULPINV, UNFL, VL, VU
635 * ..
636 * .. Local Arrays ..
637  INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
638  $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
639  $ KTYPE( MAXTYP )
640  DOUBLE PRECISION DUMMA( 1 )
641 * ..
642 * .. External Functions ..
643  INTEGER ILAENV
644  DOUBLE PRECISION DLAMCH, DLARND, DSXT1
645  EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1
646 * ..
647 * .. External Subroutines ..
648  EXTERNAL dcopy, dlabad, dlacpy, dlaset, dlasum, dlatmr,
652 * ..
653 * .. Intrinsic Functions ..
654  INTRINSIC abs, dble, int, log, max, min, sqrt
655 * ..
656 * .. Data statements ..
657  DATA ktype / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
658  $ 8, 8, 9, 9, 9, 9, 9, 10 /
659  DATA kmagn / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
660  $ 2, 3, 1, 1, 1, 2, 3, 1 /
661  DATA kmode / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
662  $ 0, 0, 4, 3, 1, 4, 4, 3 /
663 * ..
664 * .. Executable Statements ..
665 *
666 * Keep ftnchek happy
667  idumma( 1 ) = 1
668 *
669 * Check for errors
670 *
671  ntestt = 0
672  info = 0
673 *
674 * Important constants
675 *
676  badnn = .false.
677  tryrac = .true.
678  nmax = 1
679  DO 10 j = 1, nsizes
680  nmax = max( nmax, nn( j ) )
681  IF( nn( j ).LT.0 )
682  $ badnn = .true.
683  10 CONTINUE
684 *
685  nblock = ilaenv( 1, 'DSYTRD', 'L', nmax, -1, -1, -1 )
686  nblock = min( nmax, max( 1, nblock ) )
687 *
688 * Check for errors
689 *
690  IF( nsizes.LT.0 ) THEN
691  info = -1
692  ELSE IF( badnn ) THEN
693  info = -2
694  ELSE IF( ntypes.LT.0 ) THEN
695  info = -3
696  ELSE IF( lda.LT.nmax ) THEN
697  info = -9
698  ELSE IF( ldu.LT.nmax ) THEN
699  info = -23
700  ELSE IF( 2*max( 2, nmax )**2.GT.lwork ) THEN
701  info = -29
702  END IF
703 *
704  IF( info.NE.0 ) THEN
705  CALL xerbla( 'DCHKST', -info )
706  RETURN
707  END IF
708 *
709 * Quick return if possible
710 *
711  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
712  $ RETURN
713 *
714 * More Important constants
715 *
716  unfl = dlamch( 'Safe minimum' )
717  ovfl = one / unfl
718  CALL dlabad( unfl, ovfl )
719  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
720  ulpinv = one / ulp
721  log2ui = int( log( ulpinv ) / log( two ) )
722  rtunfl = sqrt( unfl )
723  rtovfl = sqrt( ovfl )
724 *
725 * Loop over sizes, types
726 *
727  DO 20 i = 1, 4
728  iseed2( i ) = iseed( i )
729  20 CONTINUE
730  nerrs = 0
731  nmats = 0
732 *
733  DO 310 jsize = 1, nsizes
734  n = nn( jsize )
735  IF( n.GT.0 ) THEN
736  lgn = int( log( dble( n ) ) / log( two ) )
737  IF( 2**lgn.LT.n )
738  $ lgn = lgn + 1
739  IF( 2**lgn.LT.n )
740  $ lgn = lgn + 1
741  lwedc = 1 + 4*n + 2*n*lgn + 4*n**2
742  liwedc = 6 + 6*n + 5*n*lgn
743  ELSE
744  lwedc = 8
745  liwedc = 12
746  END IF
747  nap = ( n*( n+1 ) ) / 2
748  aninv = one / dble( max( 1, n ) )
749 *
750  IF( nsizes.NE.1 ) THEN
751  mtypes = min( maxtyp, ntypes )
752  ELSE
753  mtypes = min( maxtyp+1, ntypes )
754  END IF
755 *
756  DO 300 jtype = 1, mtypes
757  IF( .NOT.dotype( jtype ) )
758  $ GO TO 300
759  nmats = nmats + 1
760  ntest = 0
761 *
762  DO 30 j = 1, 4
763  ioldsd( j ) = iseed( j )
764  30 CONTINUE
765 *
766 * Compute "A"
767 *
768 * Control parameters:
769 *
770 * KMAGN KMODE KTYPE
771 * =1 O(1) clustered 1 zero
772 * =2 large clustered 2 identity
773 * =3 small exponential (none)
774 * =4 arithmetic diagonal, (w/ eigenvalues)
775 * =5 random log symmetric, w/ eigenvalues
776 * =6 random (none)
777 * =7 random diagonal
778 * =8 random symmetric
779 * =9 positive definite
780 * =10 diagonally dominant tridiagonal
781 *
782  IF( mtypes.GT.maxtyp )
783  $ GO TO 100
784 *
785  itype = ktype( jtype )
786  imode = kmode( jtype )
787 *
788 * Compute norm
789 *
790  GO TO ( 40, 50, 60 )kmagn( jtype )
791 *
792  40 CONTINUE
793  anorm = one
794  GO TO 70
795 *
796  50 CONTINUE
797  anorm = ( rtovfl*ulp )*aninv
798  GO TO 70
799 *
800  60 CONTINUE
801  anorm = rtunfl*n*ulpinv
802  GO TO 70
803 *
804  70 CONTINUE
805 *
806  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
807  iinfo = 0
808  IF( jtype.LE.15 ) THEN
809  cond = ulpinv
810  ELSE
811  cond = ulpinv*aninv / ten
812  END IF
813 *
814 * Special Matrices -- Identity & Jordan block
815 *
816 * Zero
817 *
818  IF( itype.EQ.1 ) THEN
819  iinfo = 0
820 *
821  ELSE IF( itype.EQ.2 ) THEN
822 *
823 * Identity
824 *
825  DO 80 jc = 1, n
826  a( jc, jc ) = anorm
827  80 CONTINUE
828 *
829  ELSE IF( itype.EQ.4 ) THEN
830 *
831 * Diagonal Matrix, [Eigen]values Specified
832 *
833  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
834  $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
835  $ iinfo )
836 *
837 *
838  ELSE IF( itype.EQ.5 ) THEN
839 *
840 * Symmetric, eigenvalues specified
841 *
842  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
843  $ anorm, n, n, 'N', a, lda, work( n+1 ),
844  $ iinfo )
845 *
846  ELSE IF( itype.EQ.7 ) THEN
847 *
848 * Diagonal, random eigenvalues
849 *
850  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
851  $ 'T', 'N', work( n+1 ), 1, one,
852  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
853  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
854 *
855  ELSE IF( itype.EQ.8 ) THEN
856 *
857 * Symmetric, random eigenvalues
858 *
859  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
860  $ 'T', 'N', work( n+1 ), 1, one,
861  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
862  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
863 *
864  ELSE IF( itype.EQ.9 ) THEN
865 *
866 * Positive definite, eigenvalues specified.
867 *
868  CALL dlatms( n, n, 'S', iseed, 'P', work, imode, cond,
869  $ anorm, n, n, 'N', a, lda, work( n+1 ),
870  $ iinfo )
871 *
872  ELSE IF( itype.EQ.10 ) THEN
873 *
874 * Positive definite tridiagonal, eigenvalues specified.
875 *
876  CALL dlatms( n, n, 'S', iseed, 'P', work, imode, cond,
877  $ anorm, 1, 1, 'N', a, lda, work( n+1 ),
878  $ iinfo )
879  DO 90 i = 2, n
880  temp1 = abs( a( i-1, i ) ) /
881  $ sqrt( abs( a( i-1, i-1 )*a( i, i ) ) )
882  IF( temp1.GT.half ) THEN
883  a( i-1, i ) = half*sqrt( abs( a( i-1, i-1 )*a( i,
884  $ i ) ) )
885  a( i, i-1 ) = a( i-1, i )
886  END IF
887  90 CONTINUE
888 *
889  ELSE
890 *
891  iinfo = 1
892  END IF
893 *
894  IF( iinfo.NE.0 ) THEN
895  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
896  $ ioldsd
897  info = abs( iinfo )
898  RETURN
899  END IF
900 *
901  100 CONTINUE
902 *
903 * Call DSYTRD and DORGTR to compute S and U from
904 * upper triangle.
905 *
906  CALL dlacpy( 'U', n, n, a, lda, v, ldu )
907 *
908  ntest = 1
909  CALL dsytrd( 'U', n, v, ldu, sd, se, tau, work, lwork,
910  $ iinfo )
911 *
912  IF( iinfo.NE.0 ) THEN
913  WRITE( nounit, fmt = 9999 )'DSYTRD(U)', iinfo, n, jtype,
914  $ ioldsd
915  info = abs( iinfo )
916  IF( iinfo.LT.0 ) THEN
917  RETURN
918  ELSE
919  result( 1 ) = ulpinv
920  GO TO 280
921  END IF
922  END IF
923 *
924  CALL dlacpy( 'U', n, n, v, ldu, u, ldu )
925 *
926  ntest = 2
927  CALL dorgtr( 'U', n, u, ldu, tau, work, lwork, iinfo )
928  IF( iinfo.NE.0 ) THEN
929  WRITE( nounit, fmt = 9999 )'DORGTR(U)', iinfo, n, jtype,
930  $ ioldsd
931  info = abs( iinfo )
932  IF( iinfo.LT.0 ) THEN
933  RETURN
934  ELSE
935  result( 2 ) = ulpinv
936  GO TO 280
937  END IF
938  END IF
939 *
940 * Do tests 1 and 2
941 *
942  CALL dsyt21( 2, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
943  $ ldu, tau, work, result( 1 ) )
944  CALL dsyt21( 3, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
945  $ ldu, tau, work, result( 2 ) )
946 *
947 * Call DSYTRD and DORGTR to compute S and U from
948 * lower triangle, do tests.
949 *
950  CALL dlacpy( 'L', n, n, a, lda, v, ldu )
951 *
952  ntest = 3
953  CALL dsytrd( 'L', n, v, ldu, sd, se, tau, work, lwork,
954  $ iinfo )
955 *
956  IF( iinfo.NE.0 ) THEN
957  WRITE( nounit, fmt = 9999 )'DSYTRD(L)', iinfo, n, jtype,
958  $ ioldsd
959  info = abs( iinfo )
960  IF( iinfo.LT.0 ) THEN
961  RETURN
962  ELSE
963  result( 3 ) = ulpinv
964  GO TO 280
965  END IF
966  END IF
967 *
968  CALL dlacpy( 'L', n, n, v, ldu, u, ldu )
969 *
970  ntest = 4
971  CALL dorgtr( 'L', n, u, ldu, tau, work, lwork, iinfo )
972  IF( iinfo.NE.0 ) THEN
973  WRITE( nounit, fmt = 9999 )'DORGTR(L)', iinfo, n, jtype,
974  $ ioldsd
975  info = abs( iinfo )
976  IF( iinfo.LT.0 ) THEN
977  RETURN
978  ELSE
979  result( 4 ) = ulpinv
980  GO TO 280
981  END IF
982  END IF
983 *
984  CALL dsyt21( 2, 'Lower', n, 1, a, lda, sd, se, u, ldu, v,
985  $ ldu, tau, work, result( 3 ) )
986  CALL dsyt21( 3, 'Lower', n, 1, a, lda, sd, se, u, ldu, v,
987  $ ldu, tau, work, result( 4 ) )
988 *
989 * Store the upper triangle of A in AP
990 *
991  i = 0
992  DO 120 jc = 1, n
993  DO 110 jr = 1, jc
994  i = i + 1
995  ap( i ) = a( jr, jc )
996  110 CONTINUE
997  120 CONTINUE
998 *
999 * Call DSPTRD and DOPGTR to compute S and U from AP
1000 *
1001  CALL dcopy( nap, ap, 1, vp, 1 )
1002 *
1003  ntest = 5
1004  CALL dsptrd( 'U', n, vp, sd, se, tau, iinfo )
1005 *
1006  IF( iinfo.NE.0 ) THEN
1007  WRITE( nounit, fmt = 9999 )'DSPTRD(U)', iinfo, n, jtype,
1008  $ ioldsd
1009  info = abs( iinfo )
1010  IF( iinfo.LT.0 ) THEN
1011  RETURN
1012  ELSE
1013  result( 5 ) = ulpinv
1014  GO TO 280
1015  END IF
1016  END IF
1017 *
1018  ntest = 6
1019  CALL dopgtr( 'U', n, vp, tau, u, ldu, work, iinfo )
1020  IF( iinfo.NE.0 ) THEN
1021  WRITE( nounit, fmt = 9999 )'DOPGTR(U)', iinfo, n, jtype,
1022  $ ioldsd
1023  info = abs( iinfo )
1024  IF( iinfo.LT.0 ) THEN
1025  RETURN
1026  ELSE
1027  result( 6 ) = ulpinv
1028  GO TO 280
1029  END IF
1030  END IF
1031 *
1032 * Do tests 5 and 6
1033 *
1034  CALL dspt21( 2, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1035  $ work, result( 5 ) )
1036  CALL dspt21( 3, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1037  $ work, result( 6 ) )
1038 *
1039 * Store the lower triangle of A in AP
1040 *
1041  i = 0
1042  DO 140 jc = 1, n
1043  DO 130 jr = jc, n
1044  i = i + 1
1045  ap( i ) = a( jr, jc )
1046  130 CONTINUE
1047  140 CONTINUE
1048 *
1049 * Call DSPTRD and DOPGTR to compute S and U from AP
1050 *
1051  CALL dcopy( nap, ap, 1, vp, 1 )
1052 *
1053  ntest = 7
1054  CALL dsptrd( 'L', n, vp, sd, se, tau, iinfo )
1055 *
1056  IF( iinfo.NE.0 ) THEN
1057  WRITE( nounit, fmt = 9999 )'DSPTRD(L)', iinfo, n, jtype,
1058  $ ioldsd
1059  info = abs( iinfo )
1060  IF( iinfo.LT.0 ) THEN
1061  RETURN
1062  ELSE
1063  result( 7 ) = ulpinv
1064  GO TO 280
1065  END IF
1066  END IF
1067 *
1068  ntest = 8
1069  CALL dopgtr( 'L', n, vp, tau, u, ldu, work, iinfo )
1070  IF( iinfo.NE.0 ) THEN
1071  WRITE( nounit, fmt = 9999 )'DOPGTR(L)', iinfo, n, jtype,
1072  $ ioldsd
1073  info = abs( iinfo )
1074  IF( iinfo.LT.0 ) THEN
1075  RETURN
1076  ELSE
1077  result( 8 ) = ulpinv
1078  GO TO 280
1079  END IF
1080  END IF
1081 *
1082  CALL dspt21( 2, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1083  $ work, result( 7 ) )
1084  CALL dspt21( 3, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1085  $ work, result( 8 ) )
1086 *
1087 * Call DSTEQR to compute D1, D2, and Z, do tests.
1088 *
1089 * Compute D1 and Z
1090 *
1091  CALL dcopy( n, sd, 1, d1, 1 )
1092  IF( n.GT.0 )
1093  $ CALL dcopy( n-1, se, 1, work, 1 )
1094  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1095 *
1096  ntest = 9
1097  CALL dsteqr( 'V', n, d1, work, z, ldu, work( n+1 ), iinfo )
1098  IF( iinfo.NE.0 ) THEN
1099  WRITE( nounit, fmt = 9999 )'DSTEQR(V)', iinfo, n, jtype,
1100  $ ioldsd
1101  info = abs( iinfo )
1102  IF( iinfo.LT.0 ) THEN
1103  RETURN
1104  ELSE
1105  result( 9 ) = ulpinv
1106  GO TO 280
1107  END IF
1108  END IF
1109 *
1110 * Compute D2
1111 *
1112  CALL dcopy( n, sd, 1, d2, 1 )
1113  IF( n.GT.0 )
1114  $ CALL dcopy( n-1, se, 1, work, 1 )
1115 *
1116  ntest = 11
1117  CALL dsteqr( 'N', n, d2, work, work( n+1 ), ldu,
1118  $ work( n+1 ), iinfo )
1119  IF( iinfo.NE.0 ) THEN
1120  WRITE( nounit, fmt = 9999 )'DSTEQR(N)', iinfo, n, jtype,
1121  $ ioldsd
1122  info = abs( iinfo )
1123  IF( iinfo.LT.0 ) THEN
1124  RETURN
1125  ELSE
1126  result( 11 ) = ulpinv
1127  GO TO 280
1128  END IF
1129  END IF
1130 *
1131 * Compute D3 (using PWK method)
1132 *
1133  CALL dcopy( n, sd, 1, d3, 1 )
1134  IF( n.GT.0 )
1135  $ CALL dcopy( n-1, se, 1, work, 1 )
1136 *
1137  ntest = 12
1138  CALL dsterf( n, d3, work, iinfo )
1139  IF( iinfo.NE.0 ) THEN
1140  WRITE( nounit, fmt = 9999 )'DSTERF', iinfo, n, jtype,
1141  $ ioldsd
1142  info = abs( iinfo )
1143  IF( iinfo.LT.0 ) THEN
1144  RETURN
1145  ELSE
1146  result( 12 ) = ulpinv
1147  GO TO 280
1148  END IF
1149  END IF
1150 *
1151 * Do Tests 9 and 10
1152 *
1153  CALL dstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1154  $ result( 9 ) )
1155 *
1156 * Do Tests 11 and 12
1157 *
1158  temp1 = zero
1159  temp2 = zero
1160  temp3 = zero
1161  temp4 = zero
1162 *
1163  DO 150 j = 1, n
1164  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1165  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1166  temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
1167  temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
1168  150 CONTINUE
1169 *
1170  result( 11 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1171  result( 12 ) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
1172 *
1173 * Do Test 13 -- Sturm Sequence Test of Eigenvalues
1174 * Go up by factors of two until it succeeds
1175 *
1176  ntest = 13
1177  temp1 = thresh*( half-ulp )
1178 *
1179  DO 160 j = 0, log2ui
1180  CALL dstech( n, sd, se, d1, temp1, work, iinfo )
1181  IF( iinfo.EQ.0 )
1182  $ GO TO 170
1183  temp1 = temp1*two
1184  160 CONTINUE
1185 *
1186  170 CONTINUE
1187  result( 13 ) = temp1
1188 *
1189 * For positive definite matrices ( JTYPE.GT.15 ) call DPTEQR
1190 * and do tests 14, 15, and 16 .
1191 *
1192  IF( jtype.GT.15 ) THEN
1193 *
1194 * Compute D4 and Z4
1195 *
1196  CALL dcopy( n, sd, 1, d4, 1 )
1197  IF( n.GT.0 )
1198  $ CALL dcopy( n-1, se, 1, work, 1 )
1199  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1200 *
1201  ntest = 14
1202  CALL dpteqr( 'V', n, d4, work, z, ldu, work( n+1 ),
1203  $ iinfo )
1204  IF( iinfo.NE.0 ) THEN
1205  WRITE( nounit, fmt = 9999 )'DPTEQR(V)', iinfo, n,
1206  $ jtype, ioldsd
1207  info = abs( iinfo )
1208  IF( iinfo.LT.0 ) THEN
1209  RETURN
1210  ELSE
1211  result( 14 ) = ulpinv
1212  GO TO 280
1213  END IF
1214  END IF
1215 *
1216 * Do Tests 14 and 15
1217 *
1218  CALL dstt21( n, 0, sd, se, d4, dumma, z, ldu, work,
1219  $ result( 14 ) )
1220 *
1221 * Compute D5
1222 *
1223  CALL dcopy( n, sd, 1, d5, 1 )
1224  IF( n.GT.0 )
1225  $ CALL dcopy( n-1, se, 1, work, 1 )
1226 *
1227  ntest = 16
1228  CALL dpteqr( 'N', n, d5, work, z, ldu, work( n+1 ),
1229  $ iinfo )
1230  IF( iinfo.NE.0 ) THEN
1231  WRITE( nounit, fmt = 9999 )'DPTEQR(N)', iinfo, n,
1232  $ jtype, ioldsd
1233  info = abs( iinfo )
1234  IF( iinfo.LT.0 ) THEN
1235  RETURN
1236  ELSE
1237  result( 16 ) = ulpinv
1238  GO TO 280
1239  END IF
1240  END IF
1241 *
1242 * Do Test 16
1243 *
1244  temp1 = zero
1245  temp2 = zero
1246  DO 180 j = 1, n
1247  temp1 = max( temp1, abs( d4( j ) ), abs( d5( j ) ) )
1248  temp2 = max( temp2, abs( d4( j )-d5( j ) ) )
1249  180 CONTINUE
1250 *
1251  result( 16 ) = temp2 / max( unfl,
1252  $ hun*ulp*max( temp1, temp2 ) )
1253  ELSE
1254  result( 14 ) = zero
1255  result( 15 ) = zero
1256  result( 16 ) = zero
1257  END IF
1258 *
1259 * Call DSTEBZ with different options and do tests 17-18.
1260 *
1261 * If S is positive definite and diagonally dominant,
1262 * ask for all eigenvalues with high relative accuracy.
1263 *
1264  vl = zero
1265  vu = zero
1266  il = 0
1267  iu = 0
1268  IF( jtype.EQ.21 ) THEN
1269  ntest = 17
1270  abstol = unfl + unfl
1271  CALL dstebz( 'A', 'E', n, vl, vu, il, iu, abstol, sd, se,
1272  $ m, nsplit, wr, iwork( 1 ), iwork( n+1 ),
1273  $ work, iwork( 2*n+1 ), iinfo )
1274  IF( iinfo.NE.0 ) THEN
1275  WRITE( nounit, fmt = 9999 )'DSTEBZ(A,rel)', iinfo, n,
1276  $ jtype, ioldsd
1277  info = abs( iinfo )
1278  IF( iinfo.LT.0 ) THEN
1279  RETURN
1280  ELSE
1281  result( 17 ) = ulpinv
1282  GO TO 280
1283  END IF
1284  END IF
1285 *
1286 * Do test 17
1287 *
1288  temp2 = two*( two*n-one )*ulp*( one+eight*half**2 ) /
1289  $ ( one-half )**4
1290 *
1291  temp1 = zero
1292  DO 190 j = 1, n
1293  temp1 = max( temp1, abs( d4( j )-wr( n-j+1 ) ) /
1294  $ ( abstol+abs( d4( j ) ) ) )
1295  190 CONTINUE
1296 *
1297  result( 17 ) = temp1 / temp2
1298  ELSE
1299  result( 17 ) = zero
1300  END IF
1301 *
1302 * Now ask for all eigenvalues with high absolute accuracy.
1303 *
1304  ntest = 18
1305  abstol = unfl + unfl
1306  CALL dstebz( 'A', 'E', n, vl, vu, il, iu, abstol, sd, se, m,
1307  $ nsplit, wa1, iwork( 1 ), iwork( n+1 ), work,
1308  $ iwork( 2*n+1 ), iinfo )
1309  IF( iinfo.NE.0 ) THEN
1310  WRITE( nounit, fmt = 9999 )'DSTEBZ(A)', iinfo, n, jtype,
1311  $ ioldsd
1312  info = abs( iinfo )
1313  IF( iinfo.LT.0 ) THEN
1314  RETURN
1315  ELSE
1316  result( 18 ) = ulpinv
1317  GO TO 280
1318  END IF
1319  END IF
1320 *
1321 * Do test 18
1322 *
1323  temp1 = zero
1324  temp2 = zero
1325  DO 200 j = 1, n
1326  temp1 = max( temp1, abs( d3( j ) ), abs( wa1( j ) ) )
1327  temp2 = max( temp2, abs( d3( j )-wa1( j ) ) )
1328  200 CONTINUE
1329 *
1330  result( 18 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1331 *
1332 * Choose random values for IL and IU, and ask for the
1333 * IL-th through IU-th eigenvalues.
1334 *
1335  ntest = 19
1336  IF( n.LE.1 ) THEN
1337  il = 1
1338  iu = n
1339  ELSE
1340  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1341  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1342  IF( iu.LT.il ) THEN
1343  itemp = iu
1344  iu = il
1345  il = itemp
1346  END IF
1347  END IF
1348 *
1349  CALL dstebz( 'I', 'E', n, vl, vu, il, iu, abstol, sd, se,
1350  $ m2, nsplit, wa2, iwork( 1 ), iwork( n+1 ),
1351  $ work, iwork( 2*n+1 ), iinfo )
1352  IF( iinfo.NE.0 ) THEN
1353  WRITE( nounit, fmt = 9999 )'DSTEBZ(I)', iinfo, n, jtype,
1354  $ ioldsd
1355  info = abs( iinfo )
1356  IF( iinfo.LT.0 ) THEN
1357  RETURN
1358  ELSE
1359  result( 19 ) = ulpinv
1360  GO TO 280
1361  END IF
1362  END IF
1363 *
1364 * Determine the values VL and VU of the IL-th and IU-th
1365 * eigenvalues and ask for all eigenvalues in this range.
1366 *
1367  IF( n.GT.0 ) THEN
1368  IF( il.NE.1 ) THEN
1369  vl = wa1( il ) - max( half*( wa1( il )-wa1( il-1 ) ),
1370  $ ulp*anorm, two*rtunfl )
1371  ELSE
1372  vl = wa1( 1 ) - max( half*( wa1( n )-wa1( 1 ) ),
1373  $ ulp*anorm, two*rtunfl )
1374  END IF
1375  IF( iu.NE.n ) THEN
1376  vu = wa1( iu ) + max( half*( wa1( iu+1 )-wa1( iu ) ),
1377  $ ulp*anorm, two*rtunfl )
1378  ELSE
1379  vu = wa1( n ) + max( half*( wa1( n )-wa1( 1 ) ),
1380  $ ulp*anorm, two*rtunfl )
1381  END IF
1382  ELSE
1383  vl = zero
1384  vu = one
1385  END IF
1386 *
1387  CALL dstebz( 'V', 'E', n, vl, vu, il, iu, abstol, sd, se,
1388  $ m3, nsplit, wa3, iwork( 1 ), iwork( n+1 ),
1389  $ work, iwork( 2*n+1 ), iinfo )
1390  IF( iinfo.NE.0 ) THEN
1391  WRITE( nounit, fmt = 9999 )'DSTEBZ(V)', iinfo, n, jtype,
1392  $ ioldsd
1393  info = abs( iinfo )
1394  IF( iinfo.LT.0 ) THEN
1395  RETURN
1396  ELSE
1397  result( 19 ) = ulpinv
1398  GO TO 280
1399  END IF
1400  END IF
1401 *
1402  IF( m3.EQ.0 .AND. n.NE.0 ) THEN
1403  result( 19 ) = ulpinv
1404  GO TO 280
1405  END IF
1406 *
1407 * Do test 19
1408 *
1409  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1410  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1411  IF( n.GT.0 ) THEN
1412  temp3 = max( abs( wa1( n ) ), abs( wa1( 1 ) ) )
1413  ELSE
1414  temp3 = zero
1415  END IF
1416 *
1417  result( 19 ) = ( temp1+temp2 ) / max( unfl, temp3*ulp )
1418 *
1419 * Call DSTEIN to compute eigenvectors corresponding to
1420 * eigenvalues in WA1. (First call DSTEBZ again, to make sure
1421 * it returns these eigenvalues in the correct order.)
1422 *
1423  ntest = 21
1424  CALL dstebz( 'A', 'B', n, vl, vu, il, iu, abstol, sd, se, m,
1425  $ nsplit, wa1, iwork( 1 ), iwork( n+1 ), work,
1426  $ iwork( 2*n+1 ), iinfo )
1427  IF( iinfo.NE.0 ) THEN
1428  WRITE( nounit, fmt = 9999 )'DSTEBZ(A,B)', iinfo, n,
1429  $ jtype, ioldsd
1430  info = abs( iinfo )
1431  IF( iinfo.LT.0 ) THEN
1432  RETURN
1433  ELSE
1434  result( 20 ) = ulpinv
1435  result( 21 ) = ulpinv
1436  GO TO 280
1437  END IF
1438  END IF
1439 *
1440  CALL dstein( n, sd, se, m, wa1, iwork( 1 ), iwork( n+1 ), z,
1441  $ ldu, work, iwork( 2*n+1 ), iwork( 3*n+1 ),
1442  $ iinfo )
1443  IF( iinfo.NE.0 ) THEN
1444  WRITE( nounit, fmt = 9999 )'DSTEIN', iinfo, n, jtype,
1445  $ ioldsd
1446  info = abs( iinfo )
1447  IF( iinfo.LT.0 ) THEN
1448  RETURN
1449  ELSE
1450  result( 20 ) = ulpinv
1451  result( 21 ) = ulpinv
1452  GO TO 280
1453  END IF
1454  END IF
1455 *
1456 * Do tests 20 and 21
1457 *
1458  CALL dstt21( n, 0, sd, se, wa1, dumma, z, ldu, work,
1459  $ result( 20 ) )
1460 *
1461 * Call DSTEDC(I) to compute D1 and Z, do tests.
1462 *
1463 * Compute D1 and Z
1464 *
1465  CALL dcopy( n, sd, 1, d1, 1 )
1466  IF( n.GT.0 )
1467  $ CALL dcopy( n-1, se, 1, work, 1 )
1468  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1469 *
1470  ntest = 22
1471  CALL dstedc( 'I', n, d1, work, z, ldu, work( n+1 ), lwedc-n,
1472  $ iwork, liwedc, iinfo )
1473  IF( iinfo.NE.0 ) THEN
1474  WRITE( nounit, fmt = 9999 )'DSTEDC(I)', iinfo, n, jtype,
1475  $ ioldsd
1476  info = abs( iinfo )
1477  IF( iinfo.LT.0 ) THEN
1478  RETURN
1479  ELSE
1480  result( 22 ) = ulpinv
1481  GO TO 280
1482  END IF
1483  END IF
1484 *
1485 * Do Tests 22 and 23
1486 *
1487  CALL dstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1488  $ result( 22 ) )
1489 *
1490 * Call DSTEDC(V) to compute D1 and Z, do tests.
1491 *
1492 * Compute D1 and Z
1493 *
1494  CALL dcopy( n, sd, 1, d1, 1 )
1495  IF( n.GT.0 )
1496  $ CALL dcopy( n-1, se, 1, work, 1 )
1497  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1498 *
1499  ntest = 24
1500  CALL dstedc( 'V', n, d1, work, z, ldu, work( n+1 ), lwedc-n,
1501  $ iwork, liwedc, iinfo )
1502  IF( iinfo.NE.0 ) THEN
1503  WRITE( nounit, fmt = 9999 )'DSTEDC(V)', iinfo, n, jtype,
1504  $ ioldsd
1505  info = abs( iinfo )
1506  IF( iinfo.LT.0 ) THEN
1507  RETURN
1508  ELSE
1509  result( 24 ) = ulpinv
1510  GO TO 280
1511  END IF
1512  END IF
1513 *
1514 * Do Tests 24 and 25
1515 *
1516  CALL dstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1517  $ result( 24 ) )
1518 *
1519 * Call DSTEDC(N) to compute D2, do tests.
1520 *
1521 * Compute D2
1522 *
1523  CALL dcopy( n, sd, 1, d2, 1 )
1524  IF( n.GT.0 )
1525  $ CALL dcopy( n-1, se, 1, work, 1 )
1526  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1527 *
1528  ntest = 26
1529  CALL dstedc( 'N', n, d2, work, z, ldu, work( n+1 ), lwedc-n,
1530  $ iwork, liwedc, iinfo )
1531  IF( iinfo.NE.0 ) THEN
1532  WRITE( nounit, fmt = 9999 )'DSTEDC(N)', iinfo, n, jtype,
1533  $ ioldsd
1534  info = abs( iinfo )
1535  IF( iinfo.LT.0 ) THEN
1536  RETURN
1537  ELSE
1538  result( 26 ) = ulpinv
1539  GO TO 280
1540  END IF
1541  END IF
1542 *
1543 * Do Test 26
1544 *
1545  temp1 = zero
1546  temp2 = zero
1547 *
1548  DO 210 j = 1, n
1549  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1550  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1551  210 CONTINUE
1552 *
1553  result( 26 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1554 *
1555 * Only test DSTEMR if IEEE compliant
1556 *
1557  IF( ilaenv( 10, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
1558  $ ilaenv( 11, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
1559 *
1560 * Call DSTEMR, do test 27 (relative eigenvalue accuracy)
1561 *
1562 * If S is positive definite and diagonally dominant,
1563 * ask for all eigenvalues with high relative accuracy.
1564 *
1565  vl = zero
1566  vu = zero
1567  il = 0
1568  iu = 0
1569  IF( jtype.EQ.21 .AND. srel ) THEN
1570  ntest = 27
1571  abstol = unfl + unfl
1572  CALL dstemr( 'V', 'A', n, sd, se, vl, vu, il, iu,
1573  $ m, wr, z, ldu, n, iwork( 1 ), tryrac,
1574  $ work, lwork, iwork( 2*n+1 ), lwork-2*n,
1575  $ iinfo )
1576  IF( iinfo.NE.0 ) THEN
1577  WRITE( nounit, fmt = 9999 )'DSTEMR(V,A,rel)',
1578  $ iinfo, n, jtype, ioldsd
1579  info = abs( iinfo )
1580  IF( iinfo.LT.0 ) THEN
1581  RETURN
1582  ELSE
1583  result( 27 ) = ulpinv
1584  GO TO 270
1585  END IF
1586  END IF
1587 *
1588 * Do test 27
1589 *
1590  temp2 = two*( two*n-one )*ulp*( one+eight*half**2 ) /
1591  $ ( one-half )**4
1592 *
1593  temp1 = zero
1594  DO 220 j = 1, n
1595  temp1 = max( temp1, abs( d4( j )-wr( n-j+1 ) ) /
1596  $ ( abstol+abs( d4( j ) ) ) )
1597  220 CONTINUE
1598 *
1599  result( 27 ) = temp1 / temp2
1600 *
1601  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1602  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1603  IF( iu.LT.il ) THEN
1604  itemp = iu
1605  iu = il
1606  il = itemp
1607  END IF
1608 *
1609  IF( srange ) THEN
1610  ntest = 28
1611  abstol = unfl + unfl
1612  CALL dstemr( 'V', 'I', n, sd, se, vl, vu, il, iu,
1613  $ m, wr, z, ldu, n, iwork( 1 ), tryrac,
1614  $ work, lwork, iwork( 2*n+1 ),
1615  $ lwork-2*n, iinfo )
1616 *
1617  IF( iinfo.NE.0 ) THEN
1618  WRITE( nounit, fmt = 9999 )'DSTEMR(V,I,rel)',
1619  $ iinfo, n, jtype, ioldsd
1620  info = abs( iinfo )
1621  IF( iinfo.LT.0 ) THEN
1622  RETURN
1623  ELSE
1624  result( 28 ) = ulpinv
1625  GO TO 270
1626  END IF
1627  END IF
1628 *
1629 *
1630 * Do test 28
1631 *
1632  temp2 = two*( two*n-one )*ulp*
1633  $ ( one+eight*half**2 ) / ( one-half )**4
1634 *
1635  temp1 = zero
1636  DO 230 j = il, iu
1637  temp1 = max( temp1, abs( wr( j-il+1 )-d4( n-j+
1638  $ 1 ) ) / ( abstol+abs( wr( j-il+1 ) ) ) )
1639  230 CONTINUE
1640 *
1641  result( 28 ) = temp1 / temp2
1642  ELSE
1643  result( 28 ) = zero
1644  END IF
1645  ELSE
1646  result( 27 ) = zero
1647  result( 28 ) = zero
1648  END IF
1649 *
1650 * Call DSTEMR(V,I) to compute D1 and Z, do tests.
1651 *
1652 * Compute D1 and Z
1653 *
1654  CALL dcopy( n, sd, 1, d5, 1 )
1655  IF( n.GT.0 )
1656  $ CALL dcopy( n-1, se, 1, work, 1 )
1657  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1658 *
1659  IF( srange ) THEN
1660  ntest = 29
1661  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1662  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1663  IF( iu.LT.il ) THEN
1664  itemp = iu
1665  iu = il
1666  il = itemp
1667  END IF
1668  CALL dstemr( 'V', 'I', n, d5, work, vl, vu, il, iu,
1669  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1670  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1671  $ liwork-2*n, iinfo )
1672  IF( iinfo.NE.0 ) THEN
1673  WRITE( nounit, fmt = 9999 )'DSTEMR(V,I)', iinfo,
1674  $ n, jtype, ioldsd
1675  info = abs( iinfo )
1676  IF( iinfo.LT.0 ) THEN
1677  RETURN
1678  ELSE
1679  result( 29 ) = ulpinv
1680  GO TO 280
1681  END IF
1682  END IF
1683 *
1684 * Do Tests 29 and 30
1685 *
1686  CALL dstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1687  $ m, result( 29 ) )
1688 *
1689 * Call DSTEMR to compute D2, do tests.
1690 *
1691 * Compute D2
1692 *
1693  CALL dcopy( n, sd, 1, d5, 1 )
1694  IF( n.GT.0 )
1695  $ CALL dcopy( n-1, se, 1, work, 1 )
1696 *
1697  ntest = 31
1698  CALL dstemr( 'N', 'I', n, d5, work, vl, vu, il, iu,
1699  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1700  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1701  $ liwork-2*n, iinfo )
1702  IF( iinfo.NE.0 ) THEN
1703  WRITE( nounit, fmt = 9999 )'DSTEMR(N,I)', iinfo,
1704  $ n, jtype, ioldsd
1705  info = abs( iinfo )
1706  IF( iinfo.LT.0 ) THEN
1707  RETURN
1708  ELSE
1709  result( 31 ) = ulpinv
1710  GO TO 280
1711  END IF
1712  END IF
1713 *
1714 * Do Test 31
1715 *
1716  temp1 = zero
1717  temp2 = zero
1718 *
1719  DO 240 j = 1, iu - il + 1
1720  temp1 = max( temp1, abs( d1( j ) ),
1721  $ abs( d2( j ) ) )
1722  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1723  240 CONTINUE
1724 *
1725  result( 31 ) = temp2 / max( unfl,
1726  $ ulp*max( temp1, temp2 ) )
1727 *
1728 *
1729 * Call DSTEMR(V,V) to compute D1 and Z, do tests.
1730 *
1731 * Compute D1 and Z
1732 *
1733  CALL dcopy( n, sd, 1, d5, 1 )
1734  IF( n.GT.0 )
1735  $ CALL dcopy( n-1, se, 1, work, 1 )
1736  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1737 *
1738  ntest = 32
1739 *
1740  IF( n.GT.0 ) THEN
1741  IF( il.NE.1 ) THEN
1742  vl = d2( il ) - max( half*
1743  $ ( d2( il )-d2( il-1 ) ), ulp*anorm,
1744  $ two*rtunfl )
1745  ELSE
1746  vl = d2( 1 ) - max( half*( d2( n )-d2( 1 ) ),
1747  $ ulp*anorm, two*rtunfl )
1748  END IF
1749  IF( iu.NE.n ) THEN
1750  vu = d2( iu ) + max( half*
1751  $ ( d2( iu+1 )-d2( iu ) ), ulp*anorm,
1752  $ two*rtunfl )
1753  ELSE
1754  vu = d2( n ) + max( half*( d2( n )-d2( 1 ) ),
1755  $ ulp*anorm, two*rtunfl )
1756  END IF
1757  ELSE
1758  vl = zero
1759  vu = one
1760  END IF
1761 *
1762  CALL dstemr( 'V', 'V', n, d5, work, vl, vu, il, iu,
1763  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1764  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1765  $ liwork-2*n, iinfo )
1766  IF( iinfo.NE.0 ) THEN
1767  WRITE( nounit, fmt = 9999 )'DSTEMR(V,V)', iinfo,
1768  $ n, jtype, ioldsd
1769  info = abs( iinfo )
1770  IF( iinfo.LT.0 ) THEN
1771  RETURN
1772  ELSE
1773  result( 32 ) = ulpinv
1774  GO TO 280
1775  END IF
1776  END IF
1777 *
1778 * Do Tests 32 and 33
1779 *
1780  CALL dstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1781  $ m, result( 32 ) )
1782 *
1783 * Call DSTEMR to compute D2, do tests.
1784 *
1785 * Compute D2
1786 *
1787  CALL dcopy( n, sd, 1, d5, 1 )
1788  IF( n.GT.0 )
1789  $ CALL dcopy( n-1, se, 1, work, 1 )
1790 *
1791  ntest = 34
1792  CALL dstemr( 'N', 'V', n, d5, work, vl, vu, il, iu,
1793  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1794  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1795  $ liwork-2*n, iinfo )
1796  IF( iinfo.NE.0 ) THEN
1797  WRITE( nounit, fmt = 9999 )'DSTEMR(N,V)', iinfo,
1798  $ n, jtype, ioldsd
1799  info = abs( iinfo )
1800  IF( iinfo.LT.0 ) THEN
1801  RETURN
1802  ELSE
1803  result( 34 ) = ulpinv
1804  GO TO 280
1805  END IF
1806  END IF
1807 *
1808 * Do Test 34
1809 *
1810  temp1 = zero
1811  temp2 = zero
1812 *
1813  DO 250 j = 1, iu - il + 1
1814  temp1 = max( temp1, abs( d1( j ) ),
1815  $ abs( d2( j ) ) )
1816  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1817  250 CONTINUE
1818 *
1819  result( 34 ) = temp2 / max( unfl,
1820  $ ulp*max( temp1, temp2 ) )
1821  ELSE
1822  result( 29 ) = zero
1823  result( 30 ) = zero
1824  result( 31 ) = zero
1825  result( 32 ) = zero
1826  result( 33 ) = zero
1827  result( 34 ) = zero
1828  END IF
1829 *
1830 *
1831 * Call DSTEMR(V,A) to compute D1 and Z, do tests.
1832 *
1833 * Compute D1 and Z
1834 *
1835  CALL dcopy( n, sd, 1, d5, 1 )
1836  IF( n.GT.0 )
1837  $ CALL dcopy( n-1, se, 1, work, 1 )
1838 *
1839  ntest = 35
1840 *
1841  CALL dstemr( 'V', 'A', n, d5, work, vl, vu, il, iu,
1842  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1843  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1844  $ liwork-2*n, iinfo )
1845  IF( iinfo.NE.0 ) THEN
1846  WRITE( nounit, fmt = 9999 )'DSTEMR(V,A)', iinfo, n,
1847  $ jtype, ioldsd
1848  info = abs( iinfo )
1849  IF( iinfo.LT.0 ) THEN
1850  RETURN
1851  ELSE
1852  result( 35 ) = ulpinv
1853  GO TO 280
1854  END IF
1855  END IF
1856 *
1857 * Do Tests 35 and 36
1858 *
1859  CALL dstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work, m,
1860  $ result( 35 ) )
1861 *
1862 * Call DSTEMR to compute D2, do tests.
1863 *
1864 * Compute D2
1865 *
1866  CALL dcopy( n, sd, 1, d5, 1 )
1867  IF( n.GT.0 )
1868  $ CALL dcopy( n-1, se, 1, work, 1 )
1869 *
1870  ntest = 37
1871  CALL dstemr( 'N', 'A', n, d5, work, vl, vu, il, iu,
1872  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1873  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1874  $ liwork-2*n, iinfo )
1875  IF( iinfo.NE.0 ) THEN
1876  WRITE( nounit, fmt = 9999 )'DSTEMR(N,A)', iinfo, n,
1877  $ jtype, ioldsd
1878  info = abs( iinfo )
1879  IF( iinfo.LT.0 ) THEN
1880  RETURN
1881  ELSE
1882  result( 37 ) = ulpinv
1883  GO TO 280
1884  END IF
1885  END IF
1886 *
1887 * Do Test 34
1888 *
1889  temp1 = zero
1890  temp2 = zero
1891 *
1892  DO 260 j = 1, n
1893  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1894  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1895  260 CONTINUE
1896 *
1897  result( 37 ) = temp2 / max( unfl,
1898  $ ulp*max( temp1, temp2 ) )
1899  END IF
1900  270 CONTINUE
1901  280 CONTINUE
1902  ntestt = ntestt + ntest
1903 *
1904 * End of Loop -- Check for RESULT(j) > THRESH
1905 *
1906 *
1907 * Print out tests which fail.
1908 *
1909  DO 290 jr = 1, ntest
1910  IF( result( jr ).GE.thresh ) THEN
1911 *
1912 * If this is the first test to fail,
1913 * print a header to the data file.
1914 *
1915  IF( nerrs.EQ.0 ) THEN
1916  WRITE( nounit, fmt = 9998 )'DST'
1917  WRITE( nounit, fmt = 9997 )
1918  WRITE( nounit, fmt = 9996 )
1919  WRITE( nounit, fmt = 9995 )'Symmetric'
1920  WRITE( nounit, fmt = 9994 )
1921 *
1922 * Tests performed
1923 *
1924  WRITE( nounit, fmt = 9988 )
1925  END IF
1926  nerrs = nerrs + 1
1927  WRITE( nounit, fmt = 9990 )n, ioldsd, jtype, jr,
1928  $ result( jr )
1929  END IF
1930  290 CONTINUE
1931  300 CONTINUE
1932  310 CONTINUE
1933 *
1934 * Summary
1935 *
1936  CALL dlasum( 'DST', nounit, nerrs, ntestt )
1937  RETURN
1938 *
1939  9999 FORMAT( ' DCHKST: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
1940  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1941 *
1942  9998 FORMAT( / 1x, a3, ' -- Real Symmetric eigenvalue problem' )
1943  9997 FORMAT( ' Matrix types (see DCHKST for details): ' )
1944 *
1945  9996 FORMAT( / ' Special Matrices:',
1946  $ / ' 1=Zero matrix. ',
1947  $ ' 5=Diagonal: clustered entries.',
1948  $ / ' 2=Identity matrix. ',
1949  $ ' 6=Diagonal: large, evenly spaced.',
1950  $ / ' 3=Diagonal: evenly spaced entries. ',
1951  $ ' 7=Diagonal: small, evenly spaced.',
1952  $ / ' 4=Diagonal: geometr. spaced entries.' )
1953  9995 FORMAT( ' Dense ', a, ' Matrices:',
1954  $ / ' 8=Evenly spaced eigenvals. ',
1955  $ ' 12=Small, evenly spaced eigenvals.',
1956  $ / ' 9=Geometrically spaced eigenvals. ',
1957  $ ' 13=Matrix with random O(1) entries.',
1958  $ / ' 10=Clustered eigenvalues. ',
1959  $ ' 14=Matrix with large random entries.',
1960  $ / ' 11=Large, evenly spaced eigenvals. ',
1961  $ ' 15=Matrix with small random entries.' )
1962  9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues',
1963  $ / ' 17=Positive definite, geometrically spaced eigenvlaues',
1964  $ / ' 18=Positive definite, clustered eigenvalues',
1965  $ / ' 19=Positive definite, small evenly spaced eigenvalues',
1966  $ / ' 20=Positive definite, large evenly spaced eigenvalues',
1967  $ / ' 21=Diagonally dominant tridiagonal, geometrically',
1968  $ ' spaced eigenvalues' )
1969 *
1970  9990 FORMAT( ' N=', i5, ', seed=', 4( i4, ',' ), ' type ', i2,
1971  $ ', test(', i2, ')=', g10.3 )
1972 *
1973  9988 FORMAT( / 'Test performed: see DCHKST for details.', / )
1974 * End of DCHKST
1975 *
1976  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:188
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:43
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:591
subroutine dstt22(N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, LDWORK, RESULT)
DSTT22
Definition: dstt22.f:139
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:101
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:127
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:471
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:321
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:321
subroutine dsptrd(UPLO, N, AP, D, E, TAU, INFO)
DSPTRD
Definition: dsptrd.f:150
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:174
subroutine dopgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
DOPGTR
Definition: dopgtr.f:114
subroutine dorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
DORGTR
Definition: dorgtr.f:123
subroutine dpteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DPTEQR
Definition: dpteqr.f:145
subroutine dsytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
DSYTRD
Definition: dsytrd.f:192