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