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