OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
schkst2stg.f
Go to the documentation of this file.
1*> \brief \b SCHKST2STG
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 SCHKST2STG( 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* REAL THRESH
20* ..
21* .. Array Arguments ..
22* LOGICAL DOTYPE( * )
23* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
24* REAL 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*> SCHKST2STG 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*> SSYTRD. For that, we call the standard SSYTRD and compute D1 using
43*> SSTEQR, then we call the 2-stage SSYTRD_2STAGE with Upper and Lower
44*> and we compute D2 and D3 using SSTEQR 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 SCHKST in the next
48*> release when vectors and generation of Q will be implemented.
49*>
50*> SSYTRD factors A as U S U' , where ' means transpose,
51*> S is symmetric tridiagonal, and U is orthogonal.
52*> SSYTRD can use either just the lower or just the upper triangle
53*> of A; SCHKST2STG 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*> SSPTRD does the same as SSYTRD, except that A and V are stored
59*> in "packed" format.
60*>
61*> SORGTR constructs the matrix U from the contents of V and TAU.
62*>
63*> SOPGTR constructs the matrix U from the contents of VP and TAU.
64*>
65*> SSTEQR 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*> SSTERF computes D3, the matrix of eigenvalues, by the
71*> PWK method, which does not yield eigenvectors.
72*>
73*> SPTEQR 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*> SSTEBZ 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*> SSTEIN computes Y, the eigenvectors of S, given the
85*> eigenvalues.
86*>
87*> SSTEDC 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 SSYTRD/SORGTR or SSPTRD/SOPGTR ('V' option). It may
92*> also just compute eigenvalues ('N' option).
93*>
94*> SSTEMR 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). SSTEMR
97*> uses the Relatively Robust Representation whenever possible.
98*>
99*> When SCHKST2STG 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 ) SSYTRD( UPLO='U', ... )
106*>
107*> (2) | I - UV' | / ( n ulp ) SORGTR( UPLO='U', ... )
108*>
109*> (3) | A - V S V' | / ( |A| n ulp ) SSYTRD( 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*> SSYTRD_2STAGE("N", "U",....). D1 and D2 are computed
114*> via SSTEQR('N',...)
115*>
116*> (4) | I - UV' | / ( n ulp ) SORGTR( 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*> SSYTRD_2STAGE("N", "L",....). D1 and D3 are computed
121*> via SSTEQR('N',...)
122*>
123*> (5-8) Same as 1-4, but for SSPTRD and SOPGTR.
124*>
125*> (9) | S - Z D Z' | / ( |S| n ulp ) SSTEQR('V',...)
126*>
127*> (10) | I - ZZ' | / ( n ulp ) SSTEQR('V',...)
128*>
129*> (11) | D1 - D2 | / ( |D1| ulp ) SSTEQR('N',...)
130*>
131*> (12) | D1 - D3 | / ( |D1| ulp ) SSTERF
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*> SSTECH)
137*>
138*> For S positive definite,
139*>
140*> (14) | S - Z4 D4 Z4' | / ( |S| n ulp ) SPTEQR('V',...)
141*>
142*> (15) | I - Z4 Z4' | / ( n ulp ) SPTEQR('V',...)
143*>
144*> (16) | D4 - D5 | / ( 100 |D4| ulp ) SPTEQR('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*> SSTEBZ( 'A', 'E', ...)
152*>
153*> (18) | WA1 - D3 | / ( |D3| ulp ) SSTEBZ( '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*> SSTEBZ( 'I', 'E', ...)
160*>
161*> (20) | S - Y WA1 Y' | / ( |S| n ulp ) SSTEBZ, SSTEIN
162*>
163*> (21) | I - Y Y' | / ( n ulp ) SSTEBZ, SSTEIN
164*>
165*> (22) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('I')
166*>
167*> (23) | I - ZZ' | / ( n ulp ) SSTEDC('I')
168*>
169*> (24) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('V')
170*>
171*> (25) | I - ZZ' | / ( n ulp ) SSTEDC('V')
172*>
173*> (26) | D1 - D2 | / ( |D1| ulp ) SSTEDC('V') and
174*> SSTEDC('N')
175*>
176*> Test 27 is disabled at the moment because SSTEMR 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*> SSTEMR('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*> SSTEMR('V', 'I')
188*>
189*> Tests 29 through 34 are disable at present because SSTEMR
190*> does not handle partial spectrum requests.
191*>
192*> (29) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'I')
193*>
194*> (30) | I - ZZ' | / ( n ulp ) SSTEMR('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*> SSTEMR('N', 'I') vs. SSTEMR('V', 'I')
201*>
202*> (32) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'V')
203*>
204*> (33) | I - ZZ' | / ( n ulp ) SSTEMR('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*> SSTEMR('N', 'V') vs. SSTEMR('V', 'V')
211*>
212*> (35) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'A')
213*>
214*> (36) | I - ZZ' | / ( n ulp ) SSTEMR('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*> SSTEMR('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*> SCHKST2STG 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, SCHKST2STG
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 SCHKST2STG to continue the same random number
321*> sequence.
322*> \endverbatim
323*>
324*> \param[in] THRESH
325*> \verbatim
326*> THRESH is REAL
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 REAL 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 REAL 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 REAL array of
368*> dimension( max(NN) )
369*> The diagonal of the tridiagonal matrix computed by SSYTRD.
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 REAL array of
377*> dimension( max(NN) )
378*> The off-diagonal of the tridiagonal matrix computed by
379*> SSYTRD. 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 REAL array of
386*> dimension( max(NN) )
387*> The eigenvalues of A, as computed by SSTEQR 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 REAL array of
395*> dimension( max(NN) )
396*> The eigenvalues of A, as computed by SSTEQR 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 REAL array of
404*> dimension( max(NN) )
405*> The eigenvalues of A, as computed by SSTERF. 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 REAL array of
412*> dimension( max(NN) )
413*> The eigenvalues of A, as computed by SPTEQR(V).
414*> SPTEQR 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 REAL array of
421*> dimension( max(NN) )
422*> The eigenvalues of A, as computed by SPTEQR(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 REAL 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 SSTEBZ.
434*> \endverbatim
435*>
436*> \param[out] WA2
437*> \verbatim
438*> WA2 is REAL 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 SSTEBZ.
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 REAL 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 SSTEBZ.
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 REAL array of
461*> dimension( max(NN) )
462*> All eigenvalues of A, computed to high
463*> absolute accuracy, with different options.
464*> as computed by SSTEBZ.
465*> \endverbatim
466*>
467*> \param[out] U
468*> \verbatim
469*> U is REAL array of
470*> dimension( LDU, max(NN) ).
471*> The orthogonal matrix computed by SSYTRD + SORGTR.
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 REAL array of
484*> dimension( LDU, max(NN) ).
485*> The Housholder vectors computed by SSYTRD 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 SSYTRD, 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 SORGTR, 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 REAL 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 REAL array of
505*> dimension( max(NN) )
506*> The Householder factors computed by SSYTRD in reducing A
507*> to tridiagonal form.
508*> \endverbatim
509*>
510*> \param[out] Z
511*> \verbatim
512*> Z is REAL array of
513*> dimension( LDU, max(NN) ).
514*> The orthogonal matrix of eigenvectors computed by SSTEQR,
515*> SPTEQR, and SSTEIN.
516*> \endverbatim
517*>
518*> \param[out] WORK
519*> \verbatim
520*> WORK is REAL 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 REAL 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 SLATMR, SLATMS, SSYTRD, SORGTR, SSTEQR, SSTERF,
566*> or SORMC2 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 single_eig
606*
607* =====================================================================
608 SUBROUTINE schkst2stg( 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 REAL THRESH
621* ..
622* .. Array Arguments ..
623 LOGICAL DOTYPE( * )
624 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
625 REAL 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 REAL ZERO, ONE, TWO, EIGHT, TEN, HUN
636 PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0,
637 $ eight = 8.0e0, ten = 10.0e0, hun = 100.0e0 )
638 REAL 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 REAL 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 REAL DUMMA( 1 )
662* ..
663* .. External Functions ..
664 INTEGER ILAENV
665 REAL SLAMCH, SLARND, SSXT1
666 EXTERNAL ILAENV, SLAMCH, SLARND, SSXT1
667* ..
668* .. External Subroutines ..
669 EXTERNAL scopy, slabad, slacpy, slaset, slasum, slatmr,
670 $ slatms, sopgtr, sorgtr, spteqr, sspt21, ssptrd,
671 $ sstebz, sstech, sstedc, sstemr, sstein, ssteqr,
672 $ ssterf, sstt21, sstt22, ssyt21, ssytrd, xerbla,
673 $ ssytrd_2stage
674* ..
675* .. Intrinsic Functions ..
676 INTRINSIC abs, real, 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, 'SSYTRD', '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( 'SCHKST2STG', -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 = slamch( 'Safe minimum' )
739 ovfl = one / unfl
740 CALL slabad( unfl, ovfl )
741 ulp = slamch( 'Epsilon' )*slamch( '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( real( 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 / real( 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 slaset( '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 slatms( 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 slatms( 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 slatmr( 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 slatmr( 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 slatms( n, n, 'S', iseed, 'P', work, imode, cond,
891 $ anorm, n, n, 'N', a, lda, work( n+1 ),
892 $ iinfo )
893*
894 ELSE IF( itype.EQ.10 ) THEN
895*
896* Positive definite tridiagonal, eigenvalues specified.
897*
898 CALL slatms( 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 IF( temp1.GT.half ) 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 IF( iinfo.NE.0 ) 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 SSYTRD and SORGTR to compute S and U from
926* upper triangle.
927*
928 CALL slacpy( 'U', n, n, a, lda, v, ldu )
929*
930 ntest = 1
931 CALL ssytrd( 'U', n, v, ldu, sd, se, tau, work, lwork,
932 $ iinfo )
933*
934 IF( iinfo.NE.0 ) THEN
935 WRITE( nounit, fmt = 9999 )'SSYTRD(U)', iinfo, n, jtype,
936 $ ioldsd
937 info = abs( iinfo )
938 IF( iinfo.LT.0 ) THEN
939 RETURN
940 ELSE
941 result( 1 ) = ulpinv
942 GO TO 280
943 END IF
944 END IF
945*
946 CALL slacpy( 'U', n, n, v, ldu, u, ldu )
947*
948 ntest = 2
949 CALL sorgtr( 'U', n, u, ldu, tau, work, lwork, iinfo )
950 IF( iinfo.NE.0 ) THEN
951 WRITE( nounit, fmt = 9999 )'SORGTR(U)', iinfo, n, jtype,
952 $ ioldsd
953 info = abs( iinfo )
954 IF( iinfo.LT.0 ) 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 ssyt21( 2, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
965 $ ldu, tau, work, result( 1 ) )
966 CALL ssyt21( 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 scopy( n, sd, 1, d1, 1 )
977 IF( n.GT.0 )
978 $ CALL scopy( n-1, se, 1, work, 1 )
979*
980 CALL ssteqr( 'N', n, d1, work, work( n+1 ), ldu,
981 $ work( n+1 ), iinfo )
982 IF( iinfo.NE.0 ) THEN
983 WRITE( nounit, fmt = 9999 )'SSTEQR(N)', iinfo, n, jtype,
984 $ ioldsd
985 info = abs( iinfo )
986 IF( iinfo.LT.0 ) 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 slaset( 'Full', n, 1, zero, zero, sd, n )
1000 CALL slaset( 'Full', n, 1, zero, zero, se, n )
1001 CALL slacpy( "U", n, n, a, lda, v, ldu )
1002 lh = max(1, 4*n)
1003 lw = lwork - lh
1004 CALL ssytrd_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 scopy( n, sd, 1, d2, 1 )
1010 IF( n.GT.0 )
1011 $ CALL scopy( n-1, se, 1, work, 1 )
1012*
1013 CALL ssteqr( 'N', n, d2, work, work( n+1 ), ldu,
1014 $ work( n+1 ), iinfo )
1015 IF( iinfo.NE.0 ) THEN
1016 WRITE( nounit, fmt = 9999 )'SSTEQR(N)', iinfo, n, jtype,
1017 $ ioldsd
1018 info = abs( iinfo )
1019 IF( iinfo.LT.0 ) 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 slaset( 'Full', n, 1, zero, zero, sd, n )
1033 CALL slaset( 'Full', n, 1, zero, zero, se, n )
1034 CALL slacpy( "L", n, n, a, lda, v, ldu )
1035 CALL ssytrd_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 scopy( n, sd, 1, d3, 1 )
1041 IF( n.GT.0 )
1042 $ CALL scopy( n-1, se, 1, work, 1 )
1043*
1044 CALL ssteqr( 'N', n, d3, work, work( n+1 ), ldu,
1045 $ work( n+1 ), iinfo )
1046 IF( iinfo.NE.0 ) THEN
1047 WRITE( nounit, fmt = 9999 )'SSTEQR(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 SSPTRD and SOPGTR to compute S and U from AP
1088*
1089 CALL scopy( nap, ap, 1, vp, 1 )
1090*
1091 ntest = 5
1092 CALL ssptrd( 'U', n, vp, sd, se, tau, iinfo )
1093*
1094 IF( iinfo.NE.0 ) THEN
1095 WRITE( nounit, fmt = 9999 )'SSPTRD(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 sopgtr( 'U', n, vp, tau, u, ldu, work, iinfo )
1108 IF( iinfo.NE.0 ) THEN
1109 WRITE( nounit, fmt = 9999 )'SOPGTR(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 sspt21( 2, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1123 $ work, result( 5 ) )
1124 CALL sspt21( 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 SSPTRD and SOPGTR to compute S and U from AP
1138*
1139 CALL scopy( nap, ap, 1, vp, 1 )
1140*
1141 ntest = 7
1142 CALL ssptrd( 'L', n, vp, sd, se, tau, iinfo )
1143*
1144 IF( iinfo.NE.0 ) THEN
1145 WRITE( nounit, fmt = 9999 )'SSPTRD(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 sopgtr( 'L', n, vp, tau, u, ldu, work, iinfo )
1158 IF( iinfo.NE.0 ) THEN
1159 WRITE( nounit, fmt = 9999 )'SOPGTR(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 sspt21( 2, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1171 $ work, result( 7 ) )
1172 CALL sspt21( 3, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1173 $ work, result( 8 ) )
1174*
1175* Call SSTEQR to compute D1, D2, and Z, do tests.
1176*
1177* Compute D1 and Z
1178*
1179 CALL scopy( n, sd, 1, d1, 1 )
1180 IF( n.GT.0 )
1181 $ CALL scopy( n-1, se, 1, work, 1 )
1182 CALL slaset( 'Full', n, n, zero, one, z, ldu )
1183*
1184 ntest = 9
1185 CALL ssteqr( 'V', n, d1, work, z, ldu, work( n+1 ), iinfo )
1186 IF( iinfo.NE.0 ) THEN
1187 WRITE( nounit, fmt = 9999 )'SSTEQR(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 scopy( n, sd, 1, d2, 1 )
1201 IF( n.GT.0 )
1202 $ CALL scopy( n-1, se, 1, work, 1 )
1203*
1204 ntest = 11
1205 CALL ssteqr( 'N', n, d2, work, work( n+1 ), ldu,
1206 $ work( n+1 ), iinfo )
1207 IF( iinfo.NE.0 ) THEN
1208 WRITE( nounit, fmt = 9999 )'SSTEQR(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 scopy( n, sd, 1, d3, 1 )
1222 IF( n.GT.0 )
1223 $ CALL scopy( n-1, se, 1, work, 1 )
1224*
1225 ntest = 12
1226 CALL ssterf( n, d3, work, iinfo )
1227 IF( iinfo.NE.0 ) THEN
1228 WRITE( nounit, fmt = 9999 )'SSTERF', 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 sstt21( 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 sstech( 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 SPTEQR
1278* and do tests 14, 15, and 16 .
1279*
1280 IF( jtype.GT.15 ) THEN
1281*
1282* Compute D4 and Z4
1283*
1284 CALL scopy( n, sd, 1, d4, 1 )
1285 IF( n.GT.0 )
1286 $ CALL scopy( n-1, se, 1, work, 1 )
1287 CALL slaset( 'Full', n, n, zero, one, z, ldu )
1288*
1289 ntest = 14
1290 CALL spteqr( 'V', n, d4, work, z, ldu, work( n+1 ),
1291 $ iinfo )
1292 IF( iinfo.NE.0 ) THEN
1293 WRITE( nounit, fmt = 9999 )'spteqr(v)', IINFO, N,
1294 $ JTYPE, IOLDSD
1295 INFO = ABS( IINFO )
1296.LT. IF( IINFO0 ) 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 SSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK,
1307 $ RESULT( 14 ) )
1308*
1309* Compute D5
1310*
1311 CALL SCOPY( N, SD, 1, D5, 1 )
1312.GT. IF( N0 )
1313 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1314*
1315 NTEST = 16
1316 CALL SPTEQR( 'n', N, D5, WORK, Z, LDU, WORK( N+1 ),
1317 $ IINFO )
1318.NE. IF( IINFO0 ) THEN
1319 WRITE( NOUNIT, FMT = 9999 )'spteqr(n)', IINFO, N,
1320 $ JTYPE, IOLDSD
1321 INFO = ABS( IINFO )
1322.LT. IF( IINFO0 ) 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 SSTEBZ 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.EQ. IF( JTYPE21 ) THEN
1357 NTEST = 17
1358 ABSTOL = UNFL + UNFL
1359 CALL SSTEBZ( '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.NE. IF( IINFO0 ) THEN
1363 WRITE( NOUNIT, FMT = 9999 )'sstebz(a,rel)', IINFO, N,
1364 $ JTYPE, IOLDSD
1365 INFO = ABS( IINFO )
1366.LT. IF( IINFO0 ) 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 SSTEBZ( '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.NE. IF( IINFO0 ) THEN
1398 WRITE( NOUNIT, FMT = 9999 )'sstebz(a)', IINFO, N, JTYPE,
1399 $ IOLDSD
1400 INFO = ABS( IINFO )
1401.LT. IF( IINFO0 ) 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.LE. IF( N1 ) THEN
1425 IL = 1
1426 IU = N
1427 ELSE
1428 IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
1429 IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
1430.LT. IF( IUIL ) THEN
1431 ITEMP = IU
1432 IU = IL
1433 IL = ITEMP
1434 END IF
1435 END IF
1436*
1437 CALL SSTEBZ( '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.NE. IF( IINFO0 ) THEN
1441 WRITE( NOUNIT, FMT = 9999 )'sstebz(i)', IINFO, N, JTYPE,
1442 $ IOLDSD
1443 INFO = ABS( IINFO )
1444.LT. IF( IINFO0 ) 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.GT. IF( N0 ) THEN
1456.NE. IF( IL1 ) 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.NE. IF( IUN ) 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 SSTEBZ( '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.NE. IF( IINFO0 ) THEN
1479 WRITE( NOUNIT, FMT = 9999 )'sstebz(v)', IINFO, N, JTYPE,
1480 $ IOLDSD
1481 INFO = ABS( IINFO )
1482.LT. IF( IINFO0 ) THEN
1483 RETURN
1484 ELSE
1485 RESULT( 19 ) = ULPINV
1486 GO TO 280
1487 END IF
1488 END IF
1489*
1490.EQ..AND..NE. IF( M30 N0 ) THEN
1491 RESULT( 19 ) = ULPINV
1492 GO TO 280
1493 END IF
1494*
1495* Do test 19
1496*
1497 TEMP1 = SSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
1498 TEMP2 = SSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL )
1499.GT. IF( N0 ) 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 SSTEIN to compute eigenvectors corresponding to
1508* eigenvalues in WA1. (First call SSTEBZ again, to make sure
1509* it returns these eigenvalues in the correct order.)
1510*
1511 NTEST = 21
1512 CALL SSTEBZ( '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.NE. IF( IINFO0 ) THEN
1516 WRITE( NOUNIT, FMT = 9999 )'sstebz(a,b)', IINFO, N,
1517 $ JTYPE, IOLDSD
1518 INFO = ABS( IINFO )
1519.LT. IF( IINFO0 ) 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 SSTEIN( 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.NE. IF( IINFO0 ) THEN
1532 WRITE( NOUNIT, FMT = 9999 )'sstein', IINFO, N, JTYPE,
1533 $ IOLDSD
1534 INFO = ABS( IINFO )
1535.LT. IF( IINFO0 ) 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 SSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK,
1547 $ RESULT( 20 ) )
1548*
1549* Call SSTEDC(I) to compute D1 and Z, do tests.
1550*
1551* Compute D1 and Z
1552*
1553 CALL SCOPY( N, SD, 1, D1, 1 )
1554.GT. IF( N0 )
1555 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1556 CALL SLASET( 'full', N, N, ZERO, ONE, Z, LDU )
1557*
1558 NTEST = 22
1559 CALL SSTEDC( 'i', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
1560 $ IWORK, LIWEDC, IINFO )
1561.NE. IF( IINFO0 ) THEN
1562 WRITE( NOUNIT, FMT = 9999 )'sstedc(i)', IINFO, N, JTYPE,
1563 $ IOLDSD
1564 INFO = ABS( IINFO )
1565.LT. IF( IINFO0 ) 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 SSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1576 $ RESULT( 22 ) )
1577*
1578* Call SSTEDC(V) to compute D1 and Z, do tests.
1579*
1580* Compute D1 and Z
1581*
1582 CALL SCOPY( N, SD, 1, D1, 1 )
1583.GT. IF( N0 )
1584 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1585 CALL SLASET( 'full', N, N, ZERO, ONE, Z, LDU )
1586*
1587 NTEST = 24
1588 CALL SSTEDC( 'v', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
1589 $ IWORK, LIWEDC, IINFO )
1590.NE. IF( IINFO0 ) THEN
1591 WRITE( NOUNIT, FMT = 9999 )'sstedc(v)', IINFO, N, JTYPE,
1592 $ IOLDSD
1593 INFO = ABS( IINFO )
1594.LT. IF( IINFO0 ) 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 SSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1605 $ RESULT( 24 ) )
1606*
1607* Call SSTEDC(N) to compute D2, do tests.
1608*
1609* Compute D2
1610*
1611 CALL SCOPY( N, SD, 1, D2, 1 )
1612.GT. IF( N0 )
1613 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1614 CALL SLASET( 'full', N, N, ZERO, ONE, Z, LDU )
1615*
1616 NTEST = 26
1617 CALL SSTEDC( '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 )'sstedc(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 SSTEMR if IEEE compliant
1644*
1645 IF( ILAENV( 10, 'sstemr', 'va.EQ..AND.', 1, 0, 0, 0 )1
1646 $ ILAENV( 11, 'sstemr', 'va.EQ.', 1, 0, 0, 0 )1 ) THEN
1647*
1648* Call SSTEMR, 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 SSTEMR( '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 )'sstemr(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( SLARND( 1, ISEED2 ) )
1690 IU = 1 + ( N-1 )*INT( SLARND( 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 SSTEMR( '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 )'sstemr(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 SSTEMR(V,I) to compute D1 and Z, do tests.
1738*
1739* Compute D1 and Z
1740*
1741 CALL SCOPY( N, SD, 1, D5, 1 )
1742.GT. IF( N0 )
1743 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1744 CALL SLASET( 'full', N, N, ZERO, ONE, Z, LDU )
1745*
1746 IF( SRANGE ) THEN
1747 NTEST = 29
1748 IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
1749 IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
1750.LT. IF( IUIL ) THEN
1751 ITEMP = IU
1752 IU = IL
1753 IL = ITEMP
1754 END IF
1755 CALL SSTEMR( '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 )'sstemr(v,i)', IINFO,
1761 $ N, JTYPE, IOLDSD
1762 INFO = ABS( IINFO )
1763.LT. IF( IINFO0 ) 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 SSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1774 $ M, RESULT( 29 ) )
1775*
1776* Call SSTEMR to compute D2, do tests.
1777*
1778* Compute D2
1779*
1780 CALL SCOPY( N, SD, 1, D5, 1 )
1781.GT. IF( N0 )
1782 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1783*
1784 NTEST = 31
1785 CALL SSTEMR( '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.NE. IF( IINFO0 ) THEN
1790 WRITE( NOUNIT, FMT = 9999 )'sstemr(n,i)', IINFO,
1791 $ N, JTYPE, IOLDSD
1792 INFO = ABS( IINFO )
1793.LT. IF( IINFO0 ) 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 SSTEMR(V,V) to compute D1 and Z, do tests.
1816*
1817* Compute D1 and Z
1818*
1819 CALL SCOPY( N, SD, 1, D5, 1 )
1820.GT. IF( N0 )
1821 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1822 CALL SLASET( 'full', N, N, ZERO, ONE, Z, LDU )
1823*
1824 NTEST = 32
1825*
1826.GT. IF( N0 ) THEN
1827.NE. IF( IL1 ) 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.NE. IF( IUN ) 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 SSTEMR( '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.NE. IF( IINFO0 ) THEN
1853 WRITE( NOUNIT, FMT = 9999 )'sstemr(v,v)', IINFO,
1854 $ N, JTYPE, IOLDSD
1855 INFO = ABS( IINFO )
1856.LT. IF( IINFO0 ) 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 SSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1867 $ M, RESULT( 32 ) )
1868*
1869* Call SSTEMR to compute D2, do tests.
1870*
1871* Compute D2
1872*
1873 CALL SCOPY( N, SD, 1, D5, 1 )
1874.GT. IF( N0 )
1875 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1876*
1877 NTEST = 34
1878 CALL SSTEMR( '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.NE. IF( IINFO0 ) THEN
1883 WRITE( NOUNIT, FMT = 9999 )'sstemr(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 SSTEMR(V,A) to compute D1 and Z, do tests.
1917*
1918* Compute D1 and Z
1919*
1920 CALL SCOPY( N, SD, 1, D5, 1 )
1921.GT. IF( N0 )
1922 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1923*
1924 NTEST = 35
1925*
1926 CALL SSTEMR( '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 )'sstemr(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 SSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M,
1945 $ RESULT( 35 ) )
1946*
1947* Call SSTEMR to compute D2, do tests.
1948*
1949* Compute D2
1950*
1951 CALL SCOPY( N, SD, 1, D5, 1 )
1952.GT. IF( N0 )
1953 $ CALL SCOPY( N-1, SE, 1, WORK, 1 )
1954*
1955 NTEST = 37
1956 CALL SSTEMR( '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 )'sstemr(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 )'sst'
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 SLASUM( 'sst', NOUNIT, NERRS, NTESTT )
2021 RETURN
2022*
2023 9999 FORMAT( ' schkst2stg: ', 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 schkst2stg 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 schkst2stg for details.', / )
2058*
2059* End of SCHKST2STG
2060*
2061 END
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
ssteqr
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:131
sstedc
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:188
sstebz
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:273
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
sstemr
subroutine sstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
SSTEMR
Definition sstemr.f:321
sopgtr
subroutine sopgtr(uplo, n, ap, tau, q, ldq, work, info)
SOPGTR
Definition sopgtr.f:114
sorgtr
subroutine sorgtr(uplo, n, a, lda, tau, work, lwork, info)
SORGTR
Definition sorgtr.f:123
ssptrd
subroutine ssptrd(uplo, n, ap, d, e, tau, info)
SSPTRD
Definition ssptrd.f:150
sstein
subroutine sstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
SSTEIN
Definition sstein.f:174
spteqr
subroutine spteqr(compz, n, d, e, z, ldz, work, info)
SPTEQR
Definition spteqr.f:145
ssytrd_2stage
subroutine ssytrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
SSYTRD_2STAGE
Definition ssytrd_2stage.f:224
ssytrd
subroutine ssytrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
SSYTRD
Definition ssytrd.f:192
slatmr
subroutine slatmr(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)
SLATMR
Definition slatmr.f:471
slatms
subroutine slatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
SLATMS
Definition slatms.f:321
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sspt21
subroutine sspt21(itype, uplo, n, kband, ap, d, e, u, ldu, vp, tau, work, result)
SSPT21
Definition sspt21.f:221
sstt21
subroutine sstt21(n, kband, ad, ae, sd, se, u, ldu, work, result)
SSTT21
Definition sstt21.f:127
schkst2stg
subroutine schkst2stg(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)
SCHKST2STG
Definition schkst2stg.f:612
sstt22
subroutine sstt22(n, m, kband, ad, ae, sd, se, u, ldu, work, ldwork, result)
SSTT22
Definition sstt22.f:139
sstech
subroutine sstech(n, a, b, eig, tol, work, info)
SSTECH
Definition sstech.f:101
ssyt21
subroutine ssyt21(itype, uplo, n, kband, a, lda, d, e, u, ldu, v, ldv, tau, work, result)
SSYT21
Definition ssyt21.f:207
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
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339
slasum
subroutine slasum(type, iounit, ie, nrun)
SLASUM
Definition slasum.f:41