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