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