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