LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zchkst2stg.f
Go to the documentation of this file.
1*> \brief \b ZCHKST2STG
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 ZCHKST2STG( 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* DOUBLE PRECISION THRESH
21* ..
22* .. Array Arguments ..
23* LOGICAL DOTYPE( * )
24* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
25* DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ),
26* $ RESULT( * ), RWORK( * ), SD( * ), SE( * ),
27* $ WA1( * ), WA2( * ), WA3( * ), WR( * )
28* COMPLEX*16 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*> ZCHKST2STG 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*> ZHETRD. For that, we call the standard ZHETRD and compute D1 using
44*> DSTEQR, then we call the 2-stage ZHETRD_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 ZCHKST in the next
49*> release when vectors and generation of Q will be implemented.
50*>
51*> ZHETRD factors A as U S U* , where * means conjugate transpose,
52*> S is real symmetric tridiagonal, and U is unitary.
53*> ZHETRD can use either just the lower or just the upper triangle
54*> of A; ZCHKST2STG 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*> ZHPTRD does the same as ZHETRD, except that A and V are stored
60*> in "packed" format.
61*>
62*> ZUNGTR constructs the matrix U from the contents of V and TAU.
63*>
64*> ZUPGTR constructs the matrix U from the contents of VP and TAU.
65*>
66*> ZSTEQR 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*> DSTERF computes D3, the matrix of eigenvalues, by the
72*> PWK method, which does not yield eigenvectors.
73*>
74*> ZPTEQR 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*> DSTEBZ 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*> ZSTEIN computes Y, the eigenvectors of S, given the
86*> eigenvalues.
87*>
88*> ZSTEDC 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 ZHETRD/ZUNGTR or ZHPTRD/ZUPGTR ('V' option). It may
93*> also just compute eigenvalues ('N' option).
94*>
95*> ZSTEMR 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). ZSTEMR
98*> uses the Relatively Robust Representation whenever possible.
99*>
100*> When ZCHKST2STG 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 ) ZHETRD( UPLO='U', ... )
107*>
108*> (2) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='U', ... )
109*>
110*> (3) | A - V S V* | / ( |A| n ulp ) ZHETRD( 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*> ZHETRD_2STAGE("N", "U",....). D1 and D2 are computed
115*> via DSTEQR('N',...)
116*>
117*> (4) | I - UV* | / ( n ulp ) ZUNGTR( 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*> ZHETRD_2STAGE("N", "L",....). D1 and D3 are computed
122*> via DSTEQR('N',...)
123*>
124*> (5-8) Same as 1-4, but for ZHPTRD and ZUPGTR.
125*>
126*> (9) | S - Z D Z* | / ( |S| n ulp ) ZSTEQR('V',...)
127*>
128*> (10) | I - ZZ* | / ( n ulp ) ZSTEQR('V',...)
129*>
130*> (11) | D1 - D2 | / ( |D1| ulp ) ZSTEQR('N',...)
131*>
132*> (12) | D1 - D3 | / ( |D1| ulp ) DSTERF
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*> DSTECH)
138*>
139*> For S positive definite,
140*>
141*> (14) | S - Z4 D4 Z4* | / ( |S| n ulp ) ZPTEQR('V',...)
142*>
143*> (15) | I - Z4 Z4* | / ( n ulp ) ZPTEQR('V',...)
144*>
145*> (16) | D4 - D5 | / ( 100 |D4| ulp ) ZPTEQR('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*> DSTEBZ( 'A', 'E', ...)
153*>
154*> (18) | WA1 - D3 | / ( |D3| ulp ) DSTEBZ( '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*> DSTEBZ( 'I', 'E', ...)
161*>
162*> (20) | S - Y WA1 Y* | / ( |S| n ulp ) DSTEBZ, ZSTEIN
163*>
164*> (21) | I - Y Y* | / ( n ulp ) DSTEBZ, ZSTEIN
165*>
166*> (22) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('I')
167*>
168*> (23) | I - ZZ* | / ( n ulp ) ZSTEDC('I')
169*>
170*> (24) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('V')
171*>
172*> (25) | I - ZZ* | / ( n ulp ) ZSTEDC('V')
173*>
174*> (26) | D1 - D2 | / ( |D1| ulp ) ZSTEDC('V') and
175*> ZSTEDC('N')
176*>
177*> Test 27 is disabled at the moment because ZSTEMR 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*> ZSTEMR('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*> ZSTEMR('V', 'I')
189*>
190*> Tests 29 through 34 are disable at present because ZSTEMR
191*> does not handle partial spectrum requests.
192*>
193*> (29) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'I')
194*>
195*> (30) | I - ZZ* | / ( n ulp ) ZSTEMR('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*> ZSTEMR('N', 'I') vs. CSTEMR('V', 'I')
202*>
203*> (32) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'V')
204*>
205*> (33) | I - ZZ* | / ( n ulp ) ZSTEMR('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*> ZSTEMR('N', 'V') vs. CSTEMR('V', 'V')
212*>
213*> (35) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'A')
214*>
215*> (36) | I - ZZ* | / ( n ulp ) ZSTEMR('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*> ZSTEMR('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*> ZCHKST2STG 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, ZCHKST2STG
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 ZCHKST2STG to continue the same random number
322*> sequence.
323*> \endverbatim
324*>
325*> \param[in] THRESH
326*> \verbatim
327*> THRESH is DOUBLE PRECISION
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*16 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*16 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 DOUBLE PRECISION array of
369*> dimension( max(NN) )
370*> The diagonal of the tridiagonal matrix computed by ZHETRD.
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 DOUBLE PRECISION array of
378*> dimension( max(NN) )
379*> The off-diagonal of the tridiagonal matrix computed by
380*> ZHETRD. 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 DOUBLE PRECISION array of
387*> dimension( max(NN) )
388*> The eigenvalues of A, as computed by ZSTEQR 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 DOUBLE PRECISION array of
396*> dimension( max(NN) )
397*> The eigenvalues of A, as computed by ZSTEQR 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 DOUBLE PRECISION array of
405*> dimension( max(NN) )
406*> The eigenvalues of A, as computed by DSTERF. 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 DOUBLE PRECISION array of
413*> dimension( max(NN) )
414*> The eigenvalues of A, as computed by ZPTEQR(V).
415*> ZPTEQR 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 DOUBLE PRECISION array of
422*> dimension( max(NN) )
423*> The eigenvalues of A, as computed by ZPTEQR(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 DOUBLE PRECISION 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 DSTEBZ.
435*> \endverbatim
436*>
437*> \param[out] WA2
438*> \verbatim
439*> WA2 is DOUBLE PRECISION 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 DSTEBZ.
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 DOUBLE PRECISION 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 DSTEBZ.
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 DOUBLE PRECISION array of
462*> dimension( max(NN) )
463*> All eigenvalues of A, computed to high
464*> absolute accuracy, with different options.
465*> as computed by DSTEBZ.
466*> \endverbatim
467*>
468*> \param[out] U
469*> \verbatim
470*> U is COMPLEX*16 array of
471*> dimension( LDU, max(NN) ).
472*> The unitary matrix computed by ZHETRD + ZUNGTR.
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*16 array of
485*> dimension( LDU, max(NN) ).
486*> The Housholder vectors computed by ZHETRD 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 ZHETRD, 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 ZUNGTR, 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*16 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*16 array of
506*> dimension( max(NN) )
507*> The Householder factors computed by ZHETRD in reducing A
508*> to tridiagonal form.
509*> \endverbatim
510*>
511*> \param[out] Z
512*> \verbatim
513*> Z is COMPLEX*16 array of
514*> dimension( LDU, max(NN) ).
515*> The unitary matrix of eigenvectors computed by ZSTEQR,
516*> ZPTEQR, and ZSTEIN.
517*> \endverbatim
518*>
519*> \param[out] WORK
520*> \verbatim
521*> WORK is COMPLEX*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 ZLATMR, CLATMS, ZHETRD, ZUNGTR, ZSTEQR, DSTERF,
578*> or ZUNMC2 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 complex16_eig
618*
619* =====================================================================
620 SUBROUTINE zchkst2stg( 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 DOUBLE PRECISION THRESH
634* ..
635* .. Array Arguments ..
636 LOGICAL DOTYPE( * )
637 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
638 DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ),
639 $ RESULT( * ), RWORK( * ), SD( * ), SE( * ),
640 $ wa1( * ), wa2( * ), wa3( * ), wr( * )
641 COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ),
642 $ v( ldu, * ), vp( * ), work( * ), z( ldu, * )
643* ..
644*
645* =====================================================================
646*
647* .. Parameters ..
648 DOUBLE PRECISION ZERO, ONE, TWO, EIGHT, TEN, HUN
649 PARAMETER ( ZERO = 0.0d0, one = 1.0d0, two = 2.0d0,
650 $ eight = 8.0d0, ten = 10.0d0, hun = 100.0d0 )
651 COMPLEX*16 CZERO, CONE
652 parameter( czero = ( 0.0d+0, 0.0d+0 ),
653 $ cone = ( 1.0d+0, 0.0d+0 ) )
654 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION DUMMA( 1 )
679* ..
680* .. External Functions ..
681 INTEGER ILAENV
682 DOUBLE PRECISION DLAMCH, DLARND, DSXT1
683 EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1
684* ..
685* .. External Subroutines ..
686 EXTERNAL dcopy, dlasum, dstebz, dstech, dsterf, xerbla,
691* ..
692* .. Intrinsic Functions ..
693 INTRINSIC abs, dble, dconjg, 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, 'ZHETRD', '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( 'ZCHKST2STG', -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 = dlamch( 'Safe minimum' )
756 ovfl = one / unfl
757 ulp = dlamch( 'Epsilon' )*dlamch( '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( dble( 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 / dble( 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 zlaset( '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 zlatms( 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 zlatms( 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 zlatmr( 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 zlatmr( 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 zlatms( 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 zlatms( 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 ) = dconjg( 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 ZHETRD and ZUNGTR to compute S and U from
940* upper triangle.
941*
942 CALL zlacpy( 'U', n, n, a, lda, v, ldu )
943*
944 ntest = 1
945 CALL zhetrd( 'U', n, v, ldu, sd, se, tau, work, lwork,
946 $ iinfo )
947*
948 IF( iinfo.NE.0 ) THEN
949 WRITE( nounit, fmt = 9999 )'ZHETRD(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 zlacpy( 'U', n, n, v, ldu, u, ldu )
961*
962 ntest = 2
963 CALL zungtr( 'U', n, u, ldu, tau, work, lwork, iinfo )
964 IF( iinfo.NE.0 ) THEN
965 WRITE( nounit, fmt = 9999 )'ZUNGTR(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 zhet21( 2, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
979 $ ldu, tau, work, rwork, result( 1 ) )
980 CALL zhet21( 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 dcopy( n, sd, 1, d1, 1 )
991 IF( n.GT.0 )
992 $ CALL dcopy( n-1, se, 1, rwork, 1 )
993*
994 CALL zsteqr( 'N', n, d1, rwork, work, ldu, rwork( n+1 ),
995 $ iinfo )
996 IF( iinfo.NE.0 ) THEN
997 WRITE( nounit, fmt = 9999 )'ZSTEQR(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 dlaset( 'Full', n, 1, zero, zero, sd, n )
1014 CALL dlaset( 'Full', n, 1, zero, zero, se, n )
1015 CALL zlacpy( 'U', n, n, a, lda, v, ldu )
1016 lh = max(1, 4*n)
1017 lw = lwork - lh
1018 CALL zhetrd_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 dcopy( n, sd, 1, d2, 1 )
1024 IF( n.GT.0 )
1025 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1026*
1027 ntest = 3
1028 CALL zsteqr( 'N', n, d2, rwork, work, ldu, rwork( n+1 ),
1029 $ iinfo )
1030 IF( iinfo.NE.0 ) THEN
1031 WRITE( nounit, fmt = 9999 )'ZSTEQR(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 dlaset( 'Full', n, 1, zero, zero, sd, n )
1048 CALL dlaset( 'Full', n, 1, zero, zero, se, n )
1049 CALL zlacpy( 'L', n, n, a, lda, v, ldu )
1050 CALL zhetrd_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 dcopy( n, sd, 1, d3, 1 )
1056 IF( n.GT.0 )
1057 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1058*
1059 ntest = 4
1060 CALL zsteqr( 'N', n, d3, rwork, work, ldu, rwork( n+1 ),
1061 $ iinfo )
1062 IF( iinfo.NE.0 ) THEN
1063 WRITE( nounit, fmt = 9999 )'ZSTEQR(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 ZHPTRD and ZUPGTR to compute S and U from AP
1104*
1105 CALL zcopy( nap, ap, 1, vp, 1 )
1106*
1107 ntest = 5
1108 CALL zhptrd( 'U', n, vp, sd, se, tau, iinfo )
1109*
1110 IF( iinfo.NE.0 ) THEN
1111 WRITE( nounit, fmt = 9999 )'ZHPTRD(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 zupgtr( 'U', n, vp, tau, u, ldu, work, iinfo )
1124 IF( iinfo.NE.0 ) THEN
1125 WRITE( nounit, fmt = 9999 )'ZUPGTR(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 zhpt21( 2, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1139 $ work, rwork, result( 5 ) )
1140 CALL zhpt21( 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 ZHPTRD and ZUPGTR to compute S and U from AP
1154*
1155 CALL zcopy( nap, ap, 1, vp, 1 )
1156*
1157 ntest = 7
1158 CALL zhptrd( 'L', n, vp, sd, se, tau, iinfo )
1159*
1160 IF( iinfo.NE.0 ) THEN
1161 WRITE( nounit, fmt = 9999 )'ZHPTRD(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 zupgtr( 'L', n, vp, tau, u, ldu, work, iinfo )
1174 IF( iinfo.NE.0 ) THEN
1175 WRITE( nounit, fmt = 9999 )'ZUPGTR(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 zhpt21( 2, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1187 $ work, rwork, result( 7 ) )
1188 CALL zhpt21( 3, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1189 $ work, rwork, result( 8 ) )
1190*
1191* Call ZSTEQR to compute D1, D2, and Z, do tests.
1192*
1193* Compute D1 and Z
1194*
1195 CALL dcopy( n, sd, 1, d1, 1 )
1196 IF( n.GT.0 )
1197 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1198 CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1199*
1200 ntest = 9
1201 CALL zsteqr( 'V', n, d1, rwork, z, ldu, rwork( n+1 ),
1202 $ iinfo )
1203 IF( iinfo.NE.0 ) THEN
1204 WRITE( nounit, fmt = 9999 )'ZSTEQR(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 dcopy( n, sd, 1, d2, 1 )
1218 IF( n.GT.0 )
1219 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1220*
1221 ntest = 11
1222 CALL zsteqr( 'N', n, d2, rwork, work, ldu, rwork( n+1 ),
1223 $ iinfo )
1224 IF( iinfo.NE.0 ) THEN
1225 WRITE( nounit, fmt = 9999 )'ZSTEQR(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 dcopy( n, sd, 1, d3, 1 )
1239 IF( n.GT.0 )
1240 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1241*
1242 ntest = 12
1243 CALL dsterf( n, d3, rwork, iinfo )
1244 IF( iinfo.NE.0 ) THEN
1245 WRITE( nounit, fmt = 9999 )'DSTERF', 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 zstt21( 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 dstech( 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 ZPTEQR
1295* and do tests 14, 15, and 16 .
1296*
1297 IF( jtype.GT.15 ) THEN
1298*
1299* Compute D4 and Z4
1300*
1301 CALL dcopy( n, sd, 1, d4, 1 )
1302 IF( n.GT.0 )
1303 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1304 CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1305*
1306 ntest = 14
1307 CALL zpteqr( 'V', n, d4, rwork, z, ldu, rwork( n+1 ),
1308 $ iinfo )
1309 IF( iinfo.NE.0 ) THEN
1310 WRITE( nounit, fmt = 9999 )'ZPTEQR(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 zstt21( n, 0, sd, se, d4, dumma, z, ldu, work,
1324 $ rwork, result( 14 ) )
1325*
1326* Compute D5
1327*
1328 CALL dcopy( n, sd, 1, d5, 1 )
1329 IF( n.GT.0 )
1330 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1331*
1332 ntest = 16
1333 CALL zpteqr( 'N', n, d5, rwork, z, ldu, rwork( n+1 ),
1334 $ iinfo )
1335 IF( iinfo.NE.0 ) THEN
1336 WRITE( nounit, fmt = 9999 )'ZPTEQR(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 DSTEBZ 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 dstebz( '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 )'DSTEBZ(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 dstebz( '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 )'DSTEBZ(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( dlarnd( 1, iseed2 ) )
1446 iu = 1 + ( n-1 )*int( dlarnd( 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 dstebz( '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 )'DSTEBZ(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 dstebz( '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 )'DSTEBZ(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 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1515 temp2 = dsxt1( 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 ZSTEIN to compute eigenvectors corresponding to
1525* eigenvalues in WA1. (First call DSTEBZ again, to make sure
1526* it returns these eigenvalues in the correct order.)
1527*
1528 ntest = 21
1529 CALL dstebz( '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 )'DSTEBZ(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 zstein( 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 )'ZSTEIN', 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 zstt21( n, 0, sd, se, wa1, dumma, z, ldu, work, rwork,
1564 $ result( 20 ) )
1565*
1566* Call ZSTEDC(I) to compute D1 and Z, do tests.
1567*
1568* Compute D1 and Z
1569*
1570 inde = 1
1571 indrwk = inde + n
1572 CALL dcopy( n, sd, 1, d1, 1 )
1573 IF( n.GT.0 )
1574 $ CALL dcopy( n-1, se, 1, rwork( inde ), 1 )
1575 CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1576*
1577 ntest = 22
1578 CALL zstedc( '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 )'ZSTEDC(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 zstt21( n, 0, sd, se, d1, dumma, z, ldu, work, rwork,
1595 $ result( 22 ) )
1596*
1597* Call ZSTEDC(V) to compute D1 and Z, do tests.
1598*
1599* Compute D1 and Z
1600*
1601 CALL dcopy( n, sd, 1, d1, 1 )
1602 IF( n.GT.0 )
1603 $ CALL dcopy( n-1, se, 1, rwork( inde ), 1 )
1604 CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1605*
1606 ntest = 24
1607 CALL zstedc( '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 )'ZSTEDC(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 zstt21( n, 0, sd, se, d1, dumma, z, ldu, work, rwork,
1624 $ result( 24 ) )
1625*
1626* Call ZSTEDC(N) to compute D2, do tests.
1627*
1628* Compute D2
1629*
1630 CALL dcopy( n, sd, 1, d2, 1 )
1631 IF( n.GT.0 )
1632 $ CALL dcopy( n-1, se, 1, rwork( inde ), 1 )
1633 CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1634*
1635 ntest = 26
1636 CALL zstedc( '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 )'ZSTEDC(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 ZSTEMR if IEEE compliant
1663*
1664 IF( ilaenv( 10, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
1665 $ ilaenv( 11, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
1666*
1667* Call ZSTEMR, 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 zstemr( '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 )'ZSTEMR(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( dlarnd( 1, iseed2 ) )
1709 iu = 1 + ( n-1 )*int( dlarnd( 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 zstemr( '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 )'ZSTEMR(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 ZSTEMR(V,I) to compute D1 and Z, do tests.
1757*
1758* Compute D1 and Z
1759*
1760 CALL dcopy( n, sd, 1, d5, 1 )
1761 IF( n.GT.0 )
1762 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1763 CALL zlaset( 'Full', n, n, czero, cone, z, ldu )
1764*
1765 IF( crange ) THEN
1766 ntest = 29
1767 il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1768 iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1769 IF( iu.LT.il ) THEN
1770 itemp = iu
1771 iu = il
1772 il = itemp
1773 END IF
1774 CALL zstemr( '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 )'ZSTEMR(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 ZSTEMR to compute D2, do tests.
1793*
1794* Compute D2
1795*
1796 CALL dcopy( n, sd, 1, d5, 1 )
1797 IF( n.GT.0 )
1798 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1799*
1800 ntest = 31
1801 CALL zstemr( '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 )'ZSTEMR(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 ZSTEMR(V,V) to compute D1 and Z, do tests.
1832*
1833* Compute D1 and Z
1834*
1835 CALL dcopy( n, sd, 1, d5, 1 )
1836 IF( n.GT.0 )
1837 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1838 CALL zlaset( '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 zstemr( '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 )'ZSTEMR(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 zstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1883 $ m, rwork, result( 32 ) )
1884*
1885* Call ZSTEMR to compute D2, do tests.
1886*
1887* Compute D2
1888*
1889 CALL dcopy( n, sd, 1, d5, 1 )
1890 IF( n.GT.0 )
1891 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1892*
1893 ntest = 34
1894 CALL zstemr( '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 )'ZSTEMR(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 ZSTEMR(V,A) to compute D1 and Z, do tests.
1933*
1934* Compute D1 and Z
1935*
1936 CALL dcopy( n, sd, 1, d5, 1 )
1937 IF( n.GT.0 )
1938 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1939*
1940 ntest = 35
1941*
1942 CALL zstemr( '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 )'ZSTEMR(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 zstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work, m,
1961 $ rwork, result( 35 ) )
1962*
1963* Call ZSTEMR to compute D2, do tests.
1964*
1965* Compute D2
1966*
1967 CALL dcopy( n, sd, 1, d5, 1 )
1968 IF( n.GT.0 )
1969 $ CALL dcopy( n-1, se, 1, rwork, 1 )
1970*
1971 ntest = 37
1972 CALL zstemr( '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 )'ZSTEMR(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 )'ZST'
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.0d0 ) 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 dlasum( 'ZST', nounit, nerrs, ntestt )
2042 RETURN
2043*
2044 9999 FORMAT( ' ZCHKST2STG: ', 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 ZCHKST2STG 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, d10.3 )
2079*
2080 9987 FORMAT( / 'Test performed: see ZCHKST2STG for details.', / )
2081*
2082* End of ZCHKST2STG
2083*
2084 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlasum(type, iounit, ie, nrun)
DLASUM
Definition dlasum.f:43
subroutine dstech(n, a, b, eig, tol, work, info)
DSTECH
Definition dstech.f:101
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine zhetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
ZHETRD_2STAGE
subroutine zhetrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
ZHETRD
Definition zhetrd.f:192
subroutine zhptrd(uplo, n, ap, d, e, tau, info)
ZHPTRD
Definition zhptrd.f:151
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
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:110
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zpteqr(compz, n, d, e, z, ldz, work, info)
ZPTEQR
Definition zpteqr.f:145
subroutine dstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
DSTEBZ
Definition dstebz.f:273
subroutine zstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZSTEDC
Definition zstedc.f:206
subroutine zstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
ZSTEIN
Definition zstein.f:182
subroutine zstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
ZSTEMR
Definition zstemr.f:339
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine zungtr(uplo, n, a, lda, tau, work, lwork, info)
ZUNGTR
Definition zungtr.f:123
subroutine zupgtr(uplo, n, ap, tau, q, ldq, work, info)
ZUPGTR
Definition zupgtr.f:114
subroutine zchkst2stg(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)
ZCHKST2STG
Definition zchkst2stg.f:625
subroutine zhet21(itype, uplo, n, kband, a, lda, d, e, u, ldu, v, ldv, tau, work, rwork, result)
ZHET21
Definition zhet21.f:214
subroutine zhpt21(itype, uplo, n, kband, ap, d, e, u, ldu, vp, tau, work, rwork, result)
ZHPT21
Definition zhpt21.f:228
subroutine zlatmr(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)
ZLATMR
Definition zlatmr.f:490
subroutine zlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
ZLATMS
Definition zlatms.f:332
subroutine zstt21(n, kband, ad, ae, sd, se, u, ldu, work, rwork, result)
ZSTT21
Definition zstt21.f:133
subroutine zstt22(n, m, kband, ad, ae, sd, se, u, ldu, work, ldwork, rwork, result)
ZSTT22
Definition zstt22.f:145