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