OpenRadioss 2025.1.11
OpenRadioss project
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 simlutaneously
389*> with Z. On exit, the eigenvalues in D1 correspond with the
390*> matrix in A.
391*> \endverbatim
392*>
393*> \param[out] D2
394*> \verbatim
395*> D2 is 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, dlabad, dlasum, dstebz, dstech, dsterf,
687 $ xerbla, zcopy, zhet21, zhetrd, zhpt21, zhptrd,
688 $ zlacpy, zlaset, zlatmr, zlatms, zpteqr, zstedc,
689 $ zstemr, zstein, zsteqr, zstt21, zstt22, zungtr,
690 $ zupgtr, zhetrd_2stage, dlaset
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 CALL dlabad( unfl, ovfl )
758 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
759 ulpinv = one / ulp
760 log2ui = int( log( ulpinv ) / log( two ) )
761 rtunfl = sqrt( unfl )
762 rtovfl = sqrt( ovfl )
763*
764* Loop over sizes, types
765*
766 DO 20 i = 1, 4
767 iseed2( i ) = iseed( i )
768 20 CONTINUE
769 nerrs = 0
770 nmats = 0
771*
772 DO 310 jsize = 1, nsizes
773 n = nn( jsize )
774 IF( n.GT.0 ) THEN
775 lgn = int( log( dble( n ) ) / log( two ) )
776 IF( 2**lgn.LT.n )
777 $ lgn = lgn + 1
778 IF( 2**lgn.LT.n )
779 $ lgn = lgn + 1
780 lwedc = 1 + 4*n + 2*n*lgn + 4*n**2
781 lrwedc = 1 + 3*n + 2*n*lgn + 4*n**2
782 liwedc = 6 + 6*n + 5*n*lgn
783 ELSE
784 lwedc = 8
785 lrwedc = 7
786 liwedc = 12
787 END IF
788 nap = ( n*( n+1 ) ) / 2
789 aninv = one / dble( max( 1, n ) )
790*
791 IF( nsizes.NE.1 ) THEN
792 mtypes = min( maxtyp, ntypes )
793 ELSE
794 mtypes = min( maxtyp+1, ntypes )
795 END IF
796*
797 DO 300 jtype = 1, mtypes
798 IF( .NOT.dotype( jtype ) )
799 $ GO TO 300
800 nmats = nmats + 1
801 ntest = 0
802*
803 DO 30 j = 1, 4
804 ioldsd( j ) = iseed( j )
805 30 CONTINUE
806*
807* Compute "A"
808*
809* Control parameters:
810*
811* KMAGN KMODE KTYPE
812* =1 O(1) clustered 1 zero
813* =2 large clustered 2 identity
814* =3 small exponential (none)
815* =4 arithmetic diagonal, (w/ eigenvalues)
816* =5 random log Hermitian, w/ eigenvalues
817* =6 random (none)
818* =7 random diagonal
819* =8 random Hermitian
820* =9 positive definite
821* =10 diagonally dominant tridiagonal
822*
823 IF( mtypes.GT.maxtyp )
824 $ GO TO 100
825*
826 itype = ktype( jtype )
827 imode = kmode( jtype )
828*
829* Compute norm
830*
831 GO TO ( 40, 50, 60 )kmagn( jtype )
832*
833 40 CONTINUE
834 anorm = one
835 GO TO 70
836*
837 50 CONTINUE
838 anorm = ( rtovfl*ulp )*aninv
839 GO TO 70
840*
841 60 CONTINUE
842 anorm = rtunfl*n*ulpinv
843 GO TO 70
844*
845 70 CONTINUE
846*
847 CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
848 iinfo = 0
849 IF( jtype.LE.15 ) THEN
850 cond = ulpinv
851 ELSE
852 cond = ulpinv*aninv / ten
853 END IF
854*
855* Special Matrices -- Identity & Jordan block
856*
857* Zero
858*
859 IF( itype.EQ.1 ) THEN
860 iinfo = 0
861*
862 ELSE IF( itype.EQ.2 ) THEN
863*
864* Identity
865*
866 DO 80 jc = 1, n
867 a( jc, jc ) = anorm
868 80 CONTINUE
869*
870 ELSE IF( itype.EQ.4 ) THEN
871*
872* Diagonal Matrix, [Eigen]values Specified
873*
874 CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
875 $ anorm, 0, 0, 'N', a, lda, work, iinfo )
876*
877*
878 ELSE IF( itype.EQ.5 ) THEN
879*
880* Hermitian, eigenvalues specified
881*
882 CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
883 $ anorm, n, n, 'N', a, lda, work, iinfo )
884*
885 ELSE IF( itype.EQ.7 ) THEN
886*
887* Diagonal, random eigenvalues
888*
889 CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one, cone,
890 $ 'T', 'N', work( n+1 ), 1, one,
891 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
892 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
893*
894 ELSE IF( itype.EQ.8 ) THEN
895*
896* Hermitian, random eigenvalues
897*
898 CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one, cone,
899 $ 'T', 'N', work( n+1 ), 1, one,
900 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
901 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
902*
903 ELSE IF( itype.EQ.9 ) THEN
904*
905* Positive definite, eigenvalues specified.
906*
907 CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode, cond,
908 $ anorm, n, n, 'N', a, lda, work, iinfo )
909*
910 ELSE IF( itype.EQ.10 ) THEN
911*
912* Positive definite tridiagonal, eigenvalues specified.
913*
914 CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode, cond,
915 $ anorm, 1, 1, 'N', a, lda, work, iinfo )
916 DO 90 i = 2, n
917 temp1 = abs( a( i-1, i ) )
918 temp2 = sqrt( abs( a( i-1, i-1 )*a( i, i ) ) )
919 IF( temp1.GT.half*temp2 ) THEN
920 a( i-1, i ) = a( i-1, i )*
921 $ ( half*temp2 / ( unfl+temp1 ) )
922 a( i, i-1 ) = dconjg( a( i-1, i ) )
923 END IF
924 90 CONTINUE
925*
926 ELSE
927*
928 iinfo = 1
929 END IF
930*
931 IF( iinfo.NE.0 ) THEN
932 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
933 $ ioldsd
934 info = abs( iinfo )
935 RETURN
936 END IF
937*
938 100 CONTINUE
939*
940* Call ZHETRD and ZUNGTR to compute S and U from
941* upper triangle.
942*
943 CALL zlacpy( 'U', n, n, a, lda, v, ldu )
944*
945 ntest = 1
946 CALL zhetrd( 'U', n, v, ldu, sd, se, tau, work, lwork,
947 $ iinfo )
948*
949 IF( iinfo.NE.0 ) THEN
950 WRITE( nounit, fmt = 9999 )'ZHETRD(U)', iinfo, n, jtype,
951 $ ioldsd
952 info = abs( iinfo )
953 IF( iinfo.LT.0 ) THEN
954 RETURN
955 ELSE
956 result( 1 ) = ulpinv
957 GO TO 280
958 END IF
959 END IF
960*
961 CALL zlacpy( 'U', n, n, v, ldu, u, ldu )
962*
963 ntest = 2
964 CALL zungtr( 'U', n, u, ldu, tau, work, lwork, iinfo )
965 IF( iinfo.NE.0 ) THEN
966 WRITE( nounit, fmt = 9999 )'ZUNGTR(U)', iinfo, n, jtype,
967 $ ioldsd
968 info = abs( iinfo )
969 IF( iinfo.LT.0 ) THEN
970 RETURN
971 ELSE
972 result( 2 ) = ulpinv
973 GO TO 280
974 END IF
975 END IF
976*
977* Do tests 1 and 2
978*
979 CALL zhet21( 2, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
980 $ ldu, tau, work, rwork, result( 1 ) )
981 CALL zhet21( 3, 'upper', N, 1, A, LDA, SD, SE, U, LDU, V,
982 $ LDU, TAU, WORK, RWORK, RESULT( 2 ) )
983*
984* Compute D1 the eigenvalues resulting from the tridiagonal
985* form using the standard 1-stage algorithm and use it as a
986* reference to compare with the 2-stage technique
987*
988* Compute D1 from the 1-stage and used as reference for the
989* 2-stage
990*
991 CALL DCOPY( N, SD, 1, D1, 1 )
992.GT. IF( N0 )
993 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
994*
995 CALL ZSTEQR( 'n', N, D1, RWORK, WORK, LDU, RWORK( N+1 ),
996 $ IINFO )
997.NE. IF( IINFO0 ) THEN
998 WRITE( NOUNIT, FMT = 9999 )'zsteqr(n)', IINFO, N, JTYPE,
999 $ IOLDSD
1000 INFO = ABS( IINFO )
1001.LT. IF( IINFO0 ) THEN
1002 RETURN
1003 ELSE
1004 RESULT( 3 ) = ULPINV
1005 GO TO 280
1006 END IF
1007 END IF
1008*
1009* 2-STAGE TRD Upper case is used to compute D2.
1010* Note to set SD and SE to zero to be sure not reusing
1011* the one from above. Compare it with D1 computed
1012* using the 1-stage.
1013*
1014 CALL DLASET( 'full', N, 1, ZERO, ZERO, SD, N )
1015 CALL DLASET( 'full', N, 1, ZERO, ZERO, SE, N )
1016 CALL ZLACPY( 'u', N, N, A, LDA, V, LDU )
1017 LH = MAX(1, 4*N)
1018 LW = LWORK - LH
1019 CALL ZHETRD_2STAGE( 'n', "U", N, V, LDU, SD, SE, TAU,
1020 $ WORK, LH, WORK( LH+1 ), LW, IINFO )
1021*
1022* Compute D2 from the 2-stage Upper case
1023*
1024 CALL DCOPY( N, SD, 1, D2, 1 )
1025.GT. IF( N0 )
1026 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1027*
1028 NTEST = 3
1029 CALL ZSTEQR( 'n', N, D2, RWORK, WORK, LDU, RWORK( N+1 ),
1030 $ IINFO )
1031.NE. IF( IINFO0 ) THEN
1032 WRITE( NOUNIT, FMT = 9999 )'zsteqr(n)', IINFO, N, JTYPE,
1033 $ IOLDSD
1034 INFO = ABS( IINFO )
1035.LT. IF( IINFO0 ) THEN
1036 RETURN
1037 ELSE
1038 RESULT( 3 ) = ULPINV
1039 GO TO 280
1040 END IF
1041 END IF
1042*
1043* 2-STAGE TRD Lower case is used to compute D3.
1044* Note to set SD and SE to zero to be sure not reusing
1045* the one from above. Compare it with D1 computed
1046* using the 1-stage.
1047*
1048 CALL DLASET( 'full', N, 1, ZERO, ZERO, SD, N )
1049 CALL DLASET( 'full', N, 1, ZERO, ZERO, SE, N )
1050 CALL ZLACPY( 'l', N, N, A, LDA, V, LDU )
1051 CALL ZHETRD_2STAGE( 'n', "L", N, V, LDU, SD, SE, TAU,
1052 $ WORK, LH, WORK( LH+1 ), LW, IINFO )
1053*
1054* Compute D3 from the 2-stage Upper case
1055*
1056 CALL DCOPY( N, SD, 1, D3, 1 )
1057.GT. IF( N0 )
1058 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1059*
1060 NTEST = 4
1061 CALL ZSTEQR( 'n', N, D3, RWORK, WORK, LDU, RWORK( N+1 ),
1062 $ IINFO )
1063.NE. IF( IINFO0 ) THEN
1064 WRITE( NOUNIT, FMT = 9999 )'zsteqr(n)', IINFO, N, JTYPE,
1065 $ IOLDSD
1066 INFO = ABS( IINFO )
1067.LT. IF( IINFO0 ) THEN
1068 RETURN
1069 ELSE
1070 RESULT( 4 ) = ULPINV
1071 GO TO 280
1072 END IF
1073 END IF
1074*
1075* Do Tests 3 and 4 which are similar to 11 and 12 but with the
1076* D1 computed using the standard 1-stage reduction as reference
1077*
1078 NTEST = 4
1079 TEMP1 = ZERO
1080 TEMP2 = ZERO
1081 TEMP3 = ZERO
1082 TEMP4 = ZERO
1083*
1084 DO 151 J = 1, N
1085 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1086 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1087 TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) )
1088 TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) )
1089 151 CONTINUE
1090*
1091 RESULT( 3 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1092 RESULT( 4 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) )
1093*
1094* Store the upper triangle of A in AP
1095*
1096 I = 0
1097 DO 120 JC = 1, N
1098 DO 110 JR = 1, JC
1099 I = I + 1
1100 AP( I ) = A( JR, JC )
1101 110 CONTINUE
1102 120 CONTINUE
1103*
1104* Call ZHPTRD and ZUPGTR to compute S and U from AP
1105*
1106 CALL ZCOPY( NAP, AP, 1, VP, 1 )
1107*
1108 NTEST = 5
1109 CALL ZHPTRD( 'u', N, VP, SD, SE, TAU, IINFO )
1110*
1111.NE. IF( IINFO0 ) THEN
1112 WRITE( NOUNIT, FMT = 9999 )'zhptrd(u)', IINFO, N, JTYPE,
1113 $ IOLDSD
1114 INFO = ABS( IINFO )
1115.LT. IF( IINFO0 ) THEN
1116 RETURN
1117 ELSE
1118 RESULT( 5 ) = ULPINV
1119 GO TO 280
1120 END IF
1121 END IF
1122*
1123 NTEST = 6
1124 CALL ZUPGTR( 'u', N, VP, TAU, U, LDU, WORK, IINFO )
1125.NE. IF( IINFO0 ) THEN
1126 WRITE( NOUNIT, FMT = 9999 )'zupgtr(u)', IINFO, N, JTYPE,
1127 $ IOLDSD
1128 INFO = ABS( IINFO )
1129.LT. IF( IINFO0 ) THEN
1130 RETURN
1131 ELSE
1132 RESULT( 6 ) = ULPINV
1133 GO TO 280
1134 END IF
1135 END IF
1136*
1137* Do tests 5 and 6
1138*
1139 CALL ZHPT21( 2, 'upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1140 $ WORK, RWORK, RESULT( 5 ) )
1141 CALL ZHPT21( 3, 'upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1142 $ WORK, RWORK, RESULT( 6 ) )
1143*
1144* Store the lower triangle of A in AP
1145*
1146 I = 0
1147 DO 140 JC = 1, N
1148 DO 130 JR = JC, N
1149 I = I + 1
1150 AP( I ) = A( JR, JC )
1151 130 CONTINUE
1152 140 CONTINUE
1153*
1154* Call ZHPTRD and ZUPGTR to compute S and U from AP
1155*
1156 CALL ZCOPY( NAP, AP, 1, VP, 1 )
1157*
1158 NTEST = 7
1159 CALL ZHPTRD( 'l', N, VP, SD, SE, TAU, IINFO )
1160*
1161.NE. IF( IINFO0 ) THEN
1162 WRITE( NOUNIT, FMT = 9999 )'zhptrd(l)', IINFO, N, JTYPE,
1163 $ IOLDSD
1164 INFO = ABS( IINFO )
1165.LT. IF( IINFO0 ) THEN
1166 RETURN
1167 ELSE
1168 RESULT( 7 ) = ULPINV
1169 GO TO 280
1170 END IF
1171 END IF
1172*
1173 NTEST = 8
1174 CALL ZUPGTR( 'l', N, VP, TAU, U, LDU, WORK, IINFO )
1175.NE. IF( IINFO0 ) THEN
1176 WRITE( NOUNIT, FMT = 9999 )'zupgtr(l)', IINFO, N, JTYPE,
1177 $ IOLDSD
1178 INFO = ABS( IINFO )
1179.LT. IF( IINFO0 ) THEN
1180 RETURN
1181 ELSE
1182 RESULT( 8 ) = ULPINV
1183 GO TO 280
1184 END IF
1185 END IF
1186*
1187 CALL ZHPT21( 2, 'lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1188 $ WORK, RWORK, RESULT( 7 ) )
1189 CALL ZHPT21( 3, 'lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1190 $ WORK, RWORK, RESULT( 8 ) )
1191*
1192* Call ZSTEQR to compute D1, D2, and Z, do tests.
1193*
1194* Compute D1 and Z
1195*
1196 CALL DCOPY( N, SD, 1, D1, 1 )
1197.GT. IF( N0 )
1198 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1199 CALL ZLASET( 'full', N, N, CZERO, CONE, Z, LDU )
1200*
1201 NTEST = 9
1202 CALL ZSTEQR( 'v', N, D1, RWORK, Z, LDU, RWORK( N+1 ),
1203 $ IINFO )
1204.NE. IF( IINFO0 ) THEN
1205 WRITE( NOUNIT, FMT = 9999 )'zsteqr(v)', IINFO, N, JTYPE,
1206 $ IOLDSD
1207 INFO = ABS( IINFO )
1208.LT. IF( IINFO0 ) THEN
1209 RETURN
1210 ELSE
1211 RESULT( 9 ) = ULPINV
1212 GO TO 280
1213 END IF
1214 END IF
1215*
1216* Compute D2
1217*
1218 CALL DCOPY( N, SD, 1, D2, 1 )
1219.GT. IF( N0 )
1220 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1221*
1222 NTEST = 11
1223 CALL ZSTEQR( 'n', N, D2, RWORK, WORK, LDU, RWORK( N+1 ),
1224 $ IINFO )
1225.NE. IF( IINFO0 ) THEN
1226 WRITE( NOUNIT, FMT = 9999 )'zsteqr(n)', IINFO, N, JTYPE,
1227 $ IOLDSD
1228 INFO = ABS( IINFO )
1229.LT. IF( IINFO0 ) THEN
1230 RETURN
1231 ELSE
1232 RESULT( 11 ) = ULPINV
1233 GO TO 280
1234 END IF
1235 END IF
1236*
1237* Compute D3 (using PWK method)
1238*
1239 CALL DCOPY( N, SD, 1, D3, 1 )
1240.GT. IF( N0 )
1241 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1242*
1243 NTEST = 12
1244 CALL DSTERF( N, D3, RWORK, IINFO )
1245.NE. IF( IINFO0 ) THEN
1246 WRITE( NOUNIT, FMT = 9999 )'dsterf', IINFO, N, JTYPE,
1247 $ IOLDSD
1248 INFO = ABS( IINFO )
1249.LT. IF( IINFO0 ) THEN
1250 RETURN
1251 ELSE
1252 RESULT( 12 ) = ULPINV
1253 GO TO 280
1254 END IF
1255 END IF
1256*
1257* Do Tests 9 and 10
1258*
1259 CALL ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK,
1260 $ RESULT( 9 ) )
1261*
1262* Do Tests 11 and 12
1263*
1264 TEMP1 = ZERO
1265 TEMP2 = ZERO
1266 TEMP3 = ZERO
1267 TEMP4 = ZERO
1268*
1269 DO 150 J = 1, N
1270 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1271 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1272 TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) )
1273 TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) )
1274 150 CONTINUE
1275*
1276 RESULT( 11 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1277 RESULT( 12 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) )
1278*
1279* Do Test 13 -- Sturm Sequence Test of Eigenvalues
1280* Go up by factors of two until it succeeds
1281*
1282 NTEST = 13
1283 TEMP1 = THRESH*( HALF-ULP )
1284*
1285 DO 160 J = 0, LOG2UI
1286 CALL DSTECH( N, SD, SE, D1, TEMP1, RWORK, IINFO )
1287.EQ. IF( IINFO0 )
1288 $ GO TO 170
1289 TEMP1 = TEMP1*TWO
1290 160 CONTINUE
1291*
1292 170 CONTINUE
1293 RESULT( 13 ) = TEMP1
1294*
1295* For positive definite matrices ( JTYPE.GT.15 ) call ZPTEQR
1296* and do tests 14, 15, and 16 .
1297*
1298.GT. IF( JTYPE15 ) THEN
1299*
1300* Compute D4 and Z4
1301*
1302 CALL DCOPY( N, SD, 1, D4, 1 )
1303.GT. IF( N0 )
1304 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1305 CALL ZLASET( 'full', N, N, CZERO, CONE, Z, LDU )
1306*
1307 NTEST = 14
1308 CALL ZPTEQR( 'v', N, D4, RWORK, Z, LDU, RWORK( N+1 ),
1309 $ IINFO )
1310.NE. IF( IINFO0 ) THEN
1311 WRITE( NOUNIT, FMT = 9999 )'zpteqr(v)', IINFO, N,
1312 $ JTYPE, IOLDSD
1313 INFO = ABS( IINFO )
1314.LT. IF( IINFO0 ) THEN
1315 RETURN
1316 ELSE
1317 RESULT( 14 ) = ULPINV
1318 GO TO 280
1319 END IF
1320 END IF
1321*
1322* Do Tests 14 and 15
1323*
1324 CALL ZSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK,
1325 $ RWORK, RESULT( 14 ) )
1326*
1327* Compute D5
1328*
1329 CALL DCOPY( N, SD, 1, D5, 1 )
1330.GT. IF( N0 )
1331 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1332*
1333 NTEST = 16
1334 CALL ZPTEQR( 'n', N, D5, RWORK, Z, LDU, RWORK( N+1 ),
1335 $ IINFO )
1336.NE. IF( IINFO0 ) THEN
1337 WRITE( NOUNIT, FMT = 9999 )'zpteqr(n)', IINFO, N,
1338 $ JTYPE, IOLDSD
1339 INFO = ABS( IINFO )
1340.LT. IF( IINFO0 ) THEN
1341 RETURN
1342 ELSE
1343 RESULT( 16 ) = ULPINV
1344 GO TO 280
1345 END IF
1346 END IF
1347*
1348* Do Test 16
1349*
1350 TEMP1 = ZERO
1351 TEMP2 = ZERO
1352 DO 180 J = 1, N
1353 TEMP1 = MAX( TEMP1, ABS( D4( J ) ), ABS( D5( J ) ) )
1354 TEMP2 = MAX( TEMP2, ABS( D4( J )-D5( J ) ) )
1355 180 CONTINUE
1356*
1357 RESULT( 16 ) = TEMP2 / MAX( UNFL,
1358 $ HUN*ULP*MAX( TEMP1, TEMP2 ) )
1359 ELSE
1360 RESULT( 14 ) = ZERO
1361 RESULT( 15 ) = ZERO
1362 RESULT( 16 ) = ZERO
1363 END IF
1364*
1365* Call DSTEBZ with different options and do tests 17-18.
1366*
1367* If S is positive definite and diagonally dominant,
1368* ask for all eigenvalues with high relative accuracy.
1369*
1370 VL = ZERO
1371 VU = ZERO
1372 IL = 0
1373 IU = 0
1374.EQ. IF( JTYPE21 ) THEN
1375 NTEST = 17
1376 ABSTOL = UNFL + UNFL
1377 CALL DSTEBZ( 'a', 'e', N, VL, VU, IL, IU, ABSTOL, SD, SE,
1378 $ M, NSPLIT, WR, IWORK( 1 ), IWORK( N+1 ),
1379 $ RWORK, IWORK( 2*N+1 ), IINFO )
1380.NE. IF( IINFO0 ) THEN
1381 WRITE( NOUNIT, FMT = 9999 )'dstebz(a,rel)', IINFO, N,
1382 $ JTYPE, IOLDSD
1383 INFO = ABS( IINFO )
1384.LT. IF( IINFO0 ) THEN
1385 RETURN
1386 ELSE
1387 RESULT( 17 ) = ULPINV
1388 GO TO 280
1389 END IF
1390 END IF
1391*
1392* Do test 17
1393*
1394 TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
1395 $ ( ONE-HALF )**4
1396*
1397 TEMP1 = ZERO
1398 DO 190 J = 1, N
1399 TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
1400 $ ( ABSTOL+ABS( D4( J ) ) ) )
1401 190 CONTINUE
1402*
1403 RESULT( 17 ) = TEMP1 / TEMP2
1404 ELSE
1405 RESULT( 17 ) = ZERO
1406 END IF
1407*
1408* Now ask for all eigenvalues with high absolute accuracy.
1409*
1410 NTEST = 18
1411 ABSTOL = UNFL + UNFL
1412 CALL DSTEBZ( 'a', 'e', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
1413 $ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), RWORK,
1414 $ IWORK( 2*N+1 ), IINFO )
1415.NE. IF( IINFO0 ) THEN
1416 WRITE( NOUNIT, FMT = 9999 )'dstebz(a)', IINFO, N, JTYPE,
1417 $ IOLDSD
1418 INFO = ABS( IINFO )
1419.LT. IF( IINFO0 ) THEN
1420 RETURN
1421 ELSE
1422 RESULT( 18 ) = ULPINV
1423 GO TO 280
1424 END IF
1425 END IF
1426*
1427* Do test 18
1428*
1429 TEMP1 = ZERO
1430 TEMP2 = ZERO
1431 DO 200 J = 1, N
1432 TEMP1 = MAX( TEMP1, ABS( D3( J ) ), ABS( WA1( J ) ) )
1433 TEMP2 = MAX( TEMP2, ABS( D3( J )-WA1( J ) ) )
1434 200 CONTINUE
1435*
1436 RESULT( 18 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1437*
1438* Choose random values for IL and IU, and ask for the
1439* IL-th through IU-th eigenvalues.
1440*
1441 NTEST = 19
1442.LE. IF( N1 ) THEN
1443 IL = 1
1444 IU = N
1445 ELSE
1446 IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1447 IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1448.LT. IF( IUIL ) THEN
1449 ITEMP = IU
1450 IU = IL
1451 IL = ITEMP
1452 END IF
1453 END IF
1454*
1455 CALL DSTEBZ( 'i', 'e', N, VL, VU, IL, IU, ABSTOL, SD, SE,
1456 $ M2, NSPLIT, WA2, IWORK( 1 ), IWORK( N+1 ),
1457 $ RWORK, IWORK( 2*N+1 ), IINFO )
1458.NE. IF( IINFO0 ) THEN
1459 WRITE( NOUNIT, FMT = 9999 )'dstebz(i)', IINFO, N, JTYPE,
1460 $ IOLDSD
1461 INFO = ABS( IINFO )
1462.LT. IF( IINFO0 ) THEN
1463 RETURN
1464 ELSE
1465 RESULT( 19 ) = ULPINV
1466 GO TO 280
1467 END IF
1468 END IF
1469*
1470* Determine the values VL and VU of the IL-th and IU-th
1471* eigenvalues and ask for all eigenvalues in this range.
1472*
1473.GT. IF( N0 ) THEN
1474.NE. IF( IL1 ) THEN
1475 VL = WA1( IL ) - MAX( HALF*( WA1( IL )-WA1( IL-1 ) ),
1476 $ ULP*ANORM, TWO*RTUNFL )
1477 ELSE
1478 VL = WA1( 1 ) - MAX( HALF*( WA1( N )-WA1( 1 ) ),
1479 $ ULP*ANORM, TWO*RTUNFL )
1480 END IF
1481.NE. IF( IUN ) THEN
1482 VU = WA1( IU ) + MAX( HALF*( WA1( IU+1 )-WA1( IU ) ),
1483 $ ULP*ANORM, TWO*RTUNFL )
1484 ELSE
1485 VU = WA1( N ) + MAX( HALF*( WA1( N )-WA1( 1 ) ),
1486 $ ULP*ANORM, TWO*RTUNFL )
1487 END IF
1488 ELSE
1489 VL = ZERO
1490 VU = ONE
1491 END IF
1492*
1493 CALL DSTEBZ( 'v', 'e', N, VL, VU, IL, IU, ABSTOL, SD, SE,
1494 $ M3, NSPLIT, WA3, IWORK( 1 ), IWORK( N+1 ),
1495 $ RWORK, IWORK( 2*N+1 ), IINFO )
1496.NE. IF( IINFO0 ) THEN
1497 WRITE( NOUNIT, FMT = 9999 )'dstebz(v)', IINFO, N, JTYPE,
1498 $ IOLDSD
1499 INFO = ABS( IINFO )
1500.LT. IF( IINFO0 ) THEN
1501 RETURN
1502 ELSE
1503 RESULT( 19 ) = ULPINV
1504 GO TO 280
1505 END IF
1506 END IF
1507*
1508.EQ..AND..NE. IF( M30 N0 ) THEN
1509 RESULT( 19 ) = ULPINV
1510 GO TO 280
1511 END IF
1512*
1513* Do test 19
1514*
1515 TEMP1 = DSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
1516 TEMP2 = DSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL )
1517.GT. IF( N0 ) THEN
1518 TEMP3 = MAX( ABS( WA1( N ) ), ABS( WA1( 1 ) ) )
1519 ELSE
1520 TEMP3 = ZERO
1521 END IF
1522*
1523 RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP )
1524*
1525* Call ZSTEIN to compute eigenvectors corresponding to
1526* eigenvalues in WA1. (First call DSTEBZ again, to make sure
1527* it returns these eigenvalues in the correct order.)
1528*
1529 NTEST = 21
1530 CALL DSTEBZ( 'a', 'b', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
1531 $ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), RWORK,
1532 $ IWORK( 2*N+1 ), IINFO )
1533.NE. IF( IINFO0 ) THEN
1534 WRITE( NOUNIT, FMT = 9999 )'dstebz(a,b)', IINFO, N,
1535 $ JTYPE, IOLDSD
1536 INFO = ABS( IINFO )
1537.LT. IF( IINFO0 ) THEN
1538 RETURN
1539 ELSE
1540 RESULT( 20 ) = ULPINV
1541 RESULT( 21 ) = ULPINV
1542 GO TO 280
1543 END IF
1544 END IF
1545*
1546 CALL ZSTEIN( N, SD, SE, M, WA1, IWORK( 1 ), IWORK( N+1 ), Z,
1547 $ LDU, RWORK, IWORK( 2*N+1 ), IWORK( 3*N+1 ),
1548 $ IINFO )
1549.NE. IF( IINFO0 ) THEN
1550 WRITE( NOUNIT, FMT = 9999 )'zstein', IINFO, N, JTYPE,
1551 $ IOLDSD
1552 INFO = ABS( IINFO )
1553.LT. IF( IINFO0 ) THEN
1554 RETURN
1555 ELSE
1556 RESULT( 20 ) = ULPINV
1557 RESULT( 21 ) = ULPINV
1558 GO TO 280
1559 END IF
1560 END IF
1561*
1562* Do tests 20 and 21
1563*
1564 CALL ZSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK, RWORK,
1565 $ RESULT( 20 ) )
1566*
1567* Call ZSTEDC(I) to compute D1 and Z, do tests.
1568*
1569* Compute D1 and Z
1570*
1571 INDE = 1
1572 INDRWK = INDE + N
1573 CALL DCOPY( N, SD, 1, D1, 1 )
1574.GT. IF( N0 )
1575 $ CALL DCOPY( N-1, SE, 1, RWORK( INDE ), 1 )
1576 CALL ZLASET( 'full', N, N, CZERO, CONE, Z, LDU )
1577*
1578 NTEST = 22
1579 CALL ZSTEDC( 'i', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC,
1580 $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO )
1581.NE. IF( IINFO0 ) THEN
1582 WRITE( NOUNIT, FMT = 9999 )'zstedc(i)', IINFO, N, JTYPE,
1583 $ IOLDSD
1584 INFO = ABS( IINFO )
1585.LT. IF( IINFO0 ) THEN
1586 RETURN
1587 ELSE
1588 RESULT( 22 ) = ULPINV
1589 GO TO 280
1590 END IF
1591 END IF
1592*
1593* Do Tests 22 and 23
1594*
1595 CALL ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK,
1596 $ RESULT( 22 ) )
1597*
1598* Call ZSTEDC(V) to compute D1 and Z, do tests.
1599*
1600* Compute D1 and Z
1601*
1602 CALL DCOPY( N, SD, 1, D1, 1 )
1603.GT. IF( N0 )
1604 $ CALL DCOPY( N-1, SE, 1, RWORK( INDE ), 1 )
1605 CALL ZLASET( 'full', N, N, CZERO, CONE, Z, LDU )
1606*
1607 NTEST = 24
1608 CALL ZSTEDC( 'v', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC,
1609 $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO )
1610.NE. IF( IINFO0 ) THEN
1611 WRITE( NOUNIT, FMT = 9999 )'zstedc(v)', IINFO, N, JTYPE,
1612 $ IOLDSD
1613 INFO = ABS( IINFO )
1614.LT. IF( IINFO0 ) THEN
1615 RETURN
1616 ELSE
1617 RESULT( 24 ) = ULPINV
1618 GO TO 280
1619 END IF
1620 END IF
1621*
1622* Do Tests 24 and 25
1623*
1624 CALL ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK,
1625 $ RESULT( 24 ) )
1626*
1627* Call ZSTEDC(N) to compute D2, do tests.
1628*
1629* Compute D2
1630*
1631 CALL DCOPY( N, SD, 1, D2, 1 )
1632.GT. IF( N0 )
1633 $ CALL DCOPY( N-1, SE, 1, RWORK( INDE ), 1 )
1634 CALL ZLASET( 'full', N, N, CZERO, CONE, Z, LDU )
1635*
1636 NTEST = 26
1637 CALL ZSTEDC( 'n', N, D2, RWORK( INDE ), Z, LDU, WORK, LWEDC,
1638 $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO )
1639.NE. IF( IINFO0 ) THEN
1640 WRITE( NOUNIT, FMT = 9999 )'zstedc(n)', IINFO, N, JTYPE,
1641 $ IOLDSD
1642 INFO = ABS( IINFO )
1643.LT. IF( IINFO0 ) THEN
1644 RETURN
1645 ELSE
1646 RESULT( 26 ) = ULPINV
1647 GO TO 280
1648 END IF
1649 END IF
1650*
1651* Do Test 26
1652*
1653 TEMP1 = ZERO
1654 TEMP2 = ZERO
1655*
1656 DO 210 J = 1, N
1657 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1658 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1659 210 CONTINUE
1660*
1661 RESULT( 26 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1662*
1663* Only test ZSTEMR if IEEE compliant
1664*
1665 IF( ILAENV( 10, 'zstemr', 'va.EQ..AND.', 1, 0, 0, 0 )1
1666 $ ILAENV( 11, 'zstemr', 'va.EQ.', 1, 0, 0, 0 )1 ) THEN
1667*
1668* Call ZSTEMR, do test 27 (relative eigenvalue accuracy)
1669*
1670* If S is positive definite and diagonally dominant,
1671* ask for all eigenvalues with high relative accuracy.
1672*
1673 VL = ZERO
1674 VU = ZERO
1675 IL = 0
1676 IU = 0
1677.EQ..AND. IF( JTYPE21 CREL ) THEN
1678 NTEST = 27
1679 ABSTOL = UNFL + UNFL
1680 CALL ZSTEMR( 'v', 'a', N, SD, SE, VL, VU, IL, IU,
1681 $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
1682 $ RWORK, LRWORK, IWORK( 2*N+1 ), LWORK-2*N,
1683 $ IINFO )
1684.NE. IF( IINFO0 ) THEN
1685 WRITE( NOUNIT, FMT = 9999 )'zstemr(v,a,rel)',
1686 $ IINFO, N, JTYPE, IOLDSD
1687 INFO = ABS( IINFO )
1688.LT. IF( IINFO0 ) THEN
1689 RETURN
1690 ELSE
1691 RESULT( 27 ) = ULPINV
1692 GO TO 270
1693 END IF
1694 END IF
1695*
1696* Do test 27
1697*
1698 TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
1699 $ ( ONE-HALF )**4
1700*
1701 TEMP1 = ZERO
1702 DO 220 J = 1, N
1703 TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
1704 $ ( ABSTOL+ABS( D4( J ) ) ) )
1705 220 CONTINUE
1706*
1707 RESULT( 27 ) = TEMP1 / TEMP2
1708*
1709 IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1710 IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1711.LT. IF( IUIL ) THEN
1712 ITEMP = IU
1713 IU = IL
1714 IL = ITEMP
1715 END IF
1716*
1717 IF( CRANGE ) THEN
1718 NTEST = 28
1719 ABSTOL = UNFL + UNFL
1720 CALL ZSTEMR( 'v', 'i', N, SD, SE, VL, VU, IL, IU,
1721 $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
1722 $ RWORK, LRWORK, IWORK( 2*N+1 ),
1723 $ LWORK-2*N, IINFO )
1724*
1725.NE. IF( IINFO0 ) THEN
1726 WRITE( NOUNIT, FMT = 9999 )'zstemr(v,i,rel)',
1727 $ IINFO, N, JTYPE, IOLDSD
1728 INFO = ABS( IINFO )
1729.LT. IF( IINFO0 ) THEN
1730 RETURN
1731 ELSE
1732 RESULT( 28 ) = ULPINV
1733 GO TO 270
1734 END IF
1735 END IF
1736*
1737* Do test 28
1738*
1739 TEMP2 = TWO*( TWO*N-ONE )*ULP*
1740 $ ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4
1741*
1742 TEMP1 = ZERO
1743 DO 230 J = IL, IU
1744 TEMP1 = MAX( TEMP1, ABS( WR( J-IL+1 )-D4( N-J+
1745 $ 1 ) ) / ( ABSTOL+ABS( WR( J-IL+1 ) ) ) )
1746 230 CONTINUE
1747*
1748 RESULT( 28 ) = TEMP1 / TEMP2
1749 ELSE
1750 RESULT( 28 ) = ZERO
1751 END IF
1752 ELSE
1753 RESULT( 27 ) = ZERO
1754 RESULT( 28 ) = ZERO
1755 END IF
1756*
1757* Call ZSTEMR(V,I) to compute D1 and Z, do tests.
1758*
1759* Compute D1 and Z
1760*
1761 CALL DCOPY( N, SD, 1, D5, 1 )
1762.GT. IF( N0 )
1763 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1764 CALL ZLASET( 'full', N, N, CZERO, CONE, Z, LDU )
1765*
1766 IF( CRANGE ) THEN
1767 NTEST = 29
1768 IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1769 IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1770.LT. IF( IUIL ) THEN
1771 ITEMP = IU
1772 IU = IL
1773 IL = ITEMP
1774 END IF
1775 CALL ZSTEMR( 'v', 'i', N, D5, RWORK, VL, VU, IL, IU,
1776 $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
1777 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ),
1778 $ LIWORK-2*N, IINFO )
1779.NE. IF( IINFO0 ) THEN
1780 WRITE( NOUNIT, FMT = 9999 )'zstemr(v,i)', IINFO,
1781 $ N, JTYPE, IOLDSD
1782 INFO = ABS( IINFO )
1783.LT. IF( IINFO0 ) THEN
1784 RETURN
1785 ELSE
1786 RESULT( 29 ) = ULPINV
1787 GO TO 280
1788 END IF
1789 END IF
1790*
1791* Do Tests 29 and 30
1792*
1793* Call ZSTEMR to compute D2, do tests.
1794*
1795* Compute D2
1796*
1797 CALL DCOPY( N, SD, 1, D5, 1 )
1798.GT. IF( N0 )
1799 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1800*
1801 NTEST = 31
1802 CALL ZSTEMR( 'n', 'i', N, D5, RWORK, VL, VU, IL, IU,
1803 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
1804 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ),
1805 $ LIWORK-2*N, IINFO )
1806.NE. IF( IINFO0 ) THEN
1807 WRITE( NOUNIT, FMT = 9999 )'zstemr(n,i)', IINFO,
1808 $ N, JTYPE, IOLDSD
1809 INFO = ABS( IINFO )
1810.LT. IF( IINFO0 ) THEN
1811 RETURN
1812 ELSE
1813 RESULT( 31 ) = ULPINV
1814 GO TO 280
1815 END IF
1816 END IF
1817*
1818* Do Test 31
1819*
1820 TEMP1 = ZERO
1821 TEMP2 = ZERO
1822*
1823 DO 240 J = 1, IU - IL + 1
1824 TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
1825 $ ABS( D2( J ) ) )
1826 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1827 240 CONTINUE
1828*
1829 RESULT( 31 ) = TEMP2 / MAX( UNFL,
1830 $ ULP*MAX( TEMP1, TEMP2 ) )
1831*
1832* Call ZSTEMR(V,V) to compute D1 and Z, do tests.
1833*
1834* Compute D1 and Z
1835*
1836 CALL DCOPY( N, SD, 1, D5, 1 )
1837.GT. IF( N0 )
1838 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1839 CALL ZLASET( 'full', N, N, CZERO, CONE, Z, LDU )
1840*
1841 NTEST = 32
1842*
1843.GT. IF( N0 ) THEN
1844.NE. IF( IL1 ) THEN
1845 VL = D2( IL ) - MAX( HALF*
1846 $ ( D2( IL )-D2( IL-1 ) ), ULP*ANORM,
1847 $ TWO*RTUNFL )
1848 ELSE
1849 VL = D2( 1 ) - MAX( HALF*( D2( N )-D2( 1 ) ),
1850 $ ULP*ANORM, TWO*RTUNFL )
1851 END IF
1852.NE. IF( IUN ) THEN
1853 VU = D2( IU ) + MAX( HALF*
1854 $ ( D2( IU+1 )-D2( IU ) ), ULP*ANORM,
1855 $ TWO*RTUNFL )
1856 ELSE
1857 VU = D2( N ) + MAX( HALF*( D2( N )-D2( 1 ) ),
1858 $ ULP*ANORM, TWO*RTUNFL )
1859 END IF
1860 ELSE
1861 VL = ZERO
1862 VU = ONE
1863 END IF
1864*
1865 CALL ZSTEMR( 'v', 'v', N, D5, RWORK, VL, VU, IL, IU,
1866 $ M, D1, Z, LDU, M, IWORK( 1 ), TRYRAC,
1867 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ),
1868 $ LIWORK-2*N, IINFO )
1869.NE. IF( IINFO0 ) THEN
1870 WRITE( NOUNIT, FMT = 9999 )'zstemr(v,v)', IINFO,
1871 $ N, JTYPE, IOLDSD
1872 INFO = ABS( IINFO )
1873.LT. IF( IINFO0 ) THEN
1874 RETURN
1875 ELSE
1876 RESULT( 32 ) = ULPINV
1877 GO TO 280
1878 END IF
1879 END IF
1880*
1881* Do Tests 32 and 33
1882*
1883 CALL ZSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1884 $ M, RWORK, RESULT( 32 ) )
1885*
1886* Call ZSTEMR to compute D2, do tests.
1887*
1888* Compute D2
1889*
1890 CALL DCOPY( N, SD, 1, D5, 1 )
1891.GT. IF( N0 )
1892 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1893*
1894 NTEST = 34
1895 CALL ZSTEMR( 'n', 'v', N, D5, RWORK, VL, VU, IL, IU,
1896 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
1897 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ),
1898 $ LIWORK-2*N, IINFO )
1899.NE. IF( IINFO0 ) THEN
1900 WRITE( NOUNIT, FMT = 9999 )'zstemr(n,v)', IINFO,
1901 $ N, JTYPE, IOLDSD
1902 INFO = ABS( IINFO )
1903.LT. IF( IINFO0 ) THEN
1904 RETURN
1905 ELSE
1906 RESULT( 34 ) = ULPINV
1907 GO TO 280
1908 END IF
1909 END IF
1910*
1911* Do Test 34
1912*
1913 TEMP1 = ZERO
1914 TEMP2 = ZERO
1915*
1916 DO 250 J = 1, IU - IL + 1
1917 TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
1918 $ ABS( D2( J ) ) )
1919 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1920 250 CONTINUE
1921*
1922 RESULT( 34 ) = TEMP2 / MAX( UNFL,
1923 $ ULP*MAX( TEMP1, TEMP2 ) )
1924 ELSE
1925 RESULT( 29 ) = ZERO
1926 RESULT( 30 ) = ZERO
1927 RESULT( 31 ) = ZERO
1928 RESULT( 32 ) = ZERO
1929 RESULT( 33 ) = ZERO
1930 RESULT( 34 ) = ZERO
1931 END IF
1932*
1933* Call ZSTEMR(V,A) to compute D1 and Z, do tests.
1934*
1935* Compute D1 and Z
1936*
1937 CALL DCOPY( N, SD, 1, D5, 1 )
1938.GT. IF( N0 )
1939 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1940*
1941 NTEST = 35
1942*
1943 CALL ZSTEMR( 'v', 'a', N, D5, RWORK, VL, VU, IL, IU,
1944 $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
1945 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ),
1946 $ LIWORK-2*N, IINFO )
1947.NE. IF( IINFO0 ) THEN
1948 WRITE( NOUNIT, FMT = 9999 )'zstemr(v,a)', IINFO, N,
1949 $ JTYPE, IOLDSD
1950 INFO = ABS( IINFO )
1951.LT. IF( IINFO0 ) THEN
1952 RETURN
1953 ELSE
1954 RESULT( 35 ) = ULPINV
1955 GO TO 280
1956 END IF
1957 END IF
1958*
1959* Do Tests 35 and 36
1960*
1961 CALL ZSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M,
1962 $ RWORK, RESULT( 35 ) )
1963*
1964* Call ZSTEMR to compute D2, do tests.
1965*
1966* Compute D2
1967*
1968 CALL DCOPY( N, SD, 1, D5, 1 )
1969.GT. IF( N0 )
1970 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 )
1971*
1972 NTEST = 37
1973 CALL ZSTEMR( 'n', 'a', N, D5, RWORK, VL, VU, IL, IU,
1974 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
1975 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ),
1976 $ LIWORK-2*N, IINFO )
1977.NE. IF( IINFO0 ) THEN
1978 WRITE( NOUNIT, FMT = 9999 )'zstemr(n,a)', IINFO, N,
1979 $ JTYPE, IOLDSD
1980 INFO = ABS( IINFO )
1981.LT. IF( IINFO0 ) THEN
1982 RETURN
1983 ELSE
1984 RESULT( 37 ) = ULPINV
1985 GO TO 280
1986 END IF
1987 END IF
1988*
1989* Do Test 37
1990*
1991 TEMP1 = ZERO
1992 TEMP2 = ZERO
1993*
1994 DO 260 J = 1, N
1995 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1996 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1997 260 CONTINUE
1998*
1999 RESULT( 37 ) = TEMP2 / MAX( UNFL,
2000 $ ULP*MAX( TEMP1, TEMP2 ) )
2001 END IF
2002 270 CONTINUE
2003 280 CONTINUE
2004 NTESTT = NTESTT + NTEST
2005*
2006* End of Loop -- Check for RESULT(j) > THRESH
2007*
2008* Print out tests which fail.
2009*
2010 DO 290 JR = 1, NTEST
2011.GE. IF( RESULT( JR )THRESH ) THEN
2012*
2013* If this is the first test to fail,
2014* print a header to the data file.
2015*
2016.EQ. IF( NERRS0 ) THEN
2017 WRITE( NOUNIT, FMT = 9998 )'zst'
2018 WRITE( NOUNIT, FMT = 9997 )
2019 WRITE( NOUNIT, FMT = 9996 )
2020 WRITE( NOUNIT, FMT = 9995 )'hermitian'
2021 WRITE( NOUNIT, FMT = 9994 )
2022*
2023* Tests performed
2024*
2025 WRITE( NOUNIT, FMT = 9987 )
2026 END IF
2027 NERRS = NERRS + 1
2028.LT. IF( RESULT( JR )10000.0D0 ) THEN
2029 WRITE( NOUNIT, FMT = 9989 )N, JTYPE, IOLDSD, JR,
2030 $ RESULT( JR )
2031 ELSE
2032 WRITE( NOUNIT, FMT = 9988 )N, JTYPE, IOLDSD, JR,
2033 $ RESULT( JR )
2034 END IF
2035 END IF
2036 290 CONTINUE
2037 300 CONTINUE
2038 310 CONTINUE
2039*
2040* Summary
2041*
2042 CALL DLASUM( 'zst', NOUNIT, NERRS, NTESTT )
2043 RETURN
2044*
2045 9999 FORMAT( ' zchkst2stg: ', A, ' returned info=', I6, '.', / 9X,
2046 $ 'n=', I6, ', jtype=', I6, ', iseed=(', 3( I5, ',' ), I5, ')' )
2047*
2048 9998 FORMAT( / 1X, A3, ' -- Complex Hermitian eigenvalue problem' )
2049 9997 FORMAT( ' Matrix types (see ZCHKST2STG for details): ' )
2050*
2051 9996 FORMAT( / ' Special Matrices:',
2052 $ / ' 1=zero matrix. ',
2053 $ ' 5=diagonal: clustered entries.',
2054 $ / ' 2=identity matrix. ',
2055 $ ' 6=diagonal: large, evenly spaced.',
2056 $ / ' 3=diagonal: evenly spaced entries. ',
2057 $ ' 7=diagonal: small, evenly spaced.',
2058 $ / ' 4=diagonal: geometr. spaced entries.' )
2059 9995 FORMAT( ' dense ', A, ' matrices:',
2060 $ / ' 8=evenly spaced eigenvals. ',
2061 $ ' 12=small, evenly spaced eigenvals.',
2062 $ / ' 9=geometrically spaced eigenvals. ',
2063 $ ' 13=matrix with random o(1) entries.',
2064 $ / ' 10=clustered eigenvalues. ',
2065 $ ' 14=matrix with large random entries.',
2066 $ / ' 11=large, evenly spaced eigenvals. ',
2067 $ ' 15=matrix with small random entries.' )
2068 9994 FORMAT( ' 16=positive definite, evenly spaced eigenvalues',
2069 $ / ' 17=positive definite, geometrically spaced eigenvlaues',
2070 $ / ' 18=positive definite, clustered eigenvalues',
2071 $ / ' 19=positive definite, small evenly spaced eigenvalues',
2072 $ / ' 20=positive definite, large evenly spaced eigenvalues',
2073 $ / ' 21=diagonally dominant tridiagonal, geometrically',
2074 $ ' spaced eigenvalues' )
2075*
2076 9989 FORMAT( ' matrix order=', I5, ', type=', I2, ', seed=',
2077 $ 4( I4, ',' ), ' result ', I3, ' is', 0P, F8.2 )
2078 9988 FORMAT( ' matrix order=', I5, ', type=', I2, ', seed=',
2079 $ 4( I4, ',' ), ' result ', I3, ' is', 1P, D10.3 )
2080*
2081 9987 FORMAT( / 'test performed: see zchkst2stg for details.', / )
2082*
2083* End of ZCHKST2STG
2084*
2085 END
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dstebz
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
dsterf
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
zhetrd_2stage
subroutine zhetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
ZHETRD_2STAGE
Definition zhetrd_2stage.f:224
zhetrd
subroutine zhetrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
ZHETRD
Definition zhetrd.f:192
zlacpy
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
zlaset
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
zupgtr
subroutine zupgtr(uplo, n, ap, tau, q, ldq, work, info)
ZUPGTR
Definition zupgtr.f:114
zstein
subroutine zstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
ZSTEIN
Definition zstein.f:182
zsteqr
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
zstemr
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:338
zungtr
subroutine zungtr(uplo, n, a, lda, tau, work, lwork, info)
ZUNGTR
Definition zungtr.f:123
zstedc
subroutine zstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZSTEDC
Definition zstedc.f:212
zhptrd
subroutine zhptrd(uplo, n, ap, d, e, tau, info)
ZHPTRD
Definition zhptrd.f:151
zpteqr
subroutine zpteqr(compz, n, d, e, z, ldz, work, info)
ZPTEQR
Definition zpteqr.f:145
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zhet21
subroutine zhet21(itype, uplo, n, kband, a, lda, d, e, u, ldu, v, ldv, tau, work, rwork, result)
ZHET21
Definition zhet21.f:214
zstt22
subroutine zstt22(n, m, kband, ad, ae, sd, se, u, ldu, work, ldwork, rwork, result)
ZSTT22
Definition zstt22.f:145
zhpt21
subroutine zhpt21(itype, uplo, n, kband, ap, d, e, u, ldu, vp, tau, work, rwork, result)
ZHPT21
Definition zhpt21.f:228
zstt21
subroutine zstt21(n, kband, ad, ae, sd, se, u, ldu, work, rwork, result)
ZSTT21
Definition zstt21.f:133
zchkst2stg
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
zlatmr
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
zlatms
subroutine zlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
ZLATMS
Definition zlatms.f:332
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlasum
subroutine dlasum(type, iounit, ie, nrun)
DLASUM
Definition dlasum.f:43
dstech
subroutine dstech(n, a, b, eig, tol, work, info)
DSTECH
Definition dstech.f:101
seed
#define seed()
Definition macros.h:43
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21