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