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