LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
schkst.f
Go to the documentation of this file.
1 *> \brief \b SCHKST
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 SCHKST( 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 * REAL THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
24 * REAL 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 *> SCHKST checks the symmetric eigenvalue problem routines.
38 *>
39 *> SSYTRD factors A as U S U' , where ' means transpose,
40 *> S is symmetric tridiagonal, and U is orthogonal.
41 *> SSYTRD can use either just the lower or just the upper triangle
42 *> of A; SCHKST 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 *> SSPTRD does the same as SSYTRD, except that A and V are stored
48 *> in "packed" format.
49 *>
50 *> SORGTR constructs the matrix U from the contents of V and TAU.
51 *>
52 *> SOPGTR constructs the matrix U from the contents of VP and TAU.
53 *>
54 *> SSTEQR 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 *> SSTERF computes D3, the matrix of eigenvalues, by the
60 *> PWK method, which does not yield eigenvectors.
61 *>
62 *> SPTEQR 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 *> SSTEBZ 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 *> SSTEIN computes Y, the eigenvectors of S, given the
74 *> eigenvalues.
75 *>
76 *> SSTEDC 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 SSYTRD/SORGTR or SSPTRD/SOPGTR ('V' option). It may
81 *> also just compute eigenvalues ('N' option).
82 *>
83 *> SSTEMR 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). SSTEMR
86 *> uses the Relatively Robust Representation whenever possible.
87 *>
88 *> When SCHKST 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 ) SSYTRD( UPLO='U', ... )
95 *>
96 *> (2) | I - UV' | / ( n ulp ) SORGTR( UPLO='U', ... )
97 *>
98 *> (3) | A - V S V' | / ( |A| n ulp ) SSYTRD( UPLO='L', ... )
99 *>
100 *> (4) | I - UV' | / ( n ulp ) SORGTR( UPLO='L', ... )
101 *>
102 *> (5-8) Same as 1-4, but for SSPTRD and SOPGTR.
103 *>
104 *> (9) | S - Z D Z' | / ( |S| n ulp ) SSTEQR('V',...)
105 *>
106 *> (10) | I - ZZ' | / ( n ulp ) SSTEQR('V',...)
107 *>
108 *> (11) | D1 - D2 | / ( |D1| ulp ) SSTEQR('N',...)
109 *>
110 *> (12) | D1 - D3 | / ( |D1| ulp ) SSTERF
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 *> SSTECH)
116 *>
117 *> For S positive definite,
118 *>
119 *> (14) | S - Z4 D4 Z4' | / ( |S| n ulp ) SPTEQR('V',...)
120 *>
121 *> (15) | I - Z4 Z4' | / ( n ulp ) SPTEQR('V',...)
122 *>
123 *> (16) | D4 - D5 | / ( 100 |D4| ulp ) SPTEQR('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 *> SSTEBZ( 'A', 'E', ...)
131 *>
132 *> (18) | WA1 - D3 | / ( |D3| ulp ) SSTEBZ( '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 *> SSTEBZ( 'I', 'E', ...)
139 *>
140 *> (20) | S - Y WA1 Y' | / ( |S| n ulp ) SSTEBZ, SSTEIN
141 *>
142 *> (21) | I - Y Y' | / ( n ulp ) SSTEBZ, SSTEIN
143 *>
144 *> (22) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('I')
145 *>
146 *> (23) | I - ZZ' | / ( n ulp ) SSTEDC('I')
147 *>
148 *> (24) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('V')
149 *>
150 *> (25) | I - ZZ' | / ( n ulp ) SSTEDC('V')
151 *>
152 *> (26) | D1 - D2 | / ( |D1| ulp ) SSTEDC('V') and
153 *> SSTEDC('N')
154 *>
155 *> Test 27 is disabled at the moment because SSTEMR 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 *> SSTEMR('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 *> SSTEMR('V', 'I')
167 *>
168 *> Tests 29 through 34 are disable at present because SSTEMR
169 *> does not handle partial specturm requests.
170 *>
171 *> (29) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'I')
172 *>
173 *> (30) | I - ZZ' | / ( n ulp ) SSTEMR('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 *> SSTEMR('N', 'I') vs. SSTEMR('V', 'I')
180 *>
181 *> (32) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'V')
182 *>
183 *> (33) | I - ZZ' | / ( n ulp ) SSTEMR('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 *> SSTEMR('N', 'V') vs. SSTEMR('V', 'V')
190 *>
191 *> (35) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'A')
192 *>
193 *> (36) | I - ZZ' | / ( n ulp ) SSTEMR('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 *> SSTEMR('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 *> SCHKST 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, SCHKST
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 SCHKST to continue the same random number
300 *> sequence.
301 *> \endverbatim
302 *>
303 *> \param[in] THRESH
304 *> \verbatim
305 *> THRESH is REAL
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 REAL 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 REAL 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 REAL array of
347 *> dimension( max(NN) )
348 *> The diagonal of the tridiagonal matrix computed by SSYTRD.
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 REAL array of
356 *> dimension( max(NN) )
357 *> The off-diagonal of the tridiagonal matrix computed by
358 *> SSYTRD. 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 REAL array of
365 *> dimension( max(NN) )
366 *> The eigenvalues of A, as computed by SSTEQR 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 REAL array of
374 *> dimension( max(NN) )
375 *> The eigenvalues of A, as computed by SSTEQR 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 REAL array of
383 *> dimension( max(NN) )
384 *> The eigenvalues of A, as computed by SSTERF. 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 REAL array of
391 *> dimension( max(NN) )
392 *> The eigenvalues of A, as computed by SPTEQR(V).
393 *> ZPTEQR 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 REAL array of
400 *> dimension( max(NN) )
401 *> The eigenvalues of A, as computed by SPTEQR(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 REAL 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 SSTEBZ.
413 *> \endverbatim
414 *>
415 *> \param[out] WA2
416 *> \verbatim
417 *> WA2 is REAL 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 SSTEBZ.
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 REAL 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 SSTEBZ.
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 REAL array of
440 *> dimension( max(NN) )
441 *> All eigenvalues of A, computed to high
442 *> absolute accuracy, with different options.
443 *> as computed by SSTEBZ.
444 *> \endverbatim
445 *>
446 *> \param[out] U
447 *> \verbatim
448 *> U is REAL array of
449 *> dimension( LDU, max(NN) ).
450 *> The orthogonal matrix computed by SSYTRD + SORGTR.
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 REAL array of
463 *> dimension( LDU, max(NN) ).
464 *> The Housholder vectors computed by SSYTRD 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 SSYTRD, 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 SORGTR, 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 REAL 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 REAL array of
484 *> dimension( max(NN) )
485 *> The Householder factors computed by SSYTRD in reducing A
486 *> to tridiagonal form.
487 *> \endverbatim
488 *>
489 *> \param[out] Z
490 *> \verbatim
491 *> Z is REAL array of
492 *> dimension( LDU, max(NN) ).
493 *> The orthogonal matrix of eigenvectors computed by SSTEQR,
494 *> SPTEQR, and SSTEIN.
495 *> \endverbatim
496 *>
497 *> \param[out] WORK
498 *> \verbatim
499 *> WORK is REAL 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 REAL 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 SLATMR, SLATMS, SSYTRD, SORGTR, SSTEQR, SSTERF,
545 *> or SORMC2 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 single_eig
587 *
588 * =====================================================================
589  SUBROUTINE schkst( 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  REAL THRESH
603 * ..
604 * .. Array Arguments ..
605  LOGICAL DOTYPE( * )
606  INTEGER ISEED( 4 ), IWORK( * ), NN( * )
607  REAL 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  REAL ZERO, ONE, TWO, EIGHT, TEN, HUN
618  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
619  $ eight = 8.0e0, ten = 10.0e0, hun = 100.0e0 )
620  REAL 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  REAL 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  REAL DUMMA( 1 )
644 * ..
645 * .. External Functions ..
646  INTEGER ILAENV
647  REAL SLAMCH, SLARND, SSXT1
648  EXTERNAL ilaenv, slamch, slarnd, ssxt1
649 * ..
650 * .. External Subroutines ..
651  EXTERNAL scopy, slabad, slacpy, slaset, slasum, slatmr,
655 * ..
656 * .. Intrinsic Functions ..
657  INTRINSIC abs, int, log, max, min, REAL, 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, 'SSYTRD', '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( 'SCHKST', -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 = slamch( 'Safe minimum' )
720  ovfl = one / unfl
721  CALL slabad( unfl, ovfl )
722  ulp = slamch( 'Epsilon' )*slamch( '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( REAL( 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 / REAL( 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 slaset( '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 slatms( 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 slatms( 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 slatmr( 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 slatmr( 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 slatms( 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 slatms( 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 SSYTRD and SORGTR to compute S and U from
907 * upper triangle.
908 *
909  CALL slacpy( 'U', n, n, a, lda, v, ldu )
910 *
911  ntest = 1
912  CALL ssytrd( 'U', n, v, ldu, sd, se, tau, work, lwork,
913  $ iinfo )
914 *
915  IF( iinfo.NE.0 ) THEN
916  WRITE( nounit, fmt = 9999 )'SSYTRD(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 slacpy( 'U', n, n, v, ldu, u, ldu )
928 *
929  ntest = 2
930  CALL sorgtr( 'U', n, u, ldu, tau, work, lwork, iinfo )
931  IF( iinfo.NE.0 ) THEN
932  WRITE( nounit, fmt = 9999 )'SORGTR(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 ssyt21( 2, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
946  $ ldu, tau, work, result( 1 ) )
947  CALL ssyt21( 3, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
948  $ ldu, tau, work, result( 2 ) )
949 *
950 * Call SSYTRD and SORGTR to compute S and U from
951 * lower triangle, do tests.
952 *
953  CALL slacpy( 'L', n, n, a, lda, v, ldu )
954 *
955  ntest = 3
956  CALL ssytrd( 'L', n, v, ldu, sd, se, tau, work, lwork,
957  $ iinfo )
958 *
959  IF( iinfo.NE.0 ) THEN
960  WRITE( nounit, fmt = 9999 )'SSYTRD(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 slacpy( 'L', n, n, v, ldu, u, ldu )
972 *
973  ntest = 4
974  CALL sorgtr( 'L', n, u, ldu, tau, work, lwork, iinfo )
975  IF( iinfo.NE.0 ) THEN
976  WRITE( nounit, fmt = 9999 )'SORGTR(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 ssyt21( 2, 'Lower', n, 1, a, lda, sd, se, u, ldu, v,
988  $ ldu, tau, work, result( 3 ) )
989  CALL ssyt21( 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 SSPTRD and SOPGTR to compute S and U from AP
1003 *
1004  CALL scopy( nap, ap, 1, vp, 1 )
1005 *
1006  ntest = 5
1007  CALL ssptrd( 'U', n, vp, sd, se, tau, iinfo )
1008 *
1009  IF( iinfo.NE.0 ) THEN
1010  WRITE( nounit, fmt = 9999 )'SSPTRD(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 sopgtr( 'U', n, vp, tau, u, ldu, work, iinfo )
1023  IF( iinfo.NE.0 ) THEN
1024  WRITE( nounit, fmt = 9999 )'SOPGTR(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 sspt21( 2, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1038  $ work, result( 5 ) )
1039  CALL sspt21( 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 SSPTRD and SOPGTR to compute S and U from AP
1053 *
1054  CALL scopy( nap, ap, 1, vp, 1 )
1055 *
1056  ntest = 7
1057  CALL ssptrd( 'L', n, vp, sd, se, tau, iinfo )
1058 *
1059  IF( iinfo.NE.0 ) THEN
1060  WRITE( nounit, fmt = 9999 )'SSPTRD(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 sopgtr( 'L', n, vp, tau, u, ldu, work, iinfo )
1073  IF( iinfo.NE.0 ) THEN
1074  WRITE( nounit, fmt = 9999 )'SOPGTR(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 sspt21( 2, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1086  $ work, result( 7 ) )
1087  CALL sspt21( 3, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1088  $ work, result( 8 ) )
1089 *
1090 * Call SSTEQR to compute D1, D2, and Z, do tests.
1091 *
1092 * Compute D1 and Z
1093 *
1094  CALL scopy( n, sd, 1, d1, 1 )
1095  IF( n.GT.0 )
1096  $ CALL scopy( n-1, se, 1, work, 1 )
1097  CALL slaset( 'Full', n, n, zero, one, z, ldu )
1098 *
1099  ntest = 9
1100  CALL ssteqr( 'V', n, d1, work, z, ldu, work( n+1 ), iinfo )
1101  IF( iinfo.NE.0 ) THEN
1102  WRITE( nounit, fmt = 9999 )'SSTEQR(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 scopy( n, sd, 1, d2, 1 )
1116  IF( n.GT.0 )
1117  $ CALL scopy( n-1, se, 1, work, 1 )
1118 *
1119  ntest = 11
1120  CALL ssteqr( 'N', n, d2, work, work( n+1 ), ldu,
1121  $ work( n+1 ), iinfo )
1122  IF( iinfo.NE.0 ) THEN
1123  WRITE( nounit, fmt = 9999 )'SSTEQR(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 scopy( n, sd, 1, d3, 1 )
1137  IF( n.GT.0 )
1138  $ CALL scopy( n-1, se, 1, work, 1 )
1139 *
1140  ntest = 12
1141  CALL ssterf( n, d3, work, iinfo )
1142  IF( iinfo.NE.0 ) THEN
1143  WRITE( nounit, fmt = 9999 )'SSTERF', 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 sstt21( 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 sstech( 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 SPTEQR
1193 * and do tests 14, 15, and 16 .
1194 *
1195  IF( jtype.GT.15 ) THEN
1196 *
1197 * Compute D4 and Z4
1198 *
1199  CALL scopy( n, sd, 1, d4, 1 )
1200  IF( n.GT.0 )
1201  $ CALL scopy( n-1, se, 1, work, 1 )
1202  CALL slaset( 'Full', n, n, zero, one, z, ldu )
1203 *
1204  ntest = 14
1205  CALL spteqr( 'V', n, d4, work, z, ldu, work( n+1 ),
1206  $ iinfo )
1207  IF( iinfo.NE.0 ) THEN
1208  WRITE( nounit, fmt = 9999 )'SPTEQR(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 sstt21( n, 0, sd, se, d4, dumma, z, ldu, work,
1222  $ result( 14 ) )
1223 *
1224 * Compute D5
1225 *
1226  CALL scopy( n, sd, 1, d5, 1 )
1227  IF( n.GT.0 )
1228  $ CALL scopy( n-1, se, 1, work, 1 )
1229 *
1230  ntest = 16
1231  CALL spteqr( 'N', n, d5, work, z, ldu, work( n+1 ),
1232  $ iinfo )
1233  IF( iinfo.NE.0 ) THEN
1234  WRITE( nounit, fmt = 9999 )'SPTEQR(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 SSTEBZ 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 sstebz( '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 )'SSTEBZ(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 sstebz( '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 )'SSTEBZ(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( slarnd( 1, iseed2 ) )
1344  iu = 1 + ( n-1 )*int( slarnd( 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 sstebz( '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 )'SSTEBZ(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 sstebz( '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 )'SSTEBZ(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 = ssxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1413  temp2 = ssxt1( 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 SSTEIN to compute eigenvectors corresponding to
1423 * eigenvalues in WA1. (First call SSTEBZ again, to make sure
1424 * it returns these eigenvalues in the correct order.)
1425 *
1426  ntest = 21
1427  CALL sstebz( '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 )'SSTEBZ(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 sstein( 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 )'SSTEIN', 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 sstt21( n, 0, sd, se, wa1, dumma, z, ldu, work,
1462  $ result( 20 ) )
1463 *
1464 * Call SSTEDC(I) to compute D1 and Z, do tests.
1465 *
1466 * Compute D1 and Z
1467 *
1468  CALL scopy( n, sd, 1, d1, 1 )
1469  IF( n.GT.0 )
1470  $ CALL scopy( n-1, se, 1, work, 1 )
1471  CALL slaset( 'Full', n, n, zero, one, z, ldu )
1472 *
1473  ntest = 22
1474  CALL sstedc( '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 )'SSTEDC(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 sstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1491  $ result( 22 ) )
1492 *
1493 * Call SSTEDC(V) to compute D1 and Z, do tests.
1494 *
1495 * Compute D1 and Z
1496 *
1497  CALL scopy( n, sd, 1, d1, 1 )
1498  IF( n.GT.0 )
1499  $ CALL scopy( n-1, se, 1, work, 1 )
1500  CALL slaset( 'Full', n, n, zero, one, z, ldu )
1501 *
1502  ntest = 24
1503  CALL sstedc( '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 )'SSTEDC(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 sstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1520  $ result( 24 ) )
1521 *
1522 * Call SSTEDC(N) to compute D2, do tests.
1523 *
1524 * Compute D2
1525 *
1526  CALL scopy( n, sd, 1, d2, 1 )
1527  IF( n.GT.0 )
1528  $ CALL scopy( n-1, se, 1, work, 1 )
1529  CALL slaset( 'Full', n, n, zero, one, z, ldu )
1530 *
1531  ntest = 26
1532  CALL sstedc( '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 )'SSTEDC(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 SSTEMR if IEEE compliant
1559 *
1560  IF( ilaenv( 10, 'SSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
1561  $ ilaenv( 11, 'SSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
1562 *
1563 * Call SSTEMR, 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 sstemr( '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 )'SSTEMR(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( slarnd( 1, iseed2 ) )
1605  iu = 1 + ( n-1 )*int( slarnd( 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 sstemr( '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 )'SSTEMR(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 SSTEMR(V,I) to compute D1 and Z, do tests.
1654 *
1655 * Compute D1 and Z
1656 *
1657  CALL scopy( n, sd, 1, d5, 1 )
1658  IF( n.GT.0 )
1659  $ CALL scopy( n-1, se, 1, work, 1 )
1660  CALL slaset( 'Full', n, n, zero, one, z, ldu )
1661 *
1662  IF( srange ) THEN
1663  ntest = 29
1664  il = 1 + ( n-1 )*int( slarnd( 1, iseed2 ) )
1665  iu = 1 + ( n-1 )*int( slarnd( 1, iseed2 ) )
1666  IF( iu.LT.il ) THEN
1667  itemp = iu
1668  iu = il
1669  il = itemp
1670  END IF
1671  CALL sstemr( '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 )'SSTEMR(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 sstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1690  $ m, result( 29 ) )
1691 *
1692 * Call SSTEMR to compute D2, do tests.
1693 *
1694 * Compute D2
1695 *
1696  CALL scopy( n, sd, 1, d5, 1 )
1697  IF( n.GT.0 )
1698  $ CALL scopy( n-1, se, 1, work, 1 )
1699 *
1700  ntest = 31
1701  CALL sstemr( '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 )'SSTEMR(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 SSTEMR(V,V) to compute D1 and Z, do tests.
1733 *
1734 * Compute D1 and Z
1735 *
1736  CALL scopy( n, sd, 1, d5, 1 )
1737  IF( n.GT.0 )
1738  $ CALL scopy( n-1, se, 1, work, 1 )
1739  CALL slaset( '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 sstemr( '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 )'SSTEMR(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 sstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1784  $ m, result( 32 ) )
1785 *
1786 * Call SSTEMR to compute D2, do tests.
1787 *
1788 * Compute D2
1789 *
1790  CALL scopy( n, sd, 1, d5, 1 )
1791  IF( n.GT.0 )
1792  $ CALL scopy( n-1, se, 1, work, 1 )
1793 *
1794  ntest = 34
1795  CALL sstemr( '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 )'SSTEMR(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 SSTEMR(V,A) to compute D1 and Z, do tests.
1835 *
1836 * Compute D1 and Z
1837 *
1838  CALL scopy( n, sd, 1, d5, 1 )
1839  IF( n.GT.0 )
1840  $ CALL scopy( n-1, se, 1, work, 1 )
1841 *
1842  ntest = 35
1843 *
1844  CALL sstemr( '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 )'SSTEMR(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 sstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work, m,
1863  $ result( 35 ) )
1864 *
1865 * Call SSTEMR to compute D2, do tests.
1866 *
1867 * Compute D2
1868 *
1869  CALL scopy( n, sd, 1, d5, 1 )
1870  IF( n.GT.0 )
1871  $ CALL scopy( n-1, se, 1, work, 1 )
1872 *
1873  ntest = 37
1874  CALL sstemr( '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 )'SSTEMR(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 )'SST'
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 slasum( 'SST', nounit, nerrs, ntestt )
1940  RETURN
1941 *
1942  9999 FORMAT( ' SCHKST: ', 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 SCHKST 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 SCHKST for details.', / )
1977 * End of SCHKST
1978 *
1979  END
subroutine sstt21(N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RESULT)
SSTT21
Definition: sstt21.f:129
subroutine sorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
SORGTR
Definition: sorgtr.f:125
subroutine sstt22(N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, LDWORK, RESULT)
SSTT22
Definition: sstt22.f:141
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine ssptrd(UPLO, N, AP, D, E, TAU, INFO)
SSPTRD
Definition: ssptrd.f:152
subroutine sopgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
SOPGTR
Definition: sopgtr.f:116
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:323
subroutine slatmr(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)
SLATMR
Definition: slatmr.f:473
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:176
subroutine spteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SPTEQR
Definition: spteqr.f:147
subroutine ssytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
SSYTRD
Definition: ssytrd.f:194
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:190
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:265
subroutine schkst(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)
SCHKST
Definition: schkst.f:593
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine slasum(TYPE, IOUNIT, IE, NRUN)
SLASUM
Definition: slasum.f:42
subroutine ssyt21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RESULT)
SSYT21
Definition: ssyt21.f:207
subroutine sstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEMR
Definition: sstemr.f:314
subroutine sstech(N, A, B, EIG, TOL, WORK, INFO)
SSTECH
Definition: sstech.f:103
subroutine sspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
SSPT21
Definition: sspt21.f:221
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112