LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ddrvst.f
Go to the documentation of this file.
1 *> \brief \b DDRVST
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 DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * NOUNIT, A, LDA, D1, D2, D3, D4, EVEIGS, WA1,
13 * WA2, WA3, U, LDU, V, TAU, Z, WORK, LWORK,
14 * IWORK, LIWORK, RESULT, INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
18 * $ NTYPES
19 * DOUBLE PRECISION THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
24 * DOUBLE PRECISION A( LDA, * ), D1( * ), D2( * ), D3( * ),
25 * $ D4( * ), EVEIGS( * ), RESULT( * ), TAU( * ),
26 * $ U( LDU, * ), V( LDU, * ), WA1( * ), WA2( * ),
27 * $ WA3( * ), WORK( * ), Z( LDU, * )
28 * ..
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> DDRVST checks the symmetric eigenvalue problem drivers.
37 *>
38 *> DSTEV computes all eigenvalues and, optionally,
39 *> eigenvectors of a real symmetric tridiagonal matrix.
40 *>
41 *> DSTEVX computes selected eigenvalues and, optionally,
42 *> eigenvectors of a real symmetric tridiagonal matrix.
43 *>
44 *> DSTEVR computes selected eigenvalues and, optionally,
45 *> eigenvectors of a real symmetric tridiagonal matrix
46 *> using the Relatively Robust Representation where it can.
47 *>
48 *> DSYEV computes all eigenvalues and, optionally,
49 *> eigenvectors of a real symmetric matrix.
50 *>
51 *> DSYEVX computes selected eigenvalues and, optionally,
52 *> eigenvectors of a real symmetric matrix.
53 *>
54 *> DSYEVR computes selected eigenvalues and, optionally,
55 *> eigenvectors of a real symmetric matrix
56 *> using the Relatively Robust Representation where it can.
57 *>
58 *> DSPEV computes all eigenvalues and, optionally,
59 *> eigenvectors of a real symmetric matrix in packed
60 *> storage.
61 *>
62 *> DSPEVX computes selected eigenvalues and, optionally,
63 *> eigenvectors of a real symmetric matrix in packed
64 *> storage.
65 *>
66 *> DSBEV computes all eigenvalues and, optionally,
67 *> eigenvectors of a real symmetric band matrix.
68 *>
69 *> DSBEVX computes selected eigenvalues and, optionally,
70 *> eigenvectors of a real symmetric band matrix.
71 *>
72 *> DSYEVD computes all eigenvalues and, optionally,
73 *> eigenvectors of a real symmetric matrix using
74 *> a divide and conquer algorithm.
75 *>
76 *> DSPEVD computes all eigenvalues and, optionally,
77 *> eigenvectors of a real symmetric matrix in packed
78 *> storage, using a divide and conquer algorithm.
79 *>
80 *> DSBEVD computes all eigenvalues and, optionally,
81 *> eigenvectors of a real symmetric band matrix,
82 *> using a divide and conquer algorithm.
83 *>
84 *> When DDRVST is called, a number of matrix "sizes" ("n's") and a
85 *> number of matrix "types" are specified. For each size ("n")
86 *> and each type of matrix, one matrix will be generated and used
87 *> to test the appropriate drivers. For each matrix and each
88 *> driver routine called, the following tests will be performed:
89 *>
90 *> (1) | A - Z D Z' | / ( |A| n ulp )
91 *>
92 *> (2) | I - Z Z' | / ( n ulp )
93 *>
94 *> (3) | D1 - D2 | / ( |D1| ulp )
95 *>
96 *> where Z is the matrix of eigenvectors returned when the
97 *> eigenvector option is given and D1 and D2 are the eigenvalues
98 *> returned with and without the eigenvector option.
99 *>
100 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
101 *> each element NN(j) specifies one size.
102 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
103 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
104 *> Currently, the list of possible types is:
105 *>
106 *> (1) The zero matrix.
107 *> (2) The identity matrix.
108 *>
109 *> (3) A diagonal matrix with evenly spaced eigenvalues
110 *> 1, ..., ULP and random signs.
111 *> (ULP = (first number larger than 1) - 1 )
112 *> (4) A diagonal matrix with geometrically spaced eigenvalues
113 *> 1, ..., ULP and random signs.
114 *> (5) A diagonal matrix with "clustered" eigenvalues
115 *> 1, ULP, ..., ULP and random signs.
116 *>
117 *> (6) Same as (4), but multiplied by SQRT( overflow threshold )
118 *> (7) Same as (4), but multiplied by SQRT( underflow threshold )
119 *>
120 *> (8) A matrix of the form U' D U, where U is orthogonal and
121 *> D has evenly spaced entries 1, ..., ULP with random signs
122 *> on the diagonal.
123 *>
124 *> (9) A matrix of the form U' D U, where U is orthogonal and
125 *> D has geometrically spaced entries 1, ..., ULP with random
126 *> signs on the diagonal.
127 *>
128 *> (10) A matrix of the form U' D U, where U is orthogonal and
129 *> D has "clustered" entries 1, ULP,..., ULP with random
130 *> signs on the diagonal.
131 *>
132 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
133 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
134 *>
135 *> (13) Symmetric matrix with random entries chosen from (-1,1).
136 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
137 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
138 *> (16) A band matrix with half bandwidth randomly chosen between
139 *> 0 and N-1, with evenly spaced eigenvalues 1, ..., ULP
140 *> with random signs.
141 *> (17) Same as (16), but multiplied by SQRT( overflow threshold )
142 *> (18) Same as (16), but multiplied by SQRT( underflow threshold )
143 *> \endverbatim
144 *
145 * Arguments:
146 * ==========
147 *
148 *> \verbatim
149 *> NSIZES INTEGER
150 *> The number of sizes of matrices to use. If it is zero,
151 *> DDRVST does nothing. It must be at least zero.
152 *> Not modified.
153 *>
154 *> NN INTEGER array, dimension (NSIZES)
155 *> An array containing the sizes to be used for the matrices.
156 *> Zero values will be skipped. The values must be at least
157 *> zero.
158 *> Not modified.
159 *>
160 *> NTYPES INTEGER
161 *> The number of elements in DOTYPE. If it is zero, DDRVST
162 *> does nothing. It must be at least zero. If it is MAXTYP+1
163 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
164 *> defined, which is to use whatever matrix is in A. This
165 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
166 *> DOTYPE(MAXTYP+1) is .TRUE. .
167 *> Not modified.
168 *>
169 *> DOTYPE LOGICAL array, dimension (NTYPES)
170 *> If DOTYPE(j) is .TRUE., then for each size in NN a
171 *> matrix of that size and of type j will be generated.
172 *> If NTYPES is smaller than the maximum number of types
173 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
174 *> MAXTYP will not be generated. If NTYPES is larger
175 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
176 *> will be ignored.
177 *> Not modified.
178 *>
179 *> ISEED INTEGER array, dimension (4)
180 *> On entry ISEED specifies the seed of the random number
181 *> generator. The array elements should be between 0 and 4095;
182 *> if not they will be reduced mod 4096. Also, ISEED(4) must
183 *> be odd. The random number generator uses a linear
184 *> congruential sequence limited to small integers, and so
185 *> should produce machine independent random numbers. The
186 *> values of ISEED are changed on exit, and can be used in the
187 *> next call to DDRVST to continue the same random number
188 *> sequence.
189 *> Modified.
190 *>
191 *> THRESH DOUBLE PRECISION
192 *> A test will count as "failed" if the "error", computed as
193 *> described above, exceeds THRESH. Note that the error
194 *> is scaled to be O(1), so THRESH should be a reasonably
195 *> small multiple of 1, e.g., 10 or 100. In particular,
196 *> it should not depend on the precision (single vs. double)
197 *> or the size of the matrix. It must be at least zero.
198 *> Not modified.
199 *>
200 *> NOUNIT INTEGER
201 *> The FORTRAN unit number for printing out error messages
202 *> (e.g., if a routine returns IINFO not equal to 0.)
203 *> Not modified.
204 *>
205 *> A DOUBLE PRECISION array, dimension (LDA , max(NN))
206 *> Used to hold the matrix whose eigenvalues are to be
207 *> computed. On exit, A contains the last matrix actually
208 *> used.
209 *> Modified.
210 *>
211 *> LDA INTEGER
212 *> The leading dimension of A. It must be at
213 *> least 1 and at least max( NN ).
214 *> Not modified.
215 *>
216 *> D1 DOUBLE PRECISION array, dimension (max(NN))
217 *> The eigenvalues of A, as computed by DSTEQR simlutaneously
218 *> with Z. On exit, the eigenvalues in D1 correspond with the
219 *> matrix in A.
220 *> Modified.
221 *>
222 *> D2 DOUBLE PRECISION array, dimension (max(NN))
223 *> The eigenvalues of A, as computed by DSTEQR if Z is not
224 *> computed. On exit, the eigenvalues in D2 correspond with
225 *> the matrix in A.
226 *> Modified.
227 *>
228 *> D3 DOUBLE PRECISION array, dimension (max(NN))
229 *> The eigenvalues of A, as computed by DSTERF. On exit, the
230 *> eigenvalues in D3 correspond with the matrix in A.
231 *> Modified.
232 *>
233 *> D4 DOUBLE PRECISION array, dimension
234 *>
235 *> EVEIGS DOUBLE PRECISION array, dimension (max(NN))
236 *> The eigenvalues as computed by DSTEV('N', ... )
237 *> (I reserve the right to change this to the output of
238 *> whichever algorithm computes the most accurate eigenvalues).
239 *>
240 *> WA1 DOUBLE PRECISION array, dimension
241 *>
242 *> WA2 DOUBLE PRECISION array, dimension
243 *>
244 *> WA3 DOUBLE PRECISION array, dimension
245 *>
246 *> U DOUBLE PRECISION array, dimension (LDU, max(NN))
247 *> The orthogonal matrix computed by DSYTRD + DORGTR.
248 *> Modified.
249 *>
250 *> LDU INTEGER
251 *> The leading dimension of U, Z, and V. It must be at
252 *> least 1 and at least max( NN ).
253 *> Not modified.
254 *>
255 *> V DOUBLE PRECISION array, dimension (LDU, max(NN))
256 *> The Housholder vectors computed by DSYTRD in reducing A to
257 *> tridiagonal form.
258 *> Modified.
259 *>
260 *> TAU DOUBLE PRECISION array, dimension (max(NN))
261 *> The Householder factors computed by DSYTRD in reducing A
262 *> to tridiagonal form.
263 *> Modified.
264 *>
265 *> Z DOUBLE PRECISION array, dimension (LDU, max(NN))
266 *> The orthogonal matrix of eigenvectors computed by DSTEQR,
267 *> DPTEQR, and DSTEIN.
268 *> Modified.
269 *>
270 *> WORK DOUBLE PRECISION array, dimension (LWORK)
271 *> Workspace.
272 *> Modified.
273 *>
274 *> LWORK INTEGER
275 *> The number of entries in WORK. This must be at least
276 *> 1 + 4 * Nmax + 2 * Nmax * lg Nmax + 4 * Nmax**2
277 *> where Nmax = max( NN(j), 2 ) and lg = log base 2.
278 *> Not modified.
279 *>
280 *> IWORK INTEGER array,
281 *> dimension (6 + 6*Nmax + 5 * Nmax * lg Nmax )
282 *> where Nmax = max( NN(j), 2 ) and lg = log base 2.
283 *> Workspace.
284 *> Modified.
285 *>
286 *> RESULT DOUBLE PRECISION array, dimension (105)
287 *> The values computed by the tests described above.
288 *> The values are currently limited to 1/ulp, to avoid
289 *> overflow.
290 *> Modified.
291 *>
292 *> INFO INTEGER
293 *> If 0, then everything ran OK.
294 *> -1: NSIZES < 0
295 *> -2: Some NN(j) < 0
296 *> -3: NTYPES < 0
297 *> -5: THRESH < 0
298 *> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
299 *> -16: LDU < 1 or LDU < NMAX.
300 *> -21: LWORK too small.
301 *> If DLATMR, DLATMS, DSYTRD, DORGTR, DSTEQR, DSTERF,
302 *> or DORMTR returns an error code, the
303 *> absolute value of it is returned.
304 *> Modified.
305 *>
306 *>-----------------------------------------------------------------------
307 *>
308 *> Some Local Variables and Parameters:
309 *> ---- ----- --------- --- ----------
310 *> ZERO, ONE Real 0 and 1.
311 *> MAXTYP The number of types defined.
312 *> NTEST The number of tests performed, or which can
313 *> be performed so far, for the current matrix.
314 *> NTESTT The total number of tests performed so far.
315 *> NMAX Largest value in NN.
316 *> NMATS The number of matrices generated so far.
317 *> NERRS The number of tests which have exceeded THRESH
318 *> so far (computed by DLAFTS).
319 *> COND, IMODE Values to be passed to the matrix generators.
320 *> ANORM Norm of A; passed to matrix generators.
321 *>
322 *> OVFL, UNFL Overflow and underflow thresholds.
323 *> ULP, ULPINV Finest relative precision and its inverse.
324 *> RTOVFL, RTUNFL Square roots of the previous 2 values.
325 *> The following four arrays decode JTYPE:
326 *> KTYPE(j) The general type (1-10) for type "j".
327 *> KMODE(j) The MODE value to be passed to the matrix
328 *> generator for type "j".
329 *> KMAGN(j) The order of magnitude ( O(1),
330 *> O(overflow^(1/2) ), O(underflow^(1/2) )
331 *>
332 *> The tests performed are: Routine tested
333 *> 1= | A - U S U' | / ( |A| n ulp ) DSTEV('V', ... )
334 *> 2= | I - U U' | / ( n ulp ) DSTEV('V', ... )
335 *> 3= |D(with Z) - D(w/o Z)| / (|D| ulp) DSTEV('N', ... )
336 *> 4= | A - U S U' | / ( |A| n ulp ) DSTEVX('V','A', ... )
337 *> 5= | I - U U' | / ( n ulp ) DSTEVX('V','A', ... )
338 *> 6= |D(with Z) - EVEIGS| / (|D| ulp) DSTEVX('N','A', ... )
339 *> 7= | A - U S U' | / ( |A| n ulp ) DSTEVR('V','A', ... )
340 *> 8= | I - U U' | / ( n ulp ) DSTEVR('V','A', ... )
341 *> 9= |D(with Z) - EVEIGS| / (|D| ulp) DSTEVR('N','A', ... )
342 *> 10= | A - U S U' | / ( |A| n ulp ) DSTEVX('V','I', ... )
343 *> 11= | I - U U' | / ( n ulp ) DSTEVX('V','I', ... )
344 *> 12= |D(with Z) - D(w/o Z)| / (|D| ulp) DSTEVX('N','I', ... )
345 *> 13= | A - U S U' | / ( |A| n ulp ) DSTEVX('V','V', ... )
346 *> 14= | I - U U' | / ( n ulp ) DSTEVX('V','V', ... )
347 *> 15= |D(with Z) - D(w/o Z)| / (|D| ulp) DSTEVX('N','V', ... )
348 *> 16= | A - U S U' | / ( |A| n ulp ) DSTEVD('V', ... )
349 *> 17= | I - U U' | / ( n ulp ) DSTEVD('V', ... )
350 *> 18= |D(with Z) - EVEIGS| / (|D| ulp) DSTEVD('N', ... )
351 *> 19= | A - U S U' | / ( |A| n ulp ) DSTEVR('V','I', ... )
352 *> 20= | I - U U' | / ( n ulp ) DSTEVR('V','I', ... )
353 *> 21= |D(with Z) - D(w/o Z)| / (|D| ulp) DSTEVR('N','I', ... )
354 *> 22= | A - U S U' | / ( |A| n ulp ) DSTEVR('V','V', ... )
355 *> 23= | I - U U' | / ( n ulp ) DSTEVR('V','V', ... )
356 *> 24= |D(with Z) - D(w/o Z)| / (|D| ulp) DSTEVR('N','V', ... )
357 *>
358 *> 25= | A - U S U' | / ( |A| n ulp ) DSYEV('L','V', ... )
359 *> 26= | I - U U' | / ( n ulp ) DSYEV('L','V', ... )
360 *> 27= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEV('L','N', ... )
361 *> 28= | A - U S U' | / ( |A| n ulp ) DSYEVX('L','V','A', ... )
362 *> 29= | I - U U' | / ( n ulp ) DSYEVX('L','V','A', ... )
363 *> 30= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVX('L','N','A', ... )
364 *> 31= | A - U S U' | / ( |A| n ulp ) DSYEVX('L','V','I', ... )
365 *> 32= | I - U U' | / ( n ulp ) DSYEVX('L','V','I', ... )
366 *> 33= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVX('L','N','I', ... )
367 *> 34= | A - U S U' | / ( |A| n ulp ) DSYEVX('L','V','V', ... )
368 *> 35= | I - U U' | / ( n ulp ) DSYEVX('L','V','V', ... )
369 *> 36= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVX('L','N','V', ... )
370 *> 37= | A - U S U' | / ( |A| n ulp ) DSPEV('L','V', ... )
371 *> 38= | I - U U' | / ( n ulp ) DSPEV('L','V', ... )
372 *> 39= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEV('L','N', ... )
373 *> 40= | A - U S U' | / ( |A| n ulp ) DSPEVX('L','V','A', ... )
374 *> 41= | I - U U' | / ( n ulp ) DSPEVX('L','V','A', ... )
375 *> 42= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVX('L','N','A', ... )
376 *> 43= | A - U S U' | / ( |A| n ulp ) DSPEVX('L','V','I', ... )
377 *> 44= | I - U U' | / ( n ulp ) DSPEVX('L','V','I', ... )
378 *> 45= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVX('L','N','I', ... )
379 *> 46= | A - U S U' | / ( |A| n ulp ) DSPEVX('L','V','V', ... )
380 *> 47= | I - U U' | / ( n ulp ) DSPEVX('L','V','V', ... )
381 *> 48= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVX('L','N','V', ... )
382 *> 49= | A - U S U' | / ( |A| n ulp ) DSBEV('L','V', ... )
383 *> 50= | I - U U' | / ( n ulp ) DSBEV('L','V', ... )
384 *> 51= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEV('L','N', ... )
385 *> 52= | A - U S U' | / ( |A| n ulp ) DSBEVX('L','V','A', ... )
386 *> 53= | I - U U' | / ( n ulp ) DSBEVX('L','V','A', ... )
387 *> 54= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVX('L','N','A', ... )
388 *> 55= | A - U S U' | / ( |A| n ulp ) DSBEVX('L','V','I', ... )
389 *> 56= | I - U U' | / ( n ulp ) DSBEVX('L','V','I', ... )
390 *> 57= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVX('L','N','I', ... )
391 *> 58= | A - U S U' | / ( |A| n ulp ) DSBEVX('L','V','V', ... )
392 *> 59= | I - U U' | / ( n ulp ) DSBEVX('L','V','V', ... )
393 *> 60= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVX('L','N','V', ... )
394 *> 61= | A - U S U' | / ( |A| n ulp ) DSYEVD('L','V', ... )
395 *> 62= | I - U U' | / ( n ulp ) DSYEVD('L','V', ... )
396 *> 63= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVD('L','N', ... )
397 *> 64= | A - U S U' | / ( |A| n ulp ) DSPEVD('L','V', ... )
398 *> 65= | I - U U' | / ( n ulp ) DSPEVD('L','V', ... )
399 *> 66= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVD('L','N', ... )
400 *> 67= | A - U S U' | / ( |A| n ulp ) DSBEVD('L','V', ... )
401 *> 68= | I - U U' | / ( n ulp ) DSBEVD('L','V', ... )
402 *> 69= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVD('L','N', ... )
403 *> 70= | A - U S U' | / ( |A| n ulp ) DSYEVR('L','V','A', ... )
404 *> 71= | I - U U' | / ( n ulp ) DSYEVR('L','V','A', ... )
405 *> 72= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVR('L','N','A', ... )
406 *> 73= | A - U S U' | / ( |A| n ulp ) DSYEVR('L','V','I', ... )
407 *> 74= | I - U U' | / ( n ulp ) DSYEVR('L','V','I', ... )
408 *> 75= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVR('L','N','I', ... )
409 *> 76= | A - U S U' | / ( |A| n ulp ) DSYEVR('L','V','V', ... )
410 *> 77= | I - U U' | / ( n ulp ) DSYEVR('L','V','V', ... )
411 *> 78= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVR('L','N','V', ... )
412 *>
413 *> Tests 25 through 78 are repeated (as tests 79 through 132)
414 *> with UPLO='U'
415 *>
416 *> To be added in 1999
417 *>
418 *> 79= | A - U S U' | / ( |A| n ulp ) DSPEVR('L','V','A', ... )
419 *> 80= | I - U U' | / ( n ulp ) DSPEVR('L','V','A', ... )
420 *> 81= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVR('L','N','A', ... )
421 *> 82= | A - U S U' | / ( |A| n ulp ) DSPEVR('L','V','I', ... )
422 *> 83= | I - U U' | / ( n ulp ) DSPEVR('L','V','I', ... )
423 *> 84= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVR('L','N','I', ... )
424 *> 85= | A - U S U' | / ( |A| n ulp ) DSPEVR('L','V','V', ... )
425 *> 86= | I - U U' | / ( n ulp ) DSPEVR('L','V','V', ... )
426 *> 87= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVR('L','N','V', ... )
427 *> 88= | A - U S U' | / ( |A| n ulp ) DSBEVR('L','V','A', ... )
428 *> 89= | I - U U' | / ( n ulp ) DSBEVR('L','V','A', ... )
429 *> 90= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVR('L','N','A', ... )
430 *> 91= | A - U S U' | / ( |A| n ulp ) DSBEVR('L','V','I', ... )
431 *> 92= | I - U U' | / ( n ulp ) DSBEVR('L','V','I', ... )
432 *> 93= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVR('L','N','I', ... )
433 *> 94= | A - U S U' | / ( |A| n ulp ) DSBEVR('L','V','V', ... )
434 *> 95= | I - U U' | / ( n ulp ) DSBEVR('L','V','V', ... )
435 *> 96= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVR('L','N','V', ... )
436 *> \endverbatim
437 *
438 * Authors:
439 * ========
440 *
441 *> \author Univ. of Tennessee
442 *> \author Univ. of California Berkeley
443 *> \author Univ. of Colorado Denver
444 *> \author NAG Ltd.
445 *
446 *> \date November 2011
447 *
448 *> \ingroup double_eig
449 *
450 * =====================================================================
451  SUBROUTINE ddrvst( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
452  $ nounit, a, lda, d1, d2, d3, d4, eveigs, wa1,
453  $ wa2, wa3, u, ldu, v, tau, z, work, lwork,
454  $ iwork, liwork, result, info )
455 *
456 * -- LAPACK test routine (version 3.4.0) --
457 * -- LAPACK is a software package provided by Univ. of Tennessee, --
458 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
459 * November 2011
460 *
461 * .. Scalar Arguments ..
462  INTEGER info, lda, ldu, liwork, lwork, nounit, nsizes,
463  $ ntypes
464  DOUBLE PRECISION thresh
465 * ..
466 * .. Array Arguments ..
467  LOGICAL dotype( * )
468  INTEGER iseed( 4 ), iwork( * ), nn( * )
469  DOUBLE PRECISION a( lda, * ), d1( * ), d2( * ), d3( * ),
470  $ d4( * ), eveigs( * ), result( * ), tau( * ),
471  $ u( ldu, * ), v( ldu, * ), wa1( * ), wa2( * ),
472  $ wa3( * ), work( * ), z( ldu, * )
473 * ..
474 *
475 * =====================================================================
476 *
477 * .. Parameters ..
478  DOUBLE PRECISION zero, one, two, ten
479  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
480  $ ten = 10.0d0 )
481  DOUBLE PRECISION half
482  parameter( half = 0.5d0 )
483  INTEGER maxtyp
484  parameter( maxtyp = 18 )
485 * ..
486 * .. Local Scalars ..
487  LOGICAL badnn
488  CHARACTER uplo
489  INTEGER i, idiag, ihbw, iinfo, il, imode, indx, irow,
490  $ itemp, itype, iu, iuplo, j, j1, j2, jcol,
491  $ jsize, jtype, kd, lgn, liwedc, lwedc, m, m2,
492  $ m3, mtypes, n, nerrs, nmats, nmax, ntest,
493  $ ntestt
494  DOUBLE PRECISION abstol, aninv, anorm, cond, ovfl, rtovfl,
495  $ rtunfl, temp1, temp2, temp3, ulp, ulpinv, unfl,
496  $ vl, vu
497 * ..
498 * .. Local Arrays ..
499  INTEGER idumma( 1 ), ioldsd( 4 ), iseed2( 4 ),
500  $ iseed3( 4 ), kmagn( maxtyp ), kmode( maxtyp ),
501  $ ktype( maxtyp )
502 * ..
503 * .. External Functions ..
504  DOUBLE PRECISION dlamch, dlarnd, dsxt1
505  EXTERNAL dlamch, dlarnd, dsxt1
506 * ..
507 * .. External Subroutines ..
508  EXTERNAL alasvm, dlabad, dlacpy, dlafts, dlaset, dlatmr,
509  $ dlatms, dsbev, dsbevd, dsbevx, dspev, dspevd,
510  $ dspevx, dstev, dstevd, dstevr, dstevx, dstt21,
511  $ dstt22, dsyev, dsyevd, dsyevr, dsyevx, dsyt21,
512  $ dsyt22, xerbla
513 * ..
514 * .. Scalars in Common ..
515  CHARACTER*32 srnamt
516 * ..
517 * .. Common blocks ..
518  common / srnamc / srnamt
519 * ..
520 * .. Intrinsic Functions ..
521  INTRINSIC abs, dble, int, log, max, min, sqrt
522 * ..
523 * .. Data statements ..
524  DATA ktype / 1, 2, 5*4, 5*5, 3*8, 3*9 /
525  DATA kmagn / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
526  $ 2, 3, 1, 2, 3 /
527  DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
528  $ 0, 0, 4, 4, 4 /
529 * ..
530 * .. Executable Statements ..
531 *
532 * Keep ftrnchek happy
533 *
534  vl = zero
535  vu = zero
536 *
537 * 1) Check for errors
538 *
539  ntestt = 0
540  info = 0
541 *
542  badnn = .false.
543  nmax = 1
544  DO 10 j = 1, nsizes
545  nmax = max( nmax, nn( j ) )
546  IF( nn( j ).LT.0 )
547  $ badnn = .true.
548  10 continue
549 *
550 * Check for errors
551 *
552  IF( nsizes.LT.0 ) THEN
553  info = -1
554  ELSE IF( badnn ) THEN
555  info = -2
556  ELSE IF( ntypes.LT.0 ) THEN
557  info = -3
558  ELSE IF( lda.LT.nmax ) THEN
559  info = -9
560  ELSE IF( ldu.LT.nmax ) THEN
561  info = -16
562  ELSE IF( 2*max( 2, nmax )**2.GT.lwork ) THEN
563  info = -21
564  END IF
565 *
566  IF( info.NE.0 ) THEN
567  CALL xerbla( 'DDRVST', -info )
568  return
569  END IF
570 *
571 * Quick return if nothing to do
572 *
573  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
574  $ return
575 *
576 * More Important constants
577 *
578  unfl = dlamch( 'Safe minimum' )
579  ovfl = dlamch( 'Overflow' )
580  CALL dlabad( unfl, ovfl )
581  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
582  ulpinv = one / ulp
583  rtunfl = sqrt( unfl )
584  rtovfl = sqrt( ovfl )
585 *
586 * Loop over sizes, types
587 *
588  DO 20 i = 1, 4
589  iseed2( i ) = iseed( i )
590  iseed3( i ) = iseed( i )
591  20 continue
592 *
593  nerrs = 0
594  nmats = 0
595 *
596 *
597  DO 1740 jsize = 1, nsizes
598  n = nn( jsize )
599  IF( n.GT.0 ) THEN
600  lgn = int( log( dble( n ) ) / log( two ) )
601  IF( 2**lgn.LT.n )
602  $ lgn = lgn + 1
603  IF( 2**lgn.LT.n )
604  $ lgn = lgn + 1
605  lwedc = 1 + 4*n + 2*n*lgn + 4*n**2
606 c LIWEDC = 6 + 6*N + 5*N*LGN
607  liwedc = 3 + 5*n
608  ELSE
609  lwedc = 9
610 c LIWEDC = 12
611  liwedc = 8
612  END IF
613  aninv = one / dble( max( 1, n ) )
614 *
615  IF( nsizes.NE.1 ) THEN
616  mtypes = min( maxtyp, ntypes )
617  ELSE
618  mtypes = min( maxtyp+1, ntypes )
619  END IF
620 *
621  DO 1730 jtype = 1, mtypes
622 *
623  IF( .NOT.dotype( jtype ) )
624  $ go to 1730
625  nmats = nmats + 1
626  ntest = 0
627 *
628  DO 30 j = 1, 4
629  ioldsd( j ) = iseed( j )
630  30 continue
631 *
632 * 2) Compute "A"
633 *
634 * Control parameters:
635 *
636 * KMAGN KMODE KTYPE
637 * =1 O(1) clustered 1 zero
638 * =2 large clustered 2 identity
639 * =3 small exponential (none)
640 * =4 arithmetic diagonal, (w/ eigenvalues)
641 * =5 random log symmetric, w/ eigenvalues
642 * =6 random (none)
643 * =7 random diagonal
644 * =8 random symmetric
645 * =9 band symmetric, w/ eigenvalues
646 *
647  IF( mtypes.GT.maxtyp )
648  $ go to 110
649 *
650  itype = ktype( jtype )
651  imode = kmode( jtype )
652 *
653 * Compute norm
654 *
655  go to( 40, 50, 60 )kmagn( jtype )
656 *
657  40 continue
658  anorm = one
659  go to 70
660 *
661  50 continue
662  anorm = ( rtovfl*ulp )*aninv
663  go to 70
664 *
665  60 continue
666  anorm = rtunfl*n*ulpinv
667  go to 70
668 *
669  70 continue
670 *
671  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
672  iinfo = 0
673  cond = ulpinv
674 *
675 * Special Matrices -- Identity & Jordan block
676 *
677 * Zero
678 *
679  IF( itype.EQ.1 ) THEN
680  iinfo = 0
681 *
682  ELSE IF( itype.EQ.2 ) THEN
683 *
684 * Identity
685 *
686  DO 80 jcol = 1, n
687  a( jcol, jcol ) = anorm
688  80 continue
689 *
690  ELSE IF( itype.EQ.4 ) THEN
691 *
692 * Diagonal Matrix, [Eigen]values Specified
693 *
694  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
695  $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
696  $ iinfo )
697 *
698  ELSE IF( itype.EQ.5 ) THEN
699 *
700 * Symmetric, eigenvalues specified
701 *
702  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
703  $ anorm, n, n, 'N', a, lda, work( n+1 ),
704  $ iinfo )
705 *
706  ELSE IF( itype.EQ.7 ) THEN
707 *
708 * Diagonal, random eigenvalues
709 *
710  idumma( 1 ) = 1
711  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
712  $ 'T', 'N', work( n+1 ), 1, one,
713  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
714  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
715 *
716  ELSE IF( itype.EQ.8 ) THEN
717 *
718 * Symmetric, random eigenvalues
719 *
720  idumma( 1 ) = 1
721  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
722  $ 'T', 'N', work( n+1 ), 1, one,
723  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
724  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
725 *
726  ELSE IF( itype.EQ.9 ) THEN
727 *
728 * Symmetric banded, eigenvalues specified
729 *
730  ihbw = int( ( n-1 )*dlarnd( 1, iseed3 ) )
731  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
732  $ anorm, ihbw, ihbw, 'Z', u, ldu, work( n+1 ),
733  $ iinfo )
734 *
735 * Store as dense matrix for most routines.
736 *
737  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
738  DO 100 idiag = -ihbw, ihbw
739  irow = ihbw - idiag + 1
740  j1 = max( 1, idiag+1 )
741  j2 = min( n, n+idiag )
742  DO 90 j = j1, j2
743  i = j - idiag
744  a( i, j ) = u( irow, j )
745  90 continue
746  100 continue
747  ELSE
748  iinfo = 1
749  END IF
750 *
751  IF( iinfo.NE.0 ) THEN
752  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
753  $ ioldsd
754  info = abs( iinfo )
755  return
756  END IF
757 *
758  110 continue
759 *
760  abstol = unfl + unfl
761  IF( n.LE.1 ) THEN
762  il = 1
763  iu = n
764  ELSE
765  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
766  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
767  IF( il.GT.iu ) THEN
768  itemp = il
769  il = iu
770  iu = itemp
771  END IF
772  END IF
773 *
774 * 3) If matrix is tridiagonal, call DSTEV and DSTEVX.
775 *
776  IF( jtype.LE.7 ) THEN
777  ntest = 1
778  DO 120 i = 1, n
779  d1( i ) = dble( a( i, i ) )
780  120 continue
781  DO 130 i = 1, n - 1
782  d2( i ) = dble( a( i+1, i ) )
783  130 continue
784  srnamt = 'DSTEV'
785  CALL dstev( 'V', n, d1, d2, z, ldu, work, iinfo )
786  IF( iinfo.NE.0 ) THEN
787  WRITE( nounit, fmt = 9999 )'DSTEV(V)', iinfo, n,
788  $ jtype, ioldsd
789  info = abs( iinfo )
790  IF( iinfo.LT.0 ) THEN
791  return
792  ELSE
793  result( 1 ) = ulpinv
794  result( 2 ) = ulpinv
795  result( 3 ) = ulpinv
796  go to 180
797  END IF
798  END IF
799 *
800 * Do tests 1 and 2.
801 *
802  DO 140 i = 1, n
803  d3( i ) = dble( a( i, i ) )
804  140 continue
805  DO 150 i = 1, n - 1
806  d4( i ) = dble( a( i+1, i ) )
807  150 continue
808  CALL dstt21( n, 0, d3, d4, d1, d2, z, ldu, work,
809  $ result( 1 ) )
810 *
811  ntest = 3
812  DO 160 i = 1, n - 1
813  d4( i ) = dble( a( i+1, i ) )
814  160 continue
815  srnamt = 'DSTEV'
816  CALL dstev( 'N', n, d3, d4, z, ldu, work, iinfo )
817  IF( iinfo.NE.0 ) THEN
818  WRITE( nounit, fmt = 9999 )'DSTEV(N)', iinfo, n,
819  $ jtype, ioldsd
820  info = abs( iinfo )
821  IF( iinfo.LT.0 ) THEN
822  return
823  ELSE
824  result( 3 ) = ulpinv
825  go to 180
826  END IF
827  END IF
828 *
829 * Do test 3.
830 *
831  temp1 = zero
832  temp2 = zero
833  DO 170 j = 1, n
834  temp1 = max( temp1, abs( d1( j ) ), abs( d3( j ) ) )
835  temp2 = max( temp2, abs( d1( j )-d3( j ) ) )
836  170 continue
837  result( 3 ) = temp2 / max( unfl,
838  $ ulp*max( temp1, temp2 ) )
839 *
840  180 continue
841 *
842  ntest = 4
843  DO 190 i = 1, n
844  eveigs( i ) = d3( i )
845  d1( i ) = dble( a( i, i ) )
846  190 continue
847  DO 200 i = 1, n - 1
848  d2( i ) = dble( a( i+1, i ) )
849  200 continue
850  srnamt = 'DSTEVX'
851  CALL dstevx( 'V', 'A', n, d1, d2, vl, vu, il, iu, abstol,
852  $ m, wa1, z, ldu, work, iwork, iwork( 5*n+1 ),
853  $ iinfo )
854  IF( iinfo.NE.0 ) THEN
855  WRITE( nounit, fmt = 9999 )'DSTEVX(V,A)', iinfo, n,
856  $ jtype, ioldsd
857  info = abs( iinfo )
858  IF( iinfo.LT.0 ) THEN
859  return
860  ELSE
861  result( 4 ) = ulpinv
862  result( 5 ) = ulpinv
863  result( 6 ) = ulpinv
864  go to 250
865  END IF
866  END IF
867  IF( n.GT.0 ) THEN
868  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
869  ELSE
870  temp3 = zero
871  END IF
872 *
873 * Do tests 4 and 5.
874 *
875  DO 210 i = 1, n
876  d3( i ) = dble( a( i, i ) )
877  210 continue
878  DO 220 i = 1, n - 1
879  d4( i ) = dble( a( i+1, i ) )
880  220 continue
881  CALL dstt21( n, 0, d3, d4, wa1, d2, z, ldu, work,
882  $ result( 4 ) )
883 *
884  ntest = 6
885  DO 230 i = 1, n - 1
886  d4( i ) = dble( a( i+1, i ) )
887  230 continue
888  srnamt = 'DSTEVX'
889  CALL dstevx( 'N', 'A', n, d3, d4, vl, vu, il, iu, abstol,
890  $ m2, wa2, z, ldu, work, iwork,
891  $ iwork( 5*n+1 ), iinfo )
892  IF( iinfo.NE.0 ) THEN
893  WRITE( nounit, fmt = 9999 )'DSTEVX(N,A)', iinfo, n,
894  $ jtype, ioldsd
895  info = abs( iinfo )
896  IF( iinfo.LT.0 ) THEN
897  return
898  ELSE
899  result( 6 ) = ulpinv
900  go to 250
901  END IF
902  END IF
903 *
904 * Do test 6.
905 *
906  temp1 = zero
907  temp2 = zero
908  DO 240 j = 1, n
909  temp1 = max( temp1, abs( wa2( j ) ),
910  $ abs( eveigs( j ) ) )
911  temp2 = max( temp2, abs( wa2( j )-eveigs( j ) ) )
912  240 continue
913  result( 6 ) = temp2 / max( unfl,
914  $ ulp*max( temp1, temp2 ) )
915 *
916  250 continue
917 *
918  ntest = 7
919  DO 260 i = 1, n
920  d1( i ) = dble( a( i, i ) )
921  260 continue
922  DO 270 i = 1, n - 1
923  d2( i ) = dble( a( i+1, i ) )
924  270 continue
925  srnamt = 'DSTEVR'
926  CALL dstevr( 'V', 'A', n, d1, d2, vl, vu, il, iu, abstol,
927  $ m, wa1, z, ldu, iwork, work, lwork,
928  $ iwork(2*n+1), liwork-2*n, iinfo )
929  IF( iinfo.NE.0 ) THEN
930  WRITE( nounit, fmt = 9999 )'DSTEVR(V,A)', iinfo, n,
931  $ jtype, ioldsd
932  info = abs( iinfo )
933  IF( iinfo.LT.0 ) THEN
934  return
935  ELSE
936  result( 7 ) = ulpinv
937  result( 8 ) = ulpinv
938  go to 320
939  END IF
940  END IF
941  IF( n.GT.0 ) THEN
942  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
943  ELSE
944  temp3 = zero
945  END IF
946 *
947 * Do tests 7 and 8.
948 *
949  DO 280 i = 1, n
950  d3( i ) = dble( a( i, i ) )
951  280 continue
952  DO 290 i = 1, n - 1
953  d4( i ) = dble( a( i+1, i ) )
954  290 continue
955  CALL dstt21( n, 0, d3, d4, wa1, d2, z, ldu, work,
956  $ result( 7 ) )
957 *
958  ntest = 9
959  DO 300 i = 1, n - 1
960  d4( i ) = dble( a( i+1, i ) )
961  300 continue
962  srnamt = 'DSTEVR'
963  CALL dstevr( 'N', 'A', n, d3, d4, vl, vu, il, iu, abstol,
964  $ m2, wa2, z, ldu, iwork, work, lwork,
965  $ iwork(2*n+1), liwork-2*n, iinfo )
966  IF( iinfo.NE.0 ) THEN
967  WRITE( nounit, fmt = 9999 )'DSTEVR(N,A)', iinfo, n,
968  $ jtype, ioldsd
969  info = abs( iinfo )
970  IF( iinfo.LT.0 ) THEN
971  return
972  ELSE
973  result( 9 ) = ulpinv
974  go to 320
975  END IF
976  END IF
977 *
978 * Do test 9.
979 *
980  temp1 = zero
981  temp2 = zero
982  DO 310 j = 1, n
983  temp1 = max( temp1, abs( wa2( j ) ),
984  $ abs( eveigs( j ) ) )
985  temp2 = max( temp2, abs( wa2( j )-eveigs( j ) ) )
986  310 continue
987  result( 9 ) = temp2 / max( unfl,
988  $ ulp*max( temp1, temp2 ) )
989 *
990  320 continue
991 *
992 *
993  ntest = 10
994  DO 330 i = 1, n
995  d1( i ) = dble( a( i, i ) )
996  330 continue
997  DO 340 i = 1, n - 1
998  d2( i ) = dble( a( i+1, i ) )
999  340 continue
1000  srnamt = 'DSTEVX'
1001  CALL dstevx( 'V', 'I', n, d1, d2, vl, vu, il, iu, abstol,
1002  $ m2, wa2, z, ldu, work, iwork,
1003  $ iwork( 5*n+1 ), iinfo )
1004  IF( iinfo.NE.0 ) THEN
1005  WRITE( nounit, fmt = 9999 )'DSTEVX(V,I)', iinfo, n,
1006  $ jtype, ioldsd
1007  info = abs( iinfo )
1008  IF( iinfo.LT.0 ) THEN
1009  return
1010  ELSE
1011  result( 10 ) = ulpinv
1012  result( 11 ) = ulpinv
1013  result( 12 ) = ulpinv
1014  go to 380
1015  END IF
1016  END IF
1017 *
1018 * Do tests 10 and 11.
1019 *
1020  DO 350 i = 1, n
1021  d3( i ) = dble( a( i, i ) )
1022  350 continue
1023  DO 360 i = 1, n - 1
1024  d4( i ) = dble( a( i+1, i ) )
1025  360 continue
1026  CALL dstt22( n, m2, 0, d3, d4, wa2, d2, z, ldu, work,
1027  $ max( 1, m2 ), result( 10 ) )
1028 *
1029 *
1030  ntest = 12
1031  DO 370 i = 1, n - 1
1032  d4( i ) = dble( a( i+1, i ) )
1033  370 continue
1034  srnamt = 'DSTEVX'
1035  CALL dstevx( 'N', 'I', n, d3, d4, vl, vu, il, iu, abstol,
1036  $ m3, wa3, z, ldu, work, iwork,
1037  $ iwork( 5*n+1 ), iinfo )
1038  IF( iinfo.NE.0 ) THEN
1039  WRITE( nounit, fmt = 9999 )'DSTEVX(N,I)', iinfo, n,
1040  $ jtype, ioldsd
1041  info = abs( iinfo )
1042  IF( iinfo.LT.0 ) THEN
1043  return
1044  ELSE
1045  result( 12 ) = ulpinv
1046  go to 380
1047  END IF
1048  END IF
1049 *
1050 * Do test 12.
1051 *
1052  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1053  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1054  result( 12 ) = ( temp1+temp2 ) / max( unfl, ulp*temp3 )
1055 *
1056  380 continue
1057 *
1058  ntest = 12
1059  IF( n.GT.0 ) THEN
1060  IF( il.NE.1 ) THEN
1061  vl = wa1( il ) - max( half*
1062  $ ( wa1( il )-wa1( il-1 ) ), ten*ulp*temp3,
1063  $ ten*rtunfl )
1064  ELSE
1065  vl = wa1( 1 ) - max( half*( wa1( n )-wa1( 1 ) ),
1066  $ ten*ulp*temp3, ten*rtunfl )
1067  END IF
1068  IF( iu.NE.n ) THEN
1069  vu = wa1( iu ) + max( half*
1070  $ ( wa1( iu+1 )-wa1( iu ) ), ten*ulp*temp3,
1071  $ ten*rtunfl )
1072  ELSE
1073  vu = wa1( n ) + max( half*( wa1( n )-wa1( 1 ) ),
1074  $ ten*ulp*temp3, ten*rtunfl )
1075  END IF
1076  ELSE
1077  vl = zero
1078  vu = one
1079  END IF
1080 *
1081  DO 390 i = 1, n
1082  d1( i ) = dble( a( i, i ) )
1083  390 continue
1084  DO 400 i = 1, n - 1
1085  d2( i ) = dble( a( i+1, i ) )
1086  400 continue
1087  srnamt = 'DSTEVX'
1088  CALL dstevx( 'V', 'V', n, d1, d2, vl, vu, il, iu, abstol,
1089  $ m2, wa2, z, ldu, work, iwork,
1090  $ iwork( 5*n+1 ), iinfo )
1091  IF( iinfo.NE.0 ) THEN
1092  WRITE( nounit, fmt = 9999 )'DSTEVX(V,V)', iinfo, n,
1093  $ jtype, ioldsd
1094  info = abs( iinfo )
1095  IF( iinfo.LT.0 ) THEN
1096  return
1097  ELSE
1098  result( 13 ) = ulpinv
1099  result( 14 ) = ulpinv
1100  result( 15 ) = ulpinv
1101  go to 440
1102  END IF
1103  END IF
1104 *
1105  IF( m2.EQ.0 .AND. n.GT.0 ) THEN
1106  result( 13 ) = ulpinv
1107  result( 14 ) = ulpinv
1108  result( 15 ) = ulpinv
1109  go to 440
1110  END IF
1111 *
1112 * Do tests 13 and 14.
1113 *
1114  DO 410 i = 1, n
1115  d3( i ) = dble( a( i, i ) )
1116  410 continue
1117  DO 420 i = 1, n - 1
1118  d4( i ) = dble( a( i+1, i ) )
1119  420 continue
1120  CALL dstt22( n, m2, 0, d3, d4, wa2, d2, z, ldu, work,
1121  $ max( 1, m2 ), result( 13 ) )
1122 *
1123  ntest = 15
1124  DO 430 i = 1, n - 1
1125  d4( i ) = dble( a( i+1, i ) )
1126  430 continue
1127  srnamt = 'DSTEVX'
1128  CALL dstevx( 'N', 'V', n, d3, d4, vl, vu, il, iu, abstol,
1129  $ m3, wa3, z, ldu, work, iwork,
1130  $ iwork( 5*n+1 ), iinfo )
1131  IF( iinfo.NE.0 ) THEN
1132  WRITE( nounit, fmt = 9999 )'DSTEVX(N,V)', iinfo, n,
1133  $ jtype, ioldsd
1134  info = abs( iinfo )
1135  IF( iinfo.LT.0 ) THEN
1136  return
1137  ELSE
1138  result( 15 ) = ulpinv
1139  go to 440
1140  END IF
1141  END IF
1142 *
1143 * Do test 15.
1144 *
1145  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1146  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1147  result( 15 ) = ( temp1+temp2 ) / max( unfl, temp3*ulp )
1148 *
1149  440 continue
1150 *
1151  ntest = 16
1152  DO 450 i = 1, n
1153  d1( i ) = dble( a( i, i ) )
1154  450 continue
1155  DO 460 i = 1, n - 1
1156  d2( i ) = dble( a( i+1, i ) )
1157  460 continue
1158  srnamt = 'DSTEVD'
1159  CALL dstevd( 'V', n, d1, d2, z, ldu, work, lwedc, iwork,
1160  $ liwedc, iinfo )
1161  IF( iinfo.NE.0 ) THEN
1162  WRITE( nounit, fmt = 9999 )'DSTEVD(V)', iinfo, n,
1163  $ jtype, ioldsd
1164  info = abs( iinfo )
1165  IF( iinfo.LT.0 ) THEN
1166  return
1167  ELSE
1168  result( 16 ) = ulpinv
1169  result( 17 ) = ulpinv
1170  result( 18 ) = ulpinv
1171  go to 510
1172  END IF
1173  END IF
1174 *
1175 * Do tests 16 and 17.
1176 *
1177  DO 470 i = 1, n
1178  d3( i ) = dble( a( i, i ) )
1179  470 continue
1180  DO 480 i = 1, n - 1
1181  d4( i ) = dble( a( i+1, i ) )
1182  480 continue
1183  CALL dstt21( n, 0, d3, d4, d1, d2, z, ldu, work,
1184  $ result( 16 ) )
1185 *
1186  ntest = 18
1187  DO 490 i = 1, n - 1
1188  d4( i ) = dble( a( i+1, i ) )
1189  490 continue
1190  srnamt = 'DSTEVD'
1191  CALL dstevd( 'N', n, d3, d4, z, ldu, work, lwedc, iwork,
1192  $ liwedc, iinfo )
1193  IF( iinfo.NE.0 ) THEN
1194  WRITE( nounit, fmt = 9999 )'DSTEVD(N)', iinfo, n,
1195  $ jtype, ioldsd
1196  info = abs( iinfo )
1197  IF( iinfo.LT.0 ) THEN
1198  return
1199  ELSE
1200  result( 18 ) = ulpinv
1201  go to 510
1202  END IF
1203  END IF
1204 *
1205 * Do test 18.
1206 *
1207  temp1 = zero
1208  temp2 = zero
1209  DO 500 j = 1, n
1210  temp1 = max( temp1, abs( eveigs( j ) ),
1211  $ abs( d3( j ) ) )
1212  temp2 = max( temp2, abs( eveigs( j )-d3( j ) ) )
1213  500 continue
1214  result( 18 ) = temp2 / max( unfl,
1215  $ ulp*max( temp1, temp2 ) )
1216 *
1217  510 continue
1218 *
1219  ntest = 19
1220  DO 520 i = 1, n
1221  d1( i ) = dble( a( i, i ) )
1222  520 continue
1223  DO 530 i = 1, n - 1
1224  d2( i ) = dble( a( i+1, i ) )
1225  530 continue
1226  srnamt = 'DSTEVR'
1227  CALL dstevr( 'V', 'I', n, d1, d2, vl, vu, il, iu, abstol,
1228  $ m2, wa2, z, ldu, iwork, work, lwork,
1229  $ iwork(2*n+1), liwork-2*n, iinfo )
1230  IF( iinfo.NE.0 ) THEN
1231  WRITE( nounit, fmt = 9999 )'DSTEVR(V,I)', iinfo, n,
1232  $ jtype, ioldsd
1233  info = abs( iinfo )
1234  IF( iinfo.LT.0 ) THEN
1235  return
1236  ELSE
1237  result( 19 ) = ulpinv
1238  result( 20 ) = ulpinv
1239  result( 21 ) = ulpinv
1240  go to 570
1241  END IF
1242  END IF
1243 *
1244 * DO tests 19 and 20.
1245 *
1246  DO 540 i = 1, n
1247  d3( i ) = dble( a( i, i ) )
1248  540 continue
1249  DO 550 i = 1, n - 1
1250  d4( i ) = dble( a( i+1, i ) )
1251  550 continue
1252  CALL dstt22( n, m2, 0, d3, d4, wa2, d2, z, ldu, work,
1253  $ max( 1, m2 ), result( 19 ) )
1254 *
1255 *
1256  ntest = 21
1257  DO 560 i = 1, n - 1
1258  d4( i ) = dble( a( i+1, i ) )
1259  560 continue
1260  srnamt = 'DSTEVR'
1261  CALL dstevr( 'N', 'I', n, d3, d4, vl, vu, il, iu, abstol,
1262  $ m3, wa3, z, ldu, iwork, work, lwork,
1263  $ iwork(2*n+1), liwork-2*n, iinfo )
1264  IF( iinfo.NE.0 ) THEN
1265  WRITE( nounit, fmt = 9999 )'DSTEVR(N,I)', iinfo, n,
1266  $ jtype, ioldsd
1267  info = abs( iinfo )
1268  IF( iinfo.LT.0 ) THEN
1269  return
1270  ELSE
1271  result( 21 ) = ulpinv
1272  go to 570
1273  END IF
1274  END IF
1275 *
1276 * Do test 21.
1277 *
1278  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1279  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1280  result( 21 ) = ( temp1+temp2 ) / max( unfl, ulp*temp3 )
1281 *
1282  570 continue
1283 *
1284  ntest = 21
1285  IF( n.GT.0 ) THEN
1286  IF( il.NE.1 ) THEN
1287  vl = wa1( il ) - max( half*
1288  $ ( wa1( il )-wa1( il-1 ) ), ten*ulp*temp3,
1289  $ ten*rtunfl )
1290  ELSE
1291  vl = wa1( 1 ) - max( half*( wa1( n )-wa1( 1 ) ),
1292  $ ten*ulp*temp3, ten*rtunfl )
1293  END IF
1294  IF( iu.NE.n ) THEN
1295  vu = wa1( iu ) + max( half*
1296  $ ( wa1( iu+1 )-wa1( iu ) ), ten*ulp*temp3,
1297  $ ten*rtunfl )
1298  ELSE
1299  vu = wa1( n ) + max( half*( wa1( n )-wa1( 1 ) ),
1300  $ ten*ulp*temp3, ten*rtunfl )
1301  END IF
1302  ELSE
1303  vl = zero
1304  vu = one
1305  END IF
1306 *
1307  DO 580 i = 1, n
1308  d1( i ) = dble( a( i, i ) )
1309  580 continue
1310  DO 590 i = 1, n - 1
1311  d2( i ) = dble( a( i+1, i ) )
1312  590 continue
1313  srnamt = 'DSTEVR'
1314  CALL dstevr( 'V', 'V', n, d1, d2, vl, vu, il, iu, abstol,
1315  $ m2, wa2, z, ldu, iwork, work, lwork,
1316  $ iwork(2*n+1), liwork-2*n, iinfo )
1317  IF( iinfo.NE.0 ) THEN
1318  WRITE( nounit, fmt = 9999 )'DSTEVR(V,V)', iinfo, n,
1319  $ jtype, ioldsd
1320  info = abs( iinfo )
1321  IF( iinfo.LT.0 ) THEN
1322  return
1323  ELSE
1324  result( 22 ) = ulpinv
1325  result( 23 ) = ulpinv
1326  result( 24 ) = ulpinv
1327  go to 630
1328  END IF
1329  END IF
1330 *
1331  IF( m2.EQ.0 .AND. n.GT.0 ) THEN
1332  result( 22 ) = ulpinv
1333  result( 23 ) = ulpinv
1334  result( 24 ) = ulpinv
1335  go to 630
1336  END IF
1337 *
1338 * Do tests 22 and 23.
1339 *
1340  DO 600 i = 1, n
1341  d3( i ) = dble( a( i, i ) )
1342  600 continue
1343  DO 610 i = 1, n - 1
1344  d4( i ) = dble( a( i+1, i ) )
1345  610 continue
1346  CALL dstt22( n, m2, 0, d3, d4, wa2, d2, z, ldu, work,
1347  $ max( 1, m2 ), result( 22 ) )
1348 *
1349  ntest = 24
1350  DO 620 i = 1, n - 1
1351  d4( i ) = dble( a( i+1, i ) )
1352  620 continue
1353  srnamt = 'DSTEVR'
1354  CALL dstevr( 'N', 'V', n, d3, d4, vl, vu, il, iu, abstol,
1355  $ m3, wa3, z, ldu, iwork, work, lwork,
1356  $ iwork(2*n+1), liwork-2*n, iinfo )
1357  IF( iinfo.NE.0 ) THEN
1358  WRITE( nounit, fmt = 9999 )'DSTEVR(N,V)', iinfo, n,
1359  $ jtype, ioldsd
1360  info = abs( iinfo )
1361  IF( iinfo.LT.0 ) THEN
1362  return
1363  ELSE
1364  result( 24 ) = ulpinv
1365  go to 630
1366  END IF
1367  END IF
1368 *
1369 * Do test 24.
1370 *
1371  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1372  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1373  result( 24 ) = ( temp1+temp2 ) / max( unfl, temp3*ulp )
1374 *
1375  630 continue
1376 *
1377 *
1378 *
1379  ELSE
1380 *
1381  DO 640 i = 1, 24
1382  result( i ) = zero
1383  640 continue
1384  ntest = 24
1385  END IF
1386 *
1387 * Perform remaining tests storing upper or lower triangular
1388 * part of matrix.
1389 *
1390  DO 1720 iuplo = 0, 1
1391  IF( iuplo.EQ.0 ) THEN
1392  uplo = 'L'
1393  ELSE
1394  uplo = 'U'
1395  END IF
1396 *
1397 * 4) Call DSYEV and DSYEVX.
1398 *
1399  CALL dlacpy( ' ', n, n, a, lda, v, ldu )
1400 *
1401  ntest = ntest + 1
1402  srnamt = 'DSYEV'
1403  CALL dsyev( 'V', uplo, n, a, ldu, d1, work, lwork,
1404  $ iinfo )
1405  IF( iinfo.NE.0 ) THEN
1406  WRITE( nounit, fmt = 9999 )'DSYEV(V,' // uplo // ')',
1407  $ iinfo, n, jtype, ioldsd
1408  info = abs( iinfo )
1409  IF( iinfo.LT.0 ) THEN
1410  return
1411  ELSE
1412  result( ntest ) = ulpinv
1413  result( ntest+1 ) = ulpinv
1414  result( ntest+2 ) = ulpinv
1415  go to 660
1416  END IF
1417  END IF
1418 *
1419 * Do tests 25 and 26 (or +54)
1420 *
1421  CALL dsyt21( 1, uplo, n, 0, v, ldu, d1, d2, a, ldu, z,
1422  $ ldu, tau, work, result( ntest ) )
1423 *
1424  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1425 *
1426  ntest = ntest + 2
1427  srnamt = 'DSYEV'
1428  CALL dsyev( 'N', uplo, n, a, ldu, d3, work, lwork,
1429  $ iinfo )
1430  IF( iinfo.NE.0 ) THEN
1431  WRITE( nounit, fmt = 9999 )'DSYEV(N,' // uplo // ')',
1432  $ iinfo, n, jtype, ioldsd
1433  info = abs( iinfo )
1434  IF( iinfo.LT.0 ) THEN
1435  return
1436  ELSE
1437  result( ntest ) = ulpinv
1438  go to 660
1439  END IF
1440  END IF
1441 *
1442 * Do test 27 (or +54)
1443 *
1444  temp1 = zero
1445  temp2 = zero
1446  DO 650 j = 1, n
1447  temp1 = max( temp1, abs( d1( j ) ), abs( d3( j ) ) )
1448  temp2 = max( temp2, abs( d1( j )-d3( j ) ) )
1449  650 continue
1450  result( ntest ) = temp2 / max( unfl,
1451  $ ulp*max( temp1, temp2 ) )
1452 *
1453  660 continue
1454  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1455 *
1456  ntest = ntest + 1
1457 *
1458  IF( n.GT.0 ) THEN
1459  temp3 = max( abs( d1( 1 ) ), abs( d1( n ) ) )
1460  IF( il.NE.1 ) THEN
1461  vl = d1( il ) - max( half*( d1( il )-d1( il-1 ) ),
1462  $ ten*ulp*temp3, ten*rtunfl )
1463  ELSE IF( n.GT.0 ) THEN
1464  vl = d1( 1 ) - max( half*( d1( n )-d1( 1 ) ),
1465  $ ten*ulp*temp3, ten*rtunfl )
1466  END IF
1467  IF( iu.NE.n ) THEN
1468  vu = d1( iu ) + max( half*( d1( iu+1 )-d1( iu ) ),
1469  $ ten*ulp*temp3, ten*rtunfl )
1470  ELSE IF( n.GT.0 ) THEN
1471  vu = d1( n ) + max( half*( d1( n )-d1( 1 ) ),
1472  $ ten*ulp*temp3, ten*rtunfl )
1473  END IF
1474  ELSE
1475  temp3 = zero
1476  vl = zero
1477  vu = one
1478  END IF
1479 *
1480  srnamt = 'DSYEVX'
1481  CALL dsyevx( 'V', 'A', uplo, n, a, ldu, vl, vu, il, iu,
1482  $ abstol, m, wa1, z, ldu, work, lwork, iwork,
1483  $ iwork( 5*n+1 ), iinfo )
1484  IF( iinfo.NE.0 ) THEN
1485  WRITE( nounit, fmt = 9999 )'DSYEVX(V,A,' // uplo //
1486  $ ')', iinfo, n, jtype, ioldsd
1487  info = abs( iinfo )
1488  IF( iinfo.LT.0 ) THEN
1489  return
1490  ELSE
1491  result( ntest ) = ulpinv
1492  result( ntest+1 ) = ulpinv
1493  result( ntest+2 ) = ulpinv
1494  go to 680
1495  END IF
1496  END IF
1497 *
1498 * Do tests 28 and 29 (or +54)
1499 *
1500  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1501 *
1502  CALL dsyt21( 1, uplo, n, 0, a, ldu, d1, d2, z, ldu, v,
1503  $ ldu, tau, work, result( ntest ) )
1504 *
1505  ntest = ntest + 2
1506  srnamt = 'DSYEVX'
1507  CALL dsyevx( 'N', 'A', uplo, n, a, ldu, vl, vu, il, iu,
1508  $ abstol, m2, wa2, z, ldu, work, lwork, iwork,
1509  $ iwork( 5*n+1 ), iinfo )
1510  IF( iinfo.NE.0 ) THEN
1511  WRITE( nounit, fmt = 9999 )'DSYEVX(N,A,' // uplo //
1512  $ ')', iinfo, n, jtype, ioldsd
1513  info = abs( iinfo )
1514  IF( iinfo.LT.0 ) THEN
1515  return
1516  ELSE
1517  result( ntest ) = ulpinv
1518  go to 680
1519  END IF
1520  END IF
1521 *
1522 * Do test 30 (or +54)
1523 *
1524  temp1 = zero
1525  temp2 = zero
1526  DO 670 j = 1, n
1527  temp1 = max( temp1, abs( wa1( j ) ), abs( wa2( j ) ) )
1528  temp2 = max( temp2, abs( wa1( j )-wa2( j ) ) )
1529  670 continue
1530  result( ntest ) = temp2 / max( unfl,
1531  $ ulp*max( temp1, temp2 ) )
1532 *
1533  680 continue
1534 *
1535  ntest = ntest + 1
1536  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1537  srnamt = 'DSYEVX'
1538  CALL dsyevx( 'V', 'I', uplo, n, a, ldu, vl, vu, il, iu,
1539  $ abstol, m2, wa2, z, ldu, work, lwork, iwork,
1540  $ iwork( 5*n+1 ), iinfo )
1541  IF( iinfo.NE.0 ) THEN
1542  WRITE( nounit, fmt = 9999 )'DSYEVX(V,I,' // uplo //
1543  $ ')', iinfo, n, jtype, ioldsd
1544  info = abs( iinfo )
1545  IF( iinfo.LT.0 ) THEN
1546  return
1547  ELSE
1548  result( ntest ) = ulpinv
1549  result( ntest+1 ) = ulpinv
1550  result( ntest+2 ) = ulpinv
1551  go to 690
1552  END IF
1553  END IF
1554 *
1555 * Do tests 31 and 32 (or +54)
1556 *
1557  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1558 *
1559  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
1560  $ v, ldu, tau, work, result( ntest ) )
1561 *
1562  ntest = ntest + 2
1563  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1564  srnamt = 'DSYEVX'
1565  CALL dsyevx( 'N', 'I', uplo, n, a, ldu, vl, vu, il, iu,
1566  $ abstol, m3, wa3, z, ldu, work, lwork, iwork,
1567  $ iwork( 5*n+1 ), iinfo )
1568  IF( iinfo.NE.0 ) THEN
1569  WRITE( nounit, fmt = 9999 )'DSYEVX(N,I,' // uplo //
1570  $ ')', iinfo, n, jtype, ioldsd
1571  info = abs( iinfo )
1572  IF( iinfo.LT.0 ) THEN
1573  return
1574  ELSE
1575  result( ntest ) = ulpinv
1576  go to 690
1577  END IF
1578  END IF
1579 *
1580 * Do test 33 (or +54)
1581 *
1582  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1583  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1584  result( ntest ) = ( temp1+temp2 ) /
1585  $ max( unfl, ulp*temp3 )
1586  690 continue
1587 *
1588  ntest = ntest + 1
1589  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1590  srnamt = 'DSYEVX'
1591  CALL dsyevx( 'V', 'V', uplo, n, a, ldu, vl, vu, il, iu,
1592  $ abstol, m2, wa2, z, ldu, work, lwork, iwork,
1593  $ iwork( 5*n+1 ), iinfo )
1594  IF( iinfo.NE.0 ) THEN
1595  WRITE( nounit, fmt = 9999 )'DSYEVX(V,V,' // uplo //
1596  $ ')', iinfo, n, jtype, ioldsd
1597  info = abs( iinfo )
1598  IF( iinfo.LT.0 ) THEN
1599  return
1600  ELSE
1601  result( ntest ) = ulpinv
1602  result( ntest+1 ) = ulpinv
1603  result( ntest+2 ) = ulpinv
1604  go to 700
1605  END IF
1606  END IF
1607 *
1608 * Do tests 34 and 35 (or +54)
1609 *
1610  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1611 *
1612  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
1613  $ v, ldu, tau, work, result( ntest ) )
1614 *
1615  ntest = ntest + 2
1616  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1617  srnamt = 'DSYEVX'
1618  CALL dsyevx( 'N', 'V', uplo, n, a, ldu, vl, vu, il, iu,
1619  $ abstol, m3, wa3, z, ldu, work, lwork, iwork,
1620  $ iwork( 5*n+1 ), iinfo )
1621  IF( iinfo.NE.0 ) THEN
1622  WRITE( nounit, fmt = 9999 )'DSYEVX(N,V,' // uplo //
1623  $ ')', iinfo, n, jtype, ioldsd
1624  info = abs( iinfo )
1625  IF( iinfo.LT.0 ) THEN
1626  return
1627  ELSE
1628  result( ntest ) = ulpinv
1629  go to 700
1630  END IF
1631  END IF
1632 *
1633  IF( m3.EQ.0 .AND. n.GT.0 ) THEN
1634  result( ntest ) = ulpinv
1635  go to 700
1636  END IF
1637 *
1638 * Do test 36 (or +54)
1639 *
1640  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1641  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1642  IF( n.GT.0 ) THEN
1643  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
1644  ELSE
1645  temp3 = zero
1646  END IF
1647  result( ntest ) = ( temp1+temp2 ) /
1648  $ max( unfl, temp3*ulp )
1649 *
1650  700 continue
1651 *
1652 * 5) Call DSPEV and DSPEVX.
1653 *
1654  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
1655 *
1656 * Load array WORK with the upper or lower triangular
1657 * part of the matrix in packed form.
1658 *
1659  IF( iuplo.EQ.1 ) THEN
1660  indx = 1
1661  DO 720 j = 1, n
1662  DO 710 i = 1, j
1663  work( indx ) = a( i, j )
1664  indx = indx + 1
1665  710 continue
1666  720 continue
1667  ELSE
1668  indx = 1
1669  DO 740 j = 1, n
1670  DO 730 i = j, n
1671  work( indx ) = a( i, j )
1672  indx = indx + 1
1673  730 continue
1674  740 continue
1675  END IF
1676 *
1677  ntest = ntest + 1
1678  srnamt = 'DSPEV'
1679  CALL dspev( 'V', uplo, n, work, d1, z, ldu, v, iinfo )
1680  IF( iinfo.NE.0 ) THEN
1681  WRITE( nounit, fmt = 9999 )'DSPEV(V,' // uplo // ')',
1682  $ iinfo, n, jtype, ioldsd
1683  info = abs( iinfo )
1684  IF( iinfo.LT.0 ) THEN
1685  return
1686  ELSE
1687  result( ntest ) = ulpinv
1688  result( ntest+1 ) = ulpinv
1689  result( ntest+2 ) = ulpinv
1690  go to 800
1691  END IF
1692  END IF
1693 *
1694 * Do tests 37 and 38 (or +54)
1695 *
1696  CALL dsyt21( 1, uplo, n, 0, a, lda, d1, d2, z, ldu, v,
1697  $ ldu, tau, work, result( ntest ) )
1698 *
1699  IF( iuplo.EQ.1 ) THEN
1700  indx = 1
1701  DO 760 j = 1, n
1702  DO 750 i = 1, j
1703  work( indx ) = a( i, j )
1704  indx = indx + 1
1705  750 continue
1706  760 continue
1707  ELSE
1708  indx = 1
1709  DO 780 j = 1, n
1710  DO 770 i = j, n
1711  work( indx ) = a( i, j )
1712  indx = indx + 1
1713  770 continue
1714  780 continue
1715  END IF
1716 *
1717  ntest = ntest + 2
1718  srnamt = 'DSPEV'
1719  CALL dspev( 'N', uplo, n, work, d3, z, ldu, v, iinfo )
1720  IF( iinfo.NE.0 ) THEN
1721  WRITE( nounit, fmt = 9999 )'DSPEV(N,' // uplo // ')',
1722  $ iinfo, n, jtype, ioldsd
1723  info = abs( iinfo )
1724  IF( iinfo.LT.0 ) THEN
1725  return
1726  ELSE
1727  result( ntest ) = ulpinv
1728  go to 800
1729  END IF
1730  END IF
1731 *
1732 * Do test 39 (or +54)
1733 *
1734  temp1 = zero
1735  temp2 = zero
1736  DO 790 j = 1, n
1737  temp1 = max( temp1, abs( d1( j ) ), abs( d3( j ) ) )
1738  temp2 = max( temp2, abs( d1( j )-d3( j ) ) )
1739  790 continue
1740  result( ntest ) = temp2 / max( unfl,
1741  $ ulp*max( temp1, temp2 ) )
1742 *
1743 * Load array WORK with the upper or lower triangular part
1744 * of the matrix in packed form.
1745 *
1746  800 continue
1747  IF( iuplo.EQ.1 ) THEN
1748  indx = 1
1749  DO 820 j = 1, n
1750  DO 810 i = 1, j
1751  work( indx ) = a( i, j )
1752  indx = indx + 1
1753  810 continue
1754  820 continue
1755  ELSE
1756  indx = 1
1757  DO 840 j = 1, n
1758  DO 830 i = j, n
1759  work( indx ) = a( i, j )
1760  indx = indx + 1
1761  830 continue
1762  840 continue
1763  END IF
1764 *
1765  ntest = ntest + 1
1766 *
1767  IF( n.GT.0 ) THEN
1768  temp3 = max( abs( d1( 1 ) ), abs( d1( n ) ) )
1769  IF( il.NE.1 ) THEN
1770  vl = d1( il ) - max( half*( d1( il )-d1( il-1 ) ),
1771  $ ten*ulp*temp3, ten*rtunfl )
1772  ELSE IF( n.GT.0 ) THEN
1773  vl = d1( 1 ) - max( half*( d1( n )-d1( 1 ) ),
1774  $ ten*ulp*temp3, ten*rtunfl )
1775  END IF
1776  IF( iu.NE.n ) THEN
1777  vu = d1( iu ) + max( half*( d1( iu+1 )-d1( iu ) ),
1778  $ ten*ulp*temp3, ten*rtunfl )
1779  ELSE IF( n.GT.0 ) THEN
1780  vu = d1( n ) + max( half*( d1( n )-d1( 1 ) ),
1781  $ ten*ulp*temp3, ten*rtunfl )
1782  END IF
1783  ELSE
1784  temp3 = zero
1785  vl = zero
1786  vu = one
1787  END IF
1788 *
1789  srnamt = 'DSPEVX'
1790  CALL dspevx( 'V', 'A', uplo, n, work, vl, vu, il, iu,
1791  $ abstol, m, wa1, z, ldu, v, iwork,
1792  $ iwork( 5*n+1 ), iinfo )
1793  IF( iinfo.NE.0 ) THEN
1794  WRITE( nounit, fmt = 9999 )'DSPEVX(V,A,' // uplo //
1795  $ ')', iinfo, n, jtype, ioldsd
1796  info = abs( iinfo )
1797  IF( iinfo.LT.0 ) THEN
1798  return
1799  ELSE
1800  result( ntest ) = ulpinv
1801  result( ntest+1 ) = ulpinv
1802  result( ntest+2 ) = ulpinv
1803  go to 900
1804  END IF
1805  END IF
1806 *
1807 * Do tests 40 and 41 (or +54)
1808 *
1809  CALL dsyt21( 1, uplo, n, 0, a, ldu, wa1, d2, z, ldu, v,
1810  $ ldu, tau, work, result( ntest ) )
1811 *
1812  ntest = ntest + 2
1813 *
1814  IF( iuplo.EQ.1 ) THEN
1815  indx = 1
1816  DO 860 j = 1, n
1817  DO 850 i = 1, j
1818  work( indx ) = a( i, j )
1819  indx = indx + 1
1820  850 continue
1821  860 continue
1822  ELSE
1823  indx = 1
1824  DO 880 j = 1, n
1825  DO 870 i = j, n
1826  work( indx ) = a( i, j )
1827  indx = indx + 1
1828  870 continue
1829  880 continue
1830  END IF
1831 *
1832  srnamt = 'DSPEVX'
1833  CALL dspevx( 'N', 'A', uplo, n, work, vl, vu, il, iu,
1834  $ abstol, m2, wa2, z, ldu, v, iwork,
1835  $ iwork( 5*n+1 ), iinfo )
1836  IF( iinfo.NE.0 ) THEN
1837  WRITE( nounit, fmt = 9999 )'DSPEVX(N,A,' // uplo //
1838  $ ')', iinfo, n, jtype, ioldsd
1839  info = abs( iinfo )
1840  IF( iinfo.LT.0 ) THEN
1841  return
1842  ELSE
1843  result( ntest ) = ulpinv
1844  go to 900
1845  END IF
1846  END IF
1847 *
1848 * Do test 42 (or +54)
1849 *
1850  temp1 = zero
1851  temp2 = zero
1852  DO 890 j = 1, n
1853  temp1 = max( temp1, abs( wa1( j ) ), abs( wa2( j ) ) )
1854  temp2 = max( temp2, abs( wa1( j )-wa2( j ) ) )
1855  890 continue
1856  result( ntest ) = temp2 / max( unfl,
1857  $ ulp*max( temp1, temp2 ) )
1858 *
1859  900 continue
1860  IF( iuplo.EQ.1 ) THEN
1861  indx = 1
1862  DO 920 j = 1, n
1863  DO 910 i = 1, j
1864  work( indx ) = a( i, j )
1865  indx = indx + 1
1866  910 continue
1867  920 continue
1868  ELSE
1869  indx = 1
1870  DO 940 j = 1, n
1871  DO 930 i = j, n
1872  work( indx ) = a( i, j )
1873  indx = indx + 1
1874  930 continue
1875  940 continue
1876  END IF
1877 *
1878  ntest = ntest + 1
1879 *
1880  srnamt = 'DSPEVX'
1881  CALL dspevx( 'V', 'I', uplo, n, work, vl, vu, il, iu,
1882  $ abstol, m2, wa2, z, ldu, v, iwork,
1883  $ iwork( 5*n+1 ), iinfo )
1884  IF( iinfo.NE.0 ) THEN
1885  WRITE( nounit, fmt = 9999 )'DSPEVX(V,I,' // uplo //
1886  $ ')', iinfo, n, jtype, ioldsd
1887  info = abs( iinfo )
1888  IF( iinfo.LT.0 ) THEN
1889  return
1890  ELSE
1891  result( ntest ) = ulpinv
1892  result( ntest+1 ) = ulpinv
1893  result( ntest+2 ) = ulpinv
1894  go to 990
1895  END IF
1896  END IF
1897 *
1898 * Do tests 43 and 44 (or +54)
1899 *
1900  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
1901  $ v, ldu, tau, work, result( ntest ) )
1902 *
1903  ntest = ntest + 2
1904 *
1905  IF( iuplo.EQ.1 ) THEN
1906  indx = 1
1907  DO 960 j = 1, n
1908  DO 950 i = 1, j
1909  work( indx ) = a( i, j )
1910  indx = indx + 1
1911  950 continue
1912  960 continue
1913  ELSE
1914  indx = 1
1915  DO 980 j = 1, n
1916  DO 970 i = j, n
1917  work( indx ) = a( i, j )
1918  indx = indx + 1
1919  970 continue
1920  980 continue
1921  END IF
1922 *
1923  srnamt = 'DSPEVX'
1924  CALL dspevx( 'N', 'I', uplo, n, work, vl, vu, il, iu,
1925  $ abstol, m3, wa3, z, ldu, v, iwork,
1926  $ iwork( 5*n+1 ), iinfo )
1927  IF( iinfo.NE.0 ) THEN
1928  WRITE( nounit, fmt = 9999 )'DSPEVX(N,I,' // uplo //
1929  $ ')', iinfo, n, jtype, ioldsd
1930  info = abs( iinfo )
1931  IF( iinfo.LT.0 ) THEN
1932  return
1933  ELSE
1934  result( ntest ) = ulpinv
1935  go to 990
1936  END IF
1937  END IF
1938 *
1939  IF( m3.EQ.0 .AND. n.GT.0 ) THEN
1940  result( ntest ) = ulpinv
1941  go to 990
1942  END IF
1943 *
1944 * Do test 45 (or +54)
1945 *
1946  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1947  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1948  IF( n.GT.0 ) THEN
1949  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
1950  ELSE
1951  temp3 = zero
1952  END IF
1953  result( ntest ) = ( temp1+temp2 ) /
1954  $ max( unfl, temp3*ulp )
1955 *
1956  990 continue
1957  IF( iuplo.EQ.1 ) THEN
1958  indx = 1
1959  DO 1010 j = 1, n
1960  DO 1000 i = 1, j
1961  work( indx ) = a( i, j )
1962  indx = indx + 1
1963  1000 continue
1964  1010 continue
1965  ELSE
1966  indx = 1
1967  DO 1030 j = 1, n
1968  DO 1020 i = j, n
1969  work( indx ) = a( i, j )
1970  indx = indx + 1
1971  1020 continue
1972  1030 continue
1973  END IF
1974 *
1975  ntest = ntest + 1
1976 *
1977  srnamt = 'DSPEVX'
1978  CALL dspevx( 'V', 'V', uplo, n, work, vl, vu, il, iu,
1979  $ abstol, m2, wa2, z, ldu, v, iwork,
1980  $ iwork( 5*n+1 ), iinfo )
1981  IF( iinfo.NE.0 ) THEN
1982  WRITE( nounit, fmt = 9999 )'DSPEVX(V,V,' // uplo //
1983  $ ')', iinfo, n, jtype, ioldsd
1984  info = abs( iinfo )
1985  IF( iinfo.LT.0 ) THEN
1986  return
1987  ELSE
1988  result( ntest ) = ulpinv
1989  result( ntest+1 ) = ulpinv
1990  result( ntest+2 ) = ulpinv
1991  go to 1080
1992  END IF
1993  END IF
1994 *
1995 * Do tests 46 and 47 (or +54)
1996 *
1997  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
1998  $ v, ldu, tau, work, result( ntest ) )
1999 *
2000  ntest = ntest + 2
2001 *
2002  IF( iuplo.EQ.1 ) THEN
2003  indx = 1
2004  DO 1050 j = 1, n
2005  DO 1040 i = 1, j
2006  work( indx ) = a( i, j )
2007  indx = indx + 1
2008  1040 continue
2009  1050 continue
2010  ELSE
2011  indx = 1
2012  DO 1070 j = 1, n
2013  DO 1060 i = j, n
2014  work( indx ) = a( i, j )
2015  indx = indx + 1
2016  1060 continue
2017  1070 continue
2018  END IF
2019 *
2020  srnamt = 'DSPEVX'
2021  CALL dspevx( 'N', 'V', uplo, n, work, vl, vu, il, iu,
2022  $ abstol, m3, wa3, z, ldu, v, iwork,
2023  $ iwork( 5*n+1 ), iinfo )
2024  IF( iinfo.NE.0 ) THEN
2025  WRITE( nounit, fmt = 9999 )'DSPEVX(N,V,' // uplo //
2026  $ ')', iinfo, n, jtype, ioldsd
2027  info = abs( iinfo )
2028  IF( iinfo.LT.0 ) THEN
2029  return
2030  ELSE
2031  result( ntest ) = ulpinv
2032  go to 1080
2033  END IF
2034  END IF
2035 *
2036  IF( m3.EQ.0 .AND. n.GT.0 ) THEN
2037  result( ntest ) = ulpinv
2038  go to 1080
2039  END IF
2040 *
2041 * Do test 48 (or +54)
2042 *
2043  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
2044  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
2045  IF( n.GT.0 ) THEN
2046  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
2047  ELSE
2048  temp3 = zero
2049  END IF
2050  result( ntest ) = ( temp1+temp2 ) /
2051  $ max( unfl, temp3*ulp )
2052 *
2053  1080 continue
2054 *
2055 * 6) Call DSBEV and DSBEVX.
2056 *
2057  IF( jtype.LE.7 ) THEN
2058  kd = 1
2059  ELSE IF( jtype.GE.8 .AND. jtype.LE.15 ) THEN
2060  kd = max( n-1, 0 )
2061  ELSE
2062  kd = ihbw
2063  END IF
2064 *
2065 * Load array V with the upper or lower triangular part
2066 * of the matrix in band form.
2067 *
2068  IF( iuplo.EQ.1 ) THEN
2069  DO 1100 j = 1, n
2070  DO 1090 i = max( 1, j-kd ), j
2071  v( kd+1+i-j, j ) = a( i, j )
2072  1090 continue
2073  1100 continue
2074  ELSE
2075  DO 1120 j = 1, n
2076  DO 1110 i = j, min( n, j+kd )
2077  v( 1+i-j, j ) = a( i, j )
2078  1110 continue
2079  1120 continue
2080  END IF
2081 *
2082  ntest = ntest + 1
2083  srnamt = 'DSBEV'
2084  CALL dsbev( 'V', uplo, n, kd, v, ldu, d1, z, ldu, work,
2085  $ iinfo )
2086  IF( iinfo.NE.0 ) THEN
2087  WRITE( nounit, fmt = 9999 )'DSBEV(V,' // uplo // ')',
2088  $ iinfo, n, jtype, ioldsd
2089  info = abs( iinfo )
2090  IF( iinfo.LT.0 ) THEN
2091  return
2092  ELSE
2093  result( ntest ) = ulpinv
2094  result( ntest+1 ) = ulpinv
2095  result( ntest+2 ) = ulpinv
2096  go to 1180
2097  END IF
2098  END IF
2099 *
2100 * Do tests 49 and 50 (or ... )
2101 *
2102  CALL dsyt21( 1, uplo, n, 0, a, lda, d1, d2, z, ldu, v,
2103  $ ldu, tau, work, result( ntest ) )
2104 *
2105  IF( iuplo.EQ.1 ) THEN
2106  DO 1140 j = 1, n
2107  DO 1130 i = max( 1, j-kd ), j
2108  v( kd+1+i-j, j ) = a( i, j )
2109  1130 continue
2110  1140 continue
2111  ELSE
2112  DO 1160 j = 1, n
2113  DO 1150 i = j, min( n, j+kd )
2114  v( 1+i-j, j ) = a( i, j )
2115  1150 continue
2116  1160 continue
2117  END IF
2118 *
2119  ntest = ntest + 2
2120  srnamt = 'DSBEV'
2121  CALL dsbev( 'N', uplo, n, kd, v, ldu, d3, z, ldu, work,
2122  $ iinfo )
2123  IF( iinfo.NE.0 ) THEN
2124  WRITE( nounit, fmt = 9999 )'DSBEV(N,' // uplo // ')',
2125  $ iinfo, n, jtype, ioldsd
2126  info = abs( iinfo )
2127  IF( iinfo.LT.0 ) THEN
2128  return
2129  ELSE
2130  result( ntest ) = ulpinv
2131  go to 1180
2132  END IF
2133  END IF
2134 *
2135 * Do test 51 (or +54)
2136 *
2137  temp1 = zero
2138  temp2 = zero
2139  DO 1170 j = 1, n
2140  temp1 = max( temp1, abs( d1( j ) ), abs( d3( j ) ) )
2141  temp2 = max( temp2, abs( d1( j )-d3( j ) ) )
2142  1170 continue
2143  result( ntest ) = temp2 / max( unfl,
2144  $ ulp*max( temp1, temp2 ) )
2145 *
2146 * Load array V with the upper or lower triangular part
2147 * of the matrix in band form.
2148 *
2149  1180 continue
2150  IF( iuplo.EQ.1 ) THEN
2151  DO 1200 j = 1, n
2152  DO 1190 i = max( 1, j-kd ), j
2153  v( kd+1+i-j, j ) = a( i, j )
2154  1190 continue
2155  1200 continue
2156  ELSE
2157  DO 1220 j = 1, n
2158  DO 1210 i = j, min( n, j+kd )
2159  v( 1+i-j, j ) = a( i, j )
2160  1210 continue
2161  1220 continue
2162  END IF
2163 *
2164  ntest = ntest + 1
2165  srnamt = 'DSBEVX'
2166  CALL dsbevx( 'V', 'A', uplo, n, kd, v, ldu, u, ldu, vl,
2167  $ vu, il, iu, abstol, m, wa2, z, ldu, work,
2168  $ iwork, iwork( 5*n+1 ), iinfo )
2169  IF( iinfo.NE.0 ) THEN
2170  WRITE( nounit, fmt = 9999 )'DSBEVX(V,A,' // uplo //
2171  $ ')', iinfo, n, jtype, ioldsd
2172  info = abs( iinfo )
2173  IF( iinfo.LT.0 ) THEN
2174  return
2175  ELSE
2176  result( ntest ) = ulpinv
2177  result( ntest+1 ) = ulpinv
2178  result( ntest+2 ) = ulpinv
2179  go to 1280
2180  END IF
2181  END IF
2182 *
2183 * Do tests 52 and 53 (or +54)
2184 *
2185  CALL dsyt21( 1, uplo, n, 0, a, ldu, wa2, d2, z, ldu, v,
2186  $ ldu, tau, work, result( ntest ) )
2187 *
2188  ntest = ntest + 2
2189 *
2190  IF( iuplo.EQ.1 ) THEN
2191  DO 1240 j = 1, n
2192  DO 1230 i = max( 1, j-kd ), j
2193  v( kd+1+i-j, j ) = a( i, j )
2194  1230 continue
2195  1240 continue
2196  ELSE
2197  DO 1260 j = 1, n
2198  DO 1250 i = j, min( n, j+kd )
2199  v( 1+i-j, j ) = a( i, j )
2200  1250 continue
2201  1260 continue
2202  END IF
2203 *
2204  srnamt = 'DSBEVX'
2205  CALL dsbevx( 'N', 'A', uplo, n, kd, v, ldu, u, ldu, vl,
2206  $ vu, il, iu, abstol, m3, wa3, z, ldu, work,
2207  $ iwork, iwork( 5*n+1 ), iinfo )
2208  IF( iinfo.NE.0 ) THEN
2209  WRITE( nounit, fmt = 9999 )'DSBEVX(N,A,' // uplo //
2210  $ ')', iinfo, n, jtype, ioldsd
2211  info = abs( iinfo )
2212  IF( iinfo.LT.0 ) THEN
2213  return
2214  ELSE
2215  result( ntest ) = ulpinv
2216  go to 1280
2217  END IF
2218  END IF
2219 *
2220 * Do test 54 (or +54)
2221 *
2222  temp1 = zero
2223  temp2 = zero
2224  DO 1270 j = 1, n
2225  temp1 = max( temp1, abs( wa2( j ) ), abs( wa3( j ) ) )
2226  temp2 = max( temp2, abs( wa2( j )-wa3( j ) ) )
2227  1270 continue
2228  result( ntest ) = temp2 / max( unfl,
2229  $ ulp*max( temp1, temp2 ) )
2230 *
2231  1280 continue
2232  ntest = ntest + 1
2233  IF( iuplo.EQ.1 ) THEN
2234  DO 1300 j = 1, n
2235  DO 1290 i = max( 1, j-kd ), j
2236  v( kd+1+i-j, j ) = a( i, j )
2237  1290 continue
2238  1300 continue
2239  ELSE
2240  DO 1320 j = 1, n
2241  DO 1310 i = j, min( n, j+kd )
2242  v( 1+i-j, j ) = a( i, j )
2243  1310 continue
2244  1320 continue
2245  END IF
2246 *
2247  srnamt = 'DSBEVX'
2248  CALL dsbevx( 'V', 'I', uplo, n, kd, v, ldu, u, ldu, vl,
2249  $ vu, il, iu, abstol, m2, wa2, z, ldu, work,
2250  $ iwork, iwork( 5*n+1 ), iinfo )
2251  IF( iinfo.NE.0 ) THEN
2252  WRITE( nounit, fmt = 9999 )'DSBEVX(V,I,' // uplo //
2253  $ ')', iinfo, n, jtype, ioldsd
2254  info = abs( iinfo )
2255  IF( iinfo.LT.0 ) THEN
2256  return
2257  ELSE
2258  result( ntest ) = ulpinv
2259  result( ntest+1 ) = ulpinv
2260  result( ntest+2 ) = ulpinv
2261  go to 1370
2262  END IF
2263  END IF
2264 *
2265 * Do tests 55 and 56 (or +54)
2266 *
2267  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
2268  $ v, ldu, tau, work, result( ntest ) )
2269 *
2270  ntest = ntest + 2
2271 *
2272  IF( iuplo.EQ.1 ) THEN
2273  DO 1340 j = 1, n
2274  DO 1330 i = max( 1, j-kd ), j
2275  v( kd+1+i-j, j ) = a( i, j )
2276  1330 continue
2277  1340 continue
2278  ELSE
2279  DO 1360 j = 1, n
2280  DO 1350 i = j, min( n, j+kd )
2281  v( 1+i-j, j ) = a( i, j )
2282  1350 continue
2283  1360 continue
2284  END IF
2285 *
2286  srnamt = 'DSBEVX'
2287  CALL dsbevx( 'N', 'I', uplo, n, kd, v, ldu, u, ldu, vl,
2288  $ vu, il, iu, abstol, m3, wa3, z, ldu, work,
2289  $ iwork, iwork( 5*n+1 ), iinfo )
2290  IF( iinfo.NE.0 ) THEN
2291  WRITE( nounit, fmt = 9999 )'DSBEVX(N,I,' // uplo //
2292  $ ')', iinfo, n, jtype, ioldsd
2293  info = abs( iinfo )
2294  IF( iinfo.LT.0 ) THEN
2295  return
2296  ELSE
2297  result( ntest ) = ulpinv
2298  go to 1370
2299  END IF
2300  END IF
2301 *
2302 * Do test 57 (or +54)
2303 *
2304  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
2305  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
2306  IF( n.GT.0 ) THEN
2307  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
2308  ELSE
2309  temp3 = zero
2310  END IF
2311  result( ntest ) = ( temp1+temp2 ) /
2312  $ max( unfl, temp3*ulp )
2313 *
2314  1370 continue
2315  ntest = ntest + 1
2316  IF( iuplo.EQ.1 ) THEN
2317  DO 1390 j = 1, n
2318  DO 1380 i = max( 1, j-kd ), j
2319  v( kd+1+i-j, j ) = a( i, j )
2320  1380 continue
2321  1390 continue
2322  ELSE
2323  DO 1410 j = 1, n
2324  DO 1400 i = j, min( n, j+kd )
2325  v( 1+i-j, j ) = a( i, j )
2326  1400 continue
2327  1410 continue
2328  END IF
2329 *
2330  srnamt = 'DSBEVX'
2331  CALL dsbevx( 'V', 'V', uplo, n, kd, v, ldu, u, ldu, vl,
2332  $ vu, il, iu, abstol, m2, wa2, z, ldu, work,
2333  $ iwork, iwork( 5*n+1 ), iinfo )
2334  IF( iinfo.NE.0 ) THEN
2335  WRITE( nounit, fmt = 9999 )'DSBEVX(V,V,' // uplo //
2336  $ ')', iinfo, n, jtype, ioldsd
2337  info = abs( iinfo )
2338  IF( iinfo.LT.0 ) THEN
2339  return
2340  ELSE
2341  result( ntest ) = ulpinv
2342  result( ntest+1 ) = ulpinv
2343  result( ntest+2 ) = ulpinv
2344  go to 1460
2345  END IF
2346  END IF
2347 *
2348 * Do tests 58 and 59 (or +54)
2349 *
2350  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
2351  $ v, ldu, tau, work, result( ntest ) )
2352 *
2353  ntest = ntest + 2
2354 *
2355  IF( iuplo.EQ.1 ) THEN
2356  DO 1430 j = 1, n
2357  DO 1420 i = max( 1, j-kd ), j
2358  v( kd+1+i-j, j ) = a( i, j )
2359  1420 continue
2360  1430 continue
2361  ELSE
2362  DO 1450 j = 1, n
2363  DO 1440 i = j, min( n, j+kd )
2364  v( 1+i-j, j ) = a( i, j )
2365  1440 continue
2366  1450 continue
2367  END IF
2368 *
2369  srnamt = 'DSBEVX'
2370  CALL dsbevx( 'N', 'V', uplo, n, kd, v, ldu, u, ldu, vl,
2371  $ vu, il, iu, abstol, m3, wa3, z, ldu, work,
2372  $ iwork, iwork( 5*n+1 ), iinfo )
2373  IF( iinfo.NE.0 ) THEN
2374  WRITE( nounit, fmt = 9999 )'DSBEVX(N,V,' // uplo //
2375  $ ')', iinfo, n, jtype, ioldsd
2376  info = abs( iinfo )
2377  IF( iinfo.LT.0 ) THEN
2378  return
2379  ELSE
2380  result( ntest ) = ulpinv
2381  go to 1460
2382  END IF
2383  END IF
2384 *
2385  IF( m3.EQ.0 .AND. n.GT.0 ) THEN
2386  result( ntest ) = ulpinv
2387  go to 1460
2388  END IF
2389 *
2390 * Do test 60 (or +54)
2391 *
2392  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
2393  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
2394  IF( n.GT.0 ) THEN
2395  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
2396  ELSE
2397  temp3 = zero
2398  END IF
2399  result( ntest ) = ( temp1+temp2 ) /
2400  $ max( unfl, temp3*ulp )
2401 *
2402  1460 continue
2403 *
2404 * 7) Call DSYEVD
2405 *
2406  CALL dlacpy( ' ', n, n, a, lda, v, ldu )
2407 *
2408  ntest = ntest + 1
2409  srnamt = 'DSYEVD'
2410  CALL dsyevd( 'V', uplo, n, a, ldu, d1, work, lwedc,
2411  $ iwork, liwedc, iinfo )
2412  IF( iinfo.NE.0 ) THEN
2413  WRITE( nounit, fmt = 9999 )'DSYEVD(V,' // uplo //
2414  $ ')', iinfo, n, jtype, ioldsd
2415  info = abs( iinfo )
2416  IF( iinfo.LT.0 ) THEN
2417  return
2418  ELSE
2419  result( ntest ) = ulpinv
2420  result( ntest+1 ) = ulpinv
2421  result( ntest+2 ) = ulpinv
2422  go to 1480
2423  END IF
2424  END IF
2425 *
2426 * Do tests 61 and 62 (or +54)
2427 *
2428  CALL dsyt21( 1, uplo, n, 0, v, ldu, d1, d2, a, ldu, z,
2429  $ ldu, tau, work, result( ntest ) )
2430 *
2431  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2432 *
2433  ntest = ntest + 2
2434  srnamt = 'DSYEVD'
2435  CALL dsyevd( 'N', uplo, n, a, ldu, d3, work, lwedc,
2436  $ iwork, liwedc, iinfo )
2437  IF( iinfo.NE.0 ) THEN
2438  WRITE( nounit, fmt = 9999 )'DSYEVD(N,' // uplo //
2439  $ ')', iinfo, n, jtype, ioldsd
2440  info = abs( iinfo )
2441  IF( iinfo.LT.0 ) THEN
2442  return
2443  ELSE
2444  result( ntest ) = ulpinv
2445  go to 1480
2446  END IF
2447  END IF
2448 *
2449 * Do test 63 (or +54)
2450 *
2451  temp1 = zero
2452  temp2 = zero
2453  DO 1470 j = 1, n
2454  temp1 = max( temp1, abs( d1( j ) ), abs( d3( j ) ) )
2455  temp2 = max( temp2, abs( d1( j )-d3( j ) ) )
2456  1470 continue
2457  result( ntest ) = temp2 / max( unfl,
2458  $ ulp*max( temp1, temp2 ) )
2459 *
2460  1480 continue
2461 *
2462 * 8) Call DSPEVD.
2463 *
2464  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2465 *
2466 * Load array WORK with the upper or lower triangular
2467 * part of the matrix in packed form.
2468 *
2469  IF( iuplo.EQ.1 ) THEN
2470  indx = 1
2471  DO 1500 j = 1, n
2472  DO 1490 i = 1, j
2473  work( indx ) = a( i, j )
2474  indx = indx + 1
2475  1490 continue
2476  1500 continue
2477  ELSE
2478  indx = 1
2479  DO 1520 j = 1, n
2480  DO 1510 i = j, n
2481  work( indx ) = a( i, j )
2482  indx = indx + 1
2483  1510 continue
2484  1520 continue
2485  END IF
2486 *
2487  ntest = ntest + 1
2488  srnamt = 'DSPEVD'
2489  CALL dspevd( 'V', uplo, n, work, d1, z, ldu,
2490  $ work( indx ), lwedc-indx+1, iwork, liwedc,
2491  $ iinfo )
2492  IF( iinfo.NE.0 ) THEN
2493  WRITE( nounit, fmt = 9999 )'DSPEVD(V,' // uplo //
2494  $ ')', iinfo, n, jtype, ioldsd
2495  info = abs( iinfo )
2496  IF( iinfo.LT.0 ) THEN
2497  return
2498  ELSE
2499  result( ntest ) = ulpinv
2500  result( ntest+1 ) = ulpinv
2501  result( ntest+2 ) = ulpinv
2502  go to 1580
2503  END IF
2504  END IF
2505 *
2506 * Do tests 64 and 65 (or +54)
2507 *
2508  CALL dsyt21( 1, uplo, n, 0, a, lda, d1, d2, z, ldu, v,
2509  $ ldu, tau, work, result( ntest ) )
2510 *
2511  IF( iuplo.EQ.1 ) THEN
2512  indx = 1
2513  DO 1540 j = 1, n
2514  DO 1530 i = 1, j
2515 *
2516  work( indx ) = a( i, j )
2517  indx = indx + 1
2518  1530 continue
2519  1540 continue
2520  ELSE
2521  indx = 1
2522  DO 1560 j = 1, n
2523  DO 1550 i = j, n
2524  work( indx ) = a( i, j )
2525  indx = indx + 1
2526  1550 continue
2527  1560 continue
2528  END IF
2529 *
2530  ntest = ntest + 2
2531  srnamt = 'DSPEVD'
2532  CALL dspevd( 'N', uplo, n, work, d3, z, ldu,
2533  $ work( indx ), lwedc-indx+1, iwork, liwedc,
2534  $ iinfo )
2535  IF( iinfo.NE.0 ) THEN
2536  WRITE( nounit, fmt = 9999 )'DSPEVD(N,' // uplo //
2537  $ ')', iinfo, n, jtype, ioldsd
2538  info = abs( iinfo )
2539  IF( iinfo.LT.0 ) THEN
2540  return
2541  ELSE
2542  result( ntest ) = ulpinv
2543  go to 1580
2544  END IF
2545  END IF
2546 *
2547 * Do test 66 (or +54)
2548 *
2549  temp1 = zero
2550  temp2 = zero
2551  DO 1570 j = 1, n
2552  temp1 = max( temp1, abs( d1( j ) ), abs( d3( j ) ) )
2553  temp2 = max( temp2, abs( d1( j )-d3( j ) ) )
2554  1570 continue
2555  result( ntest ) = temp2 / max( unfl,
2556  $ ulp*max( temp1, temp2 ) )
2557  1580 continue
2558 *
2559 * 9) Call DSBEVD.
2560 *
2561  IF( jtype.LE.7 ) THEN
2562  kd = 1
2563  ELSE IF( jtype.GE.8 .AND. jtype.LE.15 ) THEN
2564  kd = max( n-1, 0 )
2565  ELSE
2566  kd = ihbw
2567  END IF
2568 *
2569 * Load array V with the upper or lower triangular part
2570 * of the matrix in band form.
2571 *
2572  IF( iuplo.EQ.1 ) THEN
2573  DO 1600 j = 1, n
2574  DO 1590 i = max( 1, j-kd ), j
2575  v( kd+1+i-j, j ) = a( i, j )
2576  1590 continue
2577  1600 continue
2578  ELSE
2579  DO 1620 j = 1, n
2580  DO 1610 i = j, min( n, j+kd )
2581  v( 1+i-j, j ) = a( i, j )
2582  1610 continue
2583  1620 continue
2584  END IF
2585 *
2586  ntest = ntest + 1
2587  srnamt = 'DSBEVD'
2588  CALL dsbevd( 'V', uplo, n, kd, v, ldu, d1, z, ldu, work,
2589  $ lwedc, iwork, liwedc, iinfo )
2590  IF( iinfo.NE.0 ) THEN
2591  WRITE( nounit, fmt = 9999 )'DSBEVD(V,' // uplo //
2592  $ ')', iinfo, n, jtype, ioldsd
2593  info = abs( iinfo )
2594  IF( iinfo.LT.0 ) THEN
2595  return
2596  ELSE
2597  result( ntest ) = ulpinv
2598  result( ntest+1 ) = ulpinv
2599  result( ntest+2 ) = ulpinv
2600  go to 1680
2601  END IF
2602  END IF
2603 *
2604 * Do tests 67 and 68 (or +54)
2605 *
2606  CALL dsyt21( 1, uplo, n, 0, a, lda, d1, d2, z, ldu, v,
2607  $ ldu, tau, work, result( ntest ) )
2608 *
2609  IF( iuplo.EQ.1 ) THEN
2610  DO 1640 j = 1, n
2611  DO 1630 i = max( 1, j-kd ), j
2612  v( kd+1+i-j, j ) = a( i, j )
2613  1630 continue
2614  1640 continue
2615  ELSE
2616  DO 1660 j = 1, n
2617  DO 1650 i = j, min( n, j+kd )
2618  v( 1+i-j, j ) = a( i, j )
2619  1650 continue
2620  1660 continue
2621  END IF
2622 *
2623  ntest = ntest + 2
2624  srnamt = 'DSBEVD'
2625  CALL dsbevd( 'N', uplo, n, kd, v, ldu, d3, z, ldu, work,
2626  $ lwedc, iwork, liwedc, iinfo )
2627  IF( iinfo.NE.0 ) THEN
2628  WRITE( nounit, fmt = 9999 )'DSBEVD(N,' // uplo //
2629  $ ')', iinfo, n, jtype, ioldsd
2630  info = abs( iinfo )
2631  IF( iinfo.LT.0 ) THEN
2632  return
2633  ELSE
2634  result( ntest ) = ulpinv
2635  go to 1680
2636  END IF
2637  END IF
2638 *
2639 * Do test 69 (or +54)
2640 *
2641  temp1 = zero
2642  temp2 = zero
2643  DO 1670 j = 1, n
2644  temp1 = max( temp1, abs( d1( j ) ), abs( d3( j ) ) )
2645  temp2 = max( temp2, abs( d1( j )-d3( j ) ) )
2646  1670 continue
2647  result( ntest ) = temp2 / max( unfl,
2648  $ ulp*max( temp1, temp2 ) )
2649 *
2650  1680 continue
2651 *
2652 *
2653  CALL dlacpy( ' ', n, n, a, lda, v, ldu )
2654  ntest = ntest + 1
2655  srnamt = 'DSYEVR'
2656  CALL dsyevr( 'V', 'A', uplo, n, a, ldu, vl, vu, il, iu,
2657  $ abstol, m, wa1, z, ldu, iwork, work, lwork,
2658  $ iwork(2*n+1), liwork-2*n, iinfo )
2659  IF( iinfo.NE.0 ) THEN
2660  WRITE( nounit, fmt = 9999 )'DSYEVR(V,A,' // uplo //
2661  $ ')', iinfo, n, jtype, ioldsd
2662  info = abs( iinfo )
2663  IF( iinfo.LT.0 ) THEN
2664  return
2665  ELSE
2666  result( ntest ) = ulpinv
2667  result( ntest+1 ) = ulpinv
2668  result( ntest+2 ) = ulpinv
2669  go to 1700
2670  END IF
2671  END IF
2672 *
2673 * Do tests 70 and 71 (or ... )
2674 *
2675  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2676 *
2677  CALL dsyt21( 1, uplo, n, 0, a, ldu, wa1, d2, z, ldu, v,
2678  $ ldu, tau, work, result( ntest ) )
2679 *
2680  ntest = ntest + 2
2681  srnamt = 'DSYEVR'
2682  CALL dsyevr( 'N', 'A', uplo, n, a, ldu, vl, vu, il, iu,
2683  $ abstol, m2, wa2, z, ldu, iwork, work, lwork,
2684  $ iwork(2*n+1), liwork-2*n, iinfo )
2685  IF( iinfo.NE.0 ) THEN
2686  WRITE( nounit, fmt = 9999 )'DSYEVR(N,A,' // uplo //
2687  $ ')', iinfo, n, jtype, ioldsd
2688  info = abs( iinfo )
2689  IF( iinfo.LT.0 ) THEN
2690  return
2691  ELSE
2692  result( ntest ) = ulpinv
2693  go to 1700
2694  END IF
2695  END IF
2696 *
2697 * Do test 72 (or ... )
2698 *
2699  temp1 = zero
2700  temp2 = zero
2701  DO 1690 j = 1, n
2702  temp1 = max( temp1, abs( wa1( j ) ), abs( wa2( j ) ) )
2703  temp2 = max( temp2, abs( wa1( j )-wa2( j ) ) )
2704  1690 continue
2705  result( ntest ) = temp2 / max( unfl,
2706  $ ulp*max( temp1, temp2 ) )
2707 *
2708  1700 continue
2709 *
2710  ntest = ntest + 1
2711  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2712  srnamt = 'DSYEVR'
2713  CALL dsyevr( 'V', 'I', uplo, n, a, ldu, vl, vu, il, iu,
2714  $ abstol, m2, wa2, z, ldu, iwork, work, lwork,
2715  $ iwork(2*n+1), liwork-2*n, iinfo )
2716  IF( iinfo.NE.0 ) THEN
2717  WRITE( nounit, fmt = 9999 )'DSYEVR(V,I,' // uplo //
2718  $ ')', iinfo, n, jtype, ioldsd
2719  info = abs( iinfo )
2720  IF( iinfo.LT.0 ) THEN
2721  return
2722  ELSE
2723  result( ntest ) = ulpinv
2724  result( ntest+1 ) = ulpinv
2725  result( ntest+2 ) = ulpinv
2726  go to 1710
2727  END IF
2728  END IF
2729 *
2730 * Do tests 73 and 74 (or +54)
2731 *
2732  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2733 *
2734  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
2735  $ v, ldu, tau, work, result( ntest ) )
2736 *
2737  ntest = ntest + 2
2738  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2739  srnamt = 'DSYEVR'
2740  CALL dsyevr( 'N', 'I', uplo, n, a, ldu, vl, vu, il, iu,
2741  $ abstol, m3, wa3, z, ldu, iwork, work, lwork,
2742  $ iwork(2*n+1), liwork-2*n, iinfo )
2743  IF( iinfo.NE.0 ) THEN
2744  WRITE( nounit, fmt = 9999 )'DSYEVR(N,I,' // uplo //
2745  $ ')', iinfo, n, jtype, ioldsd
2746  info = abs( iinfo )
2747  IF( iinfo.LT.0 ) THEN
2748  return
2749  ELSE
2750  result( ntest ) = ulpinv
2751  go to 1710
2752  END IF
2753  END IF
2754 *
2755 * Do test 75 (or +54)
2756 *
2757  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
2758  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
2759  result( ntest ) = ( temp1+temp2 ) /
2760  $ max( unfl, ulp*temp3 )
2761  1710 continue
2762 *
2763  ntest = ntest + 1
2764  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2765  srnamt = 'DSYEVR'
2766  CALL dsyevr( 'V', 'V', uplo, n, a, ldu, vl, vu, il, iu,
2767  $ abstol, m2, wa2, z, ldu, iwork, work, lwork,
2768  $ iwork(2*n+1), liwork-2*n, iinfo )
2769  IF( iinfo.NE.0 ) THEN
2770  WRITE( nounit, fmt = 9999 )'DSYEVR(V,V,' // uplo //
2771  $ ')', iinfo, n, jtype, ioldsd
2772  info = abs( iinfo )
2773  IF( iinfo.LT.0 ) THEN
2774  return
2775  ELSE
2776  result( ntest ) = ulpinv
2777  result( ntest+1 ) = ulpinv
2778  result( ntest+2 ) = ulpinv
2779  go to 700
2780  END IF
2781  END IF
2782 *
2783 * Do tests 76 and 77 (or +54)
2784 *
2785  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2786 *
2787  CALL dsyt22( 1, uplo, n, m2, 0, a, ldu, wa2, d2, z, ldu,
2788  $ v, ldu, tau, work, result( ntest ) )
2789 *
2790  ntest = ntest + 2
2791  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2792  srnamt = 'DSYEVR'
2793  CALL dsyevr( 'N', 'V', uplo, n, a, ldu, vl, vu, il, iu,
2794  $ abstol, m3, wa3, z, ldu, iwork, work, lwork,
2795  $ iwork(2*n+1), liwork-2*n, iinfo )
2796  IF( iinfo.NE.0 ) THEN
2797  WRITE( nounit, fmt = 9999 )'DSYEVR(N,V,' // uplo //
2798  $ ')', iinfo, n, jtype, ioldsd
2799  info = abs( iinfo )
2800  IF( iinfo.LT.0 ) THEN
2801  return
2802  ELSE
2803  result( ntest ) = ulpinv
2804  go to 700
2805  END IF
2806  END IF
2807 *
2808  IF( m3.EQ.0 .AND. n.GT.0 ) THEN
2809  result( ntest ) = ulpinv
2810  go to 700
2811  END IF
2812 *
2813 * Do test 78 (or +54)
2814 *
2815  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
2816  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
2817  IF( n.GT.0 ) THEN
2818  temp3 = max( abs( wa1( 1 ) ), abs( wa1( n ) ) )
2819  ELSE
2820  temp3 = zero
2821  END IF
2822  result( ntest ) = ( temp1+temp2 ) /
2823  $ max( unfl, temp3*ulp )
2824 *
2825  CALL dlacpy( ' ', n, n, v, ldu, a, lda )
2826 *
2827  1720 continue
2828 *
2829 * End of Loop -- Check for RESULT(j) > THRESH
2830 *
2831  ntestt = ntestt + ntest
2832 *
2833  CALL dlafts( 'DST', n, n, jtype, ntest, result, ioldsd,
2834  $ thresh, nounit, nerrs )
2835 *
2836  1730 continue
2837  1740 continue
2838 *
2839 * Summary
2840 *
2841  CALL alasvm( 'DST', nounit, nerrs, ntestt, 0 )
2842 *
2843  9999 format( ' DDRVST: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
2844  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
2845 *
2846  return
2847 *
2848 * End of DDRVST
2849 *
2850  END