LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
ddrves.f
Go to the documentation of this file.
1 *> \brief \b DDRVES
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 DDRVES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * NOUNIT, A, LDA, H, HT, WR, WI, WRT, WIT, VS,
13 * LDVS, RESULT, WORK, NWORK, IWORK, BWORK, INFO )
14 *
15 * .. Scalar Arguments ..
16 * INTEGER INFO, LDA, LDVS, NOUNIT, NSIZES, NTYPES, NWORK
17 * DOUBLE PRECISION THRESH
18 * ..
19 * .. Array Arguments ..
20 * LOGICAL BWORK( * ), DOTYPE( * )
21 * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
22 * DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
23 * \$ RESULT( 13 ), VS( LDVS, * ), WI( * ), WIT( * ),
24 * \$ WORK( * ), WR( * ), WRT( * )
25 * ..
26 *
27 *
28 *> \par Purpose:
29 * =============
30 *>
31 *> \verbatim
32 *>
33 *> DDRVES checks the nonsymmetric eigenvalue (Schur form) problem
34 *> driver DGEES.
35 *>
36 *> When DDRVES is called, a number of matrix "sizes" ("n's") and a
37 *> number of matrix "types" are specified. For each size ("n")
38 *> and each type of matrix, one matrix will be generated and used
39 *> to test the nonsymmetric eigenroutines. For each matrix, 13
40 *> tests will be performed:
41 *>
42 *> (1) 0 if T is in Schur form, 1/ulp otherwise
43 *> (no sorting of eigenvalues)
44 *>
45 *> (2) | A - VS T VS' | / ( n |A| ulp )
46 *>
47 *> Here VS is the matrix of Schur eigenvectors, and T is in Schur
48 *> form (no sorting of eigenvalues).
49 *>
50 *> (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
51 *>
52 *> (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T
53 *> 1/ulp otherwise
54 *> (no sorting of eigenvalues)
55 *>
56 *> (5) 0 if T(with VS) = T(without VS),
57 *> 1/ulp otherwise
58 *> (no sorting of eigenvalues)
59 *>
60 *> (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
61 *> 1/ulp otherwise
62 *> (no sorting of eigenvalues)
63 *>
64 *> (7) 0 if T is in Schur form, 1/ulp otherwise
65 *> (with sorting of eigenvalues)
66 *>
67 *> (8) | A - VS T VS' | / ( n |A| ulp )
68 *>
69 *> Here VS is the matrix of Schur eigenvectors, and T is in Schur
70 *> form (with sorting of eigenvalues).
71 *>
72 *> (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
73 *>
74 *> (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T
75 *> 1/ulp otherwise
76 *> (with sorting of eigenvalues)
77 *>
78 *> (11) 0 if T(with VS) = T(without VS),
79 *> 1/ulp otherwise
80 *> (with sorting of eigenvalues)
81 *>
82 *> (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
83 *> 1/ulp otherwise
84 *> (with sorting of eigenvalues)
85 *>
86 *> (13) if sorting worked and SDIM is the number of
87 *> eigenvalues which were SELECTed
88 *>
89 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
90 *> each element NN(j) specifies one size.
91 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
92 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
93 *> Currently, the list of possible types is:
94 *>
95 *> (1) The zero matrix.
96 *> (2) The identity matrix.
97 *> (3) A (transposed) Jordan block, with 1's on the diagonal.
98 *>
99 *> (4) A diagonal matrix with evenly spaced entries
100 *> 1, ..., ULP and random signs.
101 *> (ULP = (first number larger than 1) - 1 )
102 *> (5) A diagonal matrix with geometrically spaced entries
103 *> 1, ..., ULP and random signs.
104 *> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
105 *> and random signs.
106 *>
107 *> (7) Same as (4), but multiplied by a constant near
108 *> the overflow threshold
109 *> (8) Same as (4), but multiplied by a constant near
110 *> the underflow threshold
111 *>
112 *> (9) A matrix of the form U' T U, where U is orthogonal and
113 *> T has evenly spaced entries 1, ..., ULP with random signs
114 *> on the diagonal and random O(1) entries in the upper
115 *> triangle.
116 *>
117 *> (10) A matrix of the form U' T U, where U is orthogonal and
118 *> T has geometrically spaced entries 1, ..., ULP with random
119 *> signs on the diagonal and random O(1) entries in the upper
120 *> triangle.
121 *>
122 *> (11) A matrix of the form U' T U, where U is orthogonal and
123 *> T has "clustered" entries 1, ULP,..., ULP with random
124 *> signs on the diagonal and random O(1) entries in the upper
125 *> triangle.
126 *>
127 *> (12) A matrix of the form U' T U, where U is orthogonal and
128 *> T has real or complex conjugate paired eigenvalues randomly
129 *> chosen from ( ULP, 1 ) and random O(1) entries in the upper
130 *> triangle.
131 *>
132 *> (13) A matrix of the form X' T X, where X has condition
133 *> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
134 *> with random signs on the diagonal and random O(1) entries
135 *> in the upper triangle.
136 *>
137 *> (14) A matrix of the form X' T X, where X has condition
138 *> SQRT( ULP ) and T has geometrically spaced entries
139 *> 1, ..., ULP with random signs on the diagonal and random
140 *> O(1) entries in the upper triangle.
141 *>
142 *> (15) A matrix of the form X' T X, where X has condition
143 *> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
144 *> with random signs on the diagonal and random O(1) entries
145 *> in the upper triangle.
146 *>
147 *> (16) A matrix of the form X' T X, where X has condition
148 *> SQRT( ULP ) and T has real or complex conjugate paired
149 *> eigenvalues randomly chosen from ( ULP, 1 ) and random
150 *> O(1) entries in the upper triangle.
151 *>
152 *> (17) Same as (16), but multiplied by a constant
153 *> near the overflow threshold
154 *> (18) Same as (16), but multiplied by a constant
155 *> near the underflow threshold
156 *>
157 *> (19) Nonsymmetric matrix with random entries chosen from (-1,1).
158 *> If N is at least 4, all entries in first two rows and last
159 *> row, and first column and last two columns are zero.
160 *> (20) Same as (19), but multiplied by a constant
161 *> near the overflow threshold
162 *> (21) Same as (19), but multiplied by a constant
163 *> near the underflow threshold
164 *> \endverbatim
165 *
166 * Arguments:
167 * ==========
168 *
169 *> \param[in] NSIZES
170 *> \verbatim
171 *> NSIZES is INTEGER
172 *> The number of sizes of matrices to use. If it is zero,
173 *> DDRVES does nothing. It must be at least zero.
174 *> \endverbatim
175 *>
176 *> \param[in] NN
177 *> \verbatim
178 *> NN is INTEGER array, dimension (NSIZES)
179 *> An array containing the sizes to be used for the matrices.
180 *> Zero values will be skipped. The values must be at least
181 *> zero.
182 *> \endverbatim
183 *>
184 *> \param[in] NTYPES
185 *> \verbatim
186 *> NTYPES is INTEGER
187 *> The number of elements in DOTYPE. If it is zero, DDRVES
188 *> does nothing. It must be at least zero. If it is MAXTYP+1
189 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
190 *> defined, which is to use whatever matrix is in A. This
191 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
192 *> DOTYPE(MAXTYP+1) is .TRUE. .
193 *> \endverbatim
194 *>
195 *> \param[in] DOTYPE
196 *> \verbatim
197 *> DOTYPE is LOGICAL array, dimension (NTYPES)
198 *> If DOTYPE(j) is .TRUE., then for each size in NN a
199 *> matrix of that size and of type j will be generated.
200 *> If NTYPES is smaller than the maximum number of types
201 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
202 *> MAXTYP will not be generated. If NTYPES is larger
203 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
204 *> will be ignored.
205 *> \endverbatim
206 *>
207 *> \param[in,out] ISEED
208 *> \verbatim
209 *> ISEED is INTEGER array, dimension (4)
210 *> On entry ISEED specifies the seed of the random number
211 *> generator. The array elements should be between 0 and 4095;
212 *> if not they will be reduced mod 4096. Also, ISEED(4) must
213 *> be odd. The random number generator uses a linear
214 *> congruential sequence limited to small integers, and so
215 *> should produce machine independent random numbers. The
216 *> values of ISEED are changed on exit, and can be used in the
217 *> next call to DDRVES to continue the same random number
218 *> sequence.
219 *> \endverbatim
220 *>
221 *> \param[in] THRESH
222 *> \verbatim
223 *> THRESH is DOUBLE PRECISION
224 *> A test will count as "failed" if the "error", computed as
225 *> described above, exceeds THRESH. Note that the error
226 *> is scaled to be O(1), so THRESH should be a reasonably
227 *> small multiple of 1, e.g., 10 or 100. In particular,
228 *> it should not depend on the precision (single vs. double)
229 *> or the size of the matrix. It must be at least zero.
230 *> \endverbatim
231 *>
232 *> \param[in] NOUNIT
233 *> \verbatim
234 *> NOUNIT is INTEGER
235 *> The FORTRAN unit number for printing out error messages
236 *> (e.g., if a routine returns INFO not equal to 0.)
237 *> \endverbatim
238 *>
239 *> \param[out] A
240 *> \verbatim
241 *> A is DOUBLE PRECISION array, dimension (LDA, max(NN))
242 *> Used to hold the matrix whose eigenvalues are to be
243 *> computed. On exit, A contains the last matrix actually used.
244 *> \endverbatim
245 *>
246 *> \param[in] LDA
247 *> \verbatim
248 *> LDA is INTEGER
249 *> The leading dimension of A, and H. LDA must be at
250 *> least 1 and at least max(NN).
251 *> \endverbatim
252 *>
253 *> \param[out] H
254 *> \verbatim
255 *> H is DOUBLE PRECISION array, dimension (LDA, max(NN))
256 *> Another copy of the test matrix A, modified by DGEES.
257 *> \endverbatim
258 *>
259 *> \param[out] HT
260 *> \verbatim
261 *> HT is DOUBLE PRECISION array, dimension (LDA, max(NN))
262 *> Yet another copy of the test matrix A, modified by DGEES.
263 *> \endverbatim
264 *>
265 *> \param[out] WR
266 *> \verbatim
267 *> WR is DOUBLE PRECISION array, dimension (max(NN))
268 *> \endverbatim
269 *>
270 *> \param[out] WI
271 *> \verbatim
272 *> WI is DOUBLE PRECISION array, dimension (max(NN))
273 *>
274 *> The real and imaginary parts of the eigenvalues of A.
275 *> On exit, WR + WI*i are the eigenvalues of the matrix in A.
276 *> \endverbatim
277 *>
278 *> \param[out] WRT
279 *> \verbatim
280 *> WRT is DOUBLE PRECISION array, dimension (max(NN))
281 *> \endverbatim
282 *>
283 *> \param[out] WIT
284 *> \verbatim
285 *> WIT is DOUBLE PRECISION array, dimension (max(NN))
286 *>
287 *> Like WR, WI, these arrays contain the eigenvalues of A,
288 *> but those computed when DGEES only computes a partial
289 *> eigendecomposition, i.e. not Schur vectors
290 *> \endverbatim
291 *>
292 *> \param[out] VS
293 *> \verbatim
294 *> VS is DOUBLE PRECISION array, dimension (LDVS, max(NN))
295 *> VS holds the computed Schur vectors.
296 *> \endverbatim
297 *>
298 *> \param[in] LDVS
299 *> \verbatim
300 *> LDVS is INTEGER
301 *> Leading dimension of VS. Must be at least max(1,max(NN)).
302 *> \endverbatim
303 *>
304 *> \param[out] RESULT
305 *> \verbatim
306 *> RESULT is DOUBLE PRECISION array, dimension (13)
307 *> The values computed by the 13 tests described above.
308 *> The values are currently limited to 1/ulp, to avoid overflow.
309 *> \endverbatim
310 *>
311 *> \param[out] WORK
312 *> \verbatim
313 *> WORK is DOUBLE PRECISION array, dimension (NWORK)
314 *> \endverbatim
315 *>
316 *> \param[in] NWORK
317 *> \verbatim
318 *> NWORK is INTEGER
319 *> The number of entries in WORK. This must be at least
320 *> 5*NN(j)+2*NN(j)**2 for all j.
321 *> \endverbatim
322 *>
323 *> \param[out] IWORK
324 *> \verbatim
325 *> IWORK is INTEGER array, dimension (max(NN))
326 *> \endverbatim
327 *>
328 *> \param[out] BWORK
329 *> \verbatim
330 *> BWORK is LOGICAL array, dimension (max(NN))
331 *> \endverbatim
332 *>
333 *> \param[out] INFO
334 *> \verbatim
335 *> INFO is INTEGER
336 *> If 0, then everything ran OK.
337 *> -1: NSIZES < 0
338 *> -2: Some NN(j) < 0
339 *> -3: NTYPES < 0
340 *> -6: THRESH < 0
341 *> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
342 *> -17: LDVS < 1 or LDVS < NMAX, where NMAX is max( NN(j) ).
343 *> -20: NWORK too small.
344 *> If DLATMR, SLATMS, SLATME or DGEES returns an error code,
345 *> the absolute value of it is returned.
346 *>
347 *>-----------------------------------------------------------------------
348 *>
349 *> Some Local Variables and Parameters:
350 *> ---- ----- --------- --- ----------
351 *>
352 *> ZERO, ONE Real 0 and 1.
353 *> MAXTYP The number of types defined.
354 *> NMAX Largest value in NN.
355 *> NERRS The number of tests which have exceeded THRESH
356 *> COND, CONDS,
357 *> IMODE Values to be passed to the matrix generators.
358 *> ANORM Norm of A; passed to matrix generators.
359 *>
360 *> OVFL, UNFL Overflow and underflow thresholds.
361 *> ULP, ULPINV Finest relative precision and its inverse.
362 *> RTULP, RTULPI Square roots of the previous 4 values.
363 *>
364 *> The following four arrays decode JTYPE:
365 *> KTYPE(j) The general type (1-10) for type "j".
366 *> KMODE(j) The MODE value to be passed to the matrix
367 *> generator for type "j".
368 *> KMAGN(j) The order of magnitude ( O(1),
369 *> O(overflow^(1/2) ), O(underflow^(1/2) )
370 *> KCONDS(j) Selectw whether CONDS is to be 1 or
371 *> 1/sqrt(ulp). (0 means irrelevant.)
372 *> \endverbatim
373 *
374 * Authors:
375 * ========
376 *
377 *> \author Univ. of Tennessee
378 *> \author Univ. of California Berkeley
379 *> \author Univ. of Colorado Denver
380 *> \author NAG Ltd.
381 *
382 *> \date November 2011
383 *
384 *> \ingroup double_eig
385 *
386 * =====================================================================
387  SUBROUTINE ddrves( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
388  \$ nounit, a, lda, h, ht, wr, wi, wrt, wit, vs,
389  \$ ldvs, result, work, nwork, iwork, bwork, info )
390 *
391 * -- LAPACK test routine (version 3.4.0) --
392 * -- LAPACK is a software package provided by Univ. of Tennessee, --
393 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
394 * November 2011
395 *
396 * .. Scalar Arguments ..
397  INTEGER info, lda, ldvs, nounit, nsizes, ntypes, nwork
398  DOUBLE PRECISION thresh
399 * ..
400 * .. Array Arguments ..
401  LOGICAL bwork( * ), dotype( * )
402  INTEGER iseed( 4 ), iwork( * ), nn( * )
403  DOUBLE PRECISION a( lda, * ), h( lda, * ), ht( lda, * ),
404  \$ result( 13 ), vs( ldvs, * ), wi( * ), wit( * ),
405  \$ work( * ), wr( * ), wrt( * )
406 * ..
407 *
408 * =====================================================================
409 *
410 * .. Parameters ..
411  DOUBLE PRECISION zero, one
412  parameter( zero = 0.0d0, one = 1.0d0 )
413  INTEGER maxtyp
414  parameter( maxtyp = 21 )
415 * ..
416 * .. Local Scalars ..
418  CHARACTER sort
419  CHARACTER*3 path
420  INTEGER i, iinfo, imode, isort, itype, iwk, j, jcol,
421  \$ jsize, jtype, knteig, lwork, mtypes, n, nerrs,
422  \$ nfail, nmax, nnwork, ntest, ntestf, ntestt,
423  \$ rsub, sdim
424  DOUBLE PRECISION anorm, cond, conds, ovfl, rtulp, rtulpi, tmp,
425  \$ ulp, ulpinv, unfl
426 * ..
427 * .. Local Arrays ..
429  INTEGER idumma( 1 ), ioldsd( 4 ), kconds( maxtyp ),
430  \$ kmagn( maxtyp ), kmode( maxtyp ),
431  \$ ktype( maxtyp )
432  DOUBLE PRECISION res( 2 )
433 * ..
434 * .. Arrays in Common ..
435  LOGICAL selval( 20 )
436  DOUBLE PRECISION selwi( 20 ), selwr( 20 )
437 * ..
438 * .. Scalars in Common ..
439  INTEGER seldim, selopt
440 * ..
441 * .. Common blocks ..
442  common / sslct / selopt, seldim, selval, selwr, selwi
443 * ..
444 * .. External Functions ..
445  LOGICAL dslect
446  DOUBLE PRECISION dlamch
447  EXTERNAL dslect, dlamch
448 * ..
449 * .. External Subroutines ..
450  EXTERNAL dgees, dhst01, dlabad, dlacpy, dlaset, dlasum,
452 * ..
453 * .. Intrinsic Functions ..
454  INTRINSIC abs, max, sign, sqrt
455 * ..
456 * .. Data statements ..
457  DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
458  DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
459  \$ 3, 1, 2, 3 /
460  DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
461  \$ 1, 5, 5, 5, 4, 3, 1 /
462  DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
463 * ..
464 * .. Executable Statements ..
465 *
466  path( 1: 1 ) = 'Double precision'
467  path( 2: 3 ) = 'ES'
468 *
469 * Check for errors
470 *
471  ntestt = 0
472  ntestf = 0
473  info = 0
474  selopt = 0
475 *
476 * Important constants
477 *
479  nmax = 0
480  DO 10 j = 1, nsizes
481  nmax = max( nmax, nn( j ) )
482  IF( nn( j ).LT.0 )
484  10 continue
485 *
486 * Check for errors
487 *
488  IF( nsizes.LT.0 ) THEN
489  info = -1
490  ELSE IF( badnn ) THEN
491  info = -2
492  ELSE IF( ntypes.LT.0 ) THEN
493  info = -3
494  ELSE IF( thresh.LT.zero ) THEN
495  info = -6
496  ELSE IF( nounit.LE.0 ) THEN
497  info = -7
498  ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
499  info = -9
500  ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax ) THEN
501  info = -17
502  ELSE IF( 5*nmax+2*nmax**2.GT.nwork ) THEN
503  info = -20
504  END IF
505 *
506  IF( info.NE.0 ) THEN
507  CALL xerbla( 'DDRVES', -info )
508  return
509  END IF
510 *
511 * Quick return if nothing to do
512 *
513  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
514  \$ return
515 *
516 * More Important constants
517 *
518  unfl = dlamch( 'Safe minimum' )
519  ovfl = one / unfl
520  CALL dlabad( unfl, ovfl )
521  ulp = dlamch( 'Precision' )
522  ulpinv = one / ulp
523  rtulp = sqrt( ulp )
524  rtulpi = one / rtulp
525 *
526 * Loop over sizes, types
527 *
528  nerrs = 0
529 *
530  DO 270 jsize = 1, nsizes
531  n = nn( jsize )
532  mtypes = maxtyp
533  IF( nsizes.EQ.1 .AND. ntypes.EQ.maxtyp+1 )
534  \$ mtypes = mtypes + 1
535 *
536  DO 260 jtype = 1, mtypes
537  IF( .NOT.dotype( jtype ) )
538  \$ go to 260
539 *
540 * Save ISEED in case of an error.
541 *
542  DO 20 j = 1, 4
543  ioldsd( j ) = iseed( j )
544  20 continue
545 *
546 * Compute "A"
547 *
548 * Control parameters:
549 *
550 * KMAGN KCONDS KMODE KTYPE
551 * =1 O(1) 1 clustered 1 zero
552 * =2 large large clustered 2 identity
553 * =3 small exponential Jordan
554 * =4 arithmetic diagonal, (w/ eigenvalues)
555 * =5 random log symmetric, w/ eigenvalues
556 * =6 random general, w/ eigenvalues
557 * =7 random diagonal
558 * =8 random symmetric
559 * =9 random general
560 * =10 random triangular
561 *
562  IF( mtypes.GT.maxtyp )
563  \$ go to 90
564 *
565  itype = ktype( jtype )
566  imode = kmode( jtype )
567 *
568 * Compute norm
569 *
570  go to( 30, 40, 50 )kmagn( jtype )
571 *
572  30 continue
573  anorm = one
574  go to 60
575 *
576  40 continue
577  anorm = ovfl*ulp
578  go to 60
579 *
580  50 continue
581  anorm = unfl*ulpinv
582  go to 60
583 *
584  60 continue
585 *
586  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
587  iinfo = 0
588  cond = ulpinv
589 *
590 * Special Matrices -- Identity & Jordan block
591 *
592 * Zero
593 *
594  IF( itype.EQ.1 ) THEN
595  iinfo = 0
596 *
597  ELSE IF( itype.EQ.2 ) THEN
598 *
599 * Identity
600 *
601  DO 70 jcol = 1, n
602  a( jcol, jcol ) = anorm
603  70 continue
604 *
605  ELSE IF( itype.EQ.3 ) THEN
606 *
607 * Jordan Block
608 *
609  DO 80 jcol = 1, n
610  a( jcol, jcol ) = anorm
611  IF( jcol.GT.1 )
612  \$ a( jcol, jcol-1 ) = one
613  80 continue
614 *
615  ELSE IF( itype.EQ.4 ) THEN
616 *
617 * Diagonal Matrix, [Eigen]values Specified
618 *
619  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
620  \$ anorm, 0, 0, 'N', a, lda, work( n+1 ),
621  \$ iinfo )
622 *
623  ELSE IF( itype.EQ.5 ) THEN
624 *
625 * Symmetric, eigenvalues specified
626 *
627  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
628  \$ anorm, n, n, 'N', a, lda, work( n+1 ),
629  \$ iinfo )
630 *
631  ELSE IF( itype.EQ.6 ) THEN
632 *
633 * General, eigenvalues specified
634 *
635  IF( kconds( jtype ).EQ.1 ) THEN
636  conds = one
637  ELSE IF( kconds( jtype ).EQ.2 ) THEN
638  conds = rtulpi
639  ELSE
640  conds = zero
641  END IF
642 *
643  adumma( 1 ) = ' '
644  CALL dlatme( n, 'S', iseed, work, imode, cond, one,
645  \$ adumma, 'T', 'T', 'T', work( n+1 ), 4,
646  \$ conds, n, n, anorm, a, lda, work( 2*n+1 ),
647  \$ iinfo )
648 *
649  ELSE IF( itype.EQ.7 ) THEN
650 *
651 * Diagonal, random eigenvalues
652 *
653  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
654  \$ 'T', 'N', work( n+1 ), 1, one,
655  \$ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
656  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
657 *
658  ELSE IF( itype.EQ.8 ) THEN
659 *
660 * Symmetric, random eigenvalues
661 *
662  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
663  \$ 'T', 'N', work( n+1 ), 1, one,
664  \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
665  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
666 *
667  ELSE IF( itype.EQ.9 ) THEN
668 *
669 * General, random eigenvalues
670 *
671  CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
672  \$ 'T', 'N', work( n+1 ), 1, one,
673  \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
674  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
675  IF( n.GE.4 ) THEN
676  CALL dlaset( 'Full', 2, n, zero, zero, a, lda )
677  CALL dlaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
678  \$ lda )
679  CALL dlaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
680  \$ lda )
681  CALL dlaset( 'Full', 1, n, zero, zero, a( n, 1 ),
682  \$ lda )
683  END IF
684 *
685  ELSE IF( itype.EQ.10 ) THEN
686 *
687 * Triangular, random eigenvalues
688 *
689  CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
690  \$ 'T', 'N', work( n+1 ), 1, one,
691  \$ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
692  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
693 *
694  ELSE
695 *
696  iinfo = 1
697  END IF
698 *
699  IF( iinfo.NE.0 ) THEN
700  WRITE( nounit, fmt = 9992 )'Generator', iinfo, n, jtype,
701  \$ ioldsd
702  info = abs( iinfo )
703  return
704  END IF
705 *
706  90 continue
707 *
708 * Test for minimal and generous workspace
709 *
710  DO 250 iwk = 1, 2
711  IF( iwk.EQ.1 ) THEN
712  nnwork = 3*n
713  ELSE
714  nnwork = 5*n + 2*n**2
715  END IF
716  nnwork = max( nnwork, 1 )
717 *
718 * Initialize RESULT
719 *
720  DO 100 j = 1, 13
721  result( j ) = -one
722  100 continue
723 *
724 * Test with and without sorting of eigenvalues
725 *
726  DO 210 isort = 0, 1
727  IF( isort.EQ.0 ) THEN
728  sort = 'N'
729  rsub = 0
730  ELSE
731  sort = 'S'
732  rsub = 6
733  END IF
734 *
735 * Compute Schur form and Schur vectors, and test them
736 *
737  CALL dlacpy( 'F', n, n, a, lda, h, lda )
738  CALL dgees( 'V', sort, dslect, n, h, lda, sdim, wr,
739  \$ wi, vs, ldvs, work, nnwork, bwork, iinfo )
740  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
741  result( 1+rsub ) = ulpinv
742  WRITE( nounit, fmt = 9992 )'DGEES1', iinfo, n,
743  \$ jtype, ioldsd
744  info = abs( iinfo )
745  go to 220
746  END IF
747 *
748 * Do Test (1) or Test (7)
749 *
750  result( 1+rsub ) = zero
751  DO 120 j = 1, n - 2
752  DO 110 i = j + 2, n
753  IF( h( i, j ).NE.zero )
754  \$ result( 1+rsub ) = ulpinv
755  110 continue
756  120 continue
757  DO 130 i = 1, n - 2
758  IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.
759  \$ zero )result( 1+rsub ) = ulpinv
760  130 continue
761  DO 140 i = 1, n - 1
762  IF( h( i+1, i ).NE.zero ) THEN
763  IF( h( i, i ).NE.h( i+1, i+1 ) .OR.
764  \$ h( i, i+1 ).EQ.zero .OR.
765  \$ sign( one, h( i+1, i ) ).EQ.
766  \$ sign( one, h( i, i+1 ) ) )result( 1+rsub )
767  \$ = ulpinv
768  END IF
769  140 continue
770 *
771 * Do Tests (2) and (3) or Tests (8) and (9)
772 *
773  lwork = max( 1, 2*n*n )
774  CALL dhst01( n, 1, n, a, lda, h, lda, vs, ldvs, work,
775  \$ lwork, res )
776  result( 2+rsub ) = res( 1 )
777  result( 3+rsub ) = res( 2 )
778 *
779 * Do Test (4) or Test (10)
780 *
781  result( 4+rsub ) = zero
782  DO 150 i = 1, n
783  IF( h( i, i ).NE.wr( i ) )
784  \$ result( 4+rsub ) = ulpinv
785  150 continue
786  IF( n.GT.1 ) THEN
787  IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
788  \$ result( 4+rsub ) = ulpinv
789  IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
790  \$ result( 4+rsub ) = ulpinv
791  END IF
792  DO 160 i = 1, n - 1
793  IF( h( i+1, i ).NE.zero ) THEN
794  tmp = sqrt( abs( h( i+1, i ) ) )*
795  \$ sqrt( abs( h( i, i+1 ) ) )
796  result( 4+rsub ) = max( result( 4+rsub ),
797  \$ abs( wi( i )-tmp ) /
798  \$ max( ulp*tmp, unfl ) )
799  result( 4+rsub ) = max( result( 4+rsub ),
800  \$ abs( wi( i+1 )+tmp ) /
801  \$ max( ulp*tmp, unfl ) )
802  ELSE IF( i.GT.1 ) THEN
803  IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.
804  \$ zero .AND. wi( i ).NE.zero )result( 4+rsub )
805  \$ = ulpinv
806  END IF
807  160 continue
808 *
809 * Do Test (5) or Test (11)
810 *
811  CALL dlacpy( 'F', n, n, a, lda, ht, lda )
812  CALL dgees( 'N', sort, dslect, n, ht, lda, sdim, wrt,
813  \$ wit, vs, ldvs, work, nnwork, bwork,
814  \$ iinfo )
815  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
816  result( 5+rsub ) = ulpinv
817  WRITE( nounit, fmt = 9992 )'DGEES2', iinfo, n,
818  \$ jtype, ioldsd
819  info = abs( iinfo )
820  go to 220
821  END IF
822 *
823  result( 5+rsub ) = zero
824  DO 180 j = 1, n
825  DO 170 i = 1, n
826  IF( h( i, j ).NE.ht( i, j ) )
827  \$ result( 5+rsub ) = ulpinv
828  170 continue
829  180 continue
830 *
831 * Do Test (6) or Test (12)
832 *
833  result( 6+rsub ) = zero
834  DO 190 i = 1, n
835  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
836  \$ result( 6+rsub ) = ulpinv
837  190 continue
838 *
839 * Do Test (13)
840 *
841  IF( isort.EQ.1 ) THEN
842  result( 13 ) = zero
843  knteig = 0
844  DO 200 i = 1, n
845  IF( dslect( wr( i ), wi( i ) ) .OR.
846  \$ dslect( wr( i ), -wi( i ) ) )
847  \$ knteig = knteig + 1
848  IF( i.LT.n ) THEN
849  IF( ( dslect( wr( i+1 ),
850  \$ wi( i+1 ) ) .OR. dslect( wr( i+1 ),
851  \$ -wi( i+1 ) ) ) .AND.
852  \$ ( .NOT.( dslect( wr( i ),
853  \$ wi( i ) ) .OR. dslect( wr( i ),
854  \$ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )
855  \$ result( 13 ) = ulpinv
856  END IF
857  200 continue
858  IF( sdim.NE.knteig ) THEN
859  result( 13 ) = ulpinv
860  END IF
861  END IF
862 *
863  210 continue
864 *
865 * End of Loop -- Check for RESULT(j) > THRESH
866 *
867  220 continue
868 *
869  ntest = 0
870  nfail = 0
871  DO 230 j = 1, 13
872  IF( result( j ).GE.zero )
873  \$ ntest = ntest + 1
874  IF( result( j ).GE.thresh )
875  \$ nfail = nfail + 1
876  230 continue
877 *
878  IF( nfail.GT.0 )
879  \$ ntestf = ntestf + 1
880  IF( ntestf.EQ.1 ) THEN
881  WRITE( nounit, fmt = 9999 )path
882  WRITE( nounit, fmt = 9998 )
883  WRITE( nounit, fmt = 9997 )
884  WRITE( nounit, fmt = 9996 )
885  WRITE( nounit, fmt = 9995 )thresh
886  WRITE( nounit, fmt = 9994 )
887  ntestf = 2
888  END IF
889 *
890  DO 240 j = 1, 13
891  IF( result( j ).GE.thresh ) THEN
892  WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
893  \$ j, result( j )
894  END IF
895  240 continue
896 *
897  nerrs = nerrs + nfail
898  ntestt = ntestt + ntest
899 *
900  250 continue
901  260 continue
902  270 continue
903 *
904 * Summary
905 *
906  CALL dlasum( path, nounit, nerrs, ntestt )
907 *
908  9999 format( / 1x, a3, ' -- Real Schur Form Decomposition Driver',
909  \$ / ' Matrix types (see DDRVES for details): ' )
910 *
911  9998 format( / ' Special Matrices:', / ' 1=Zero matrix. ',
912  \$ ' ', ' 5=Diagonal: geometr. spaced entries.',
913  \$ / ' 2=Identity matrix. ', ' 6=Diagona',
914  \$ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
915  \$ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
916  \$ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
917  \$ 'mall, evenly spaced.' )
918  9997 format( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
919  \$ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
920  \$ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
921  \$ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
922  \$ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
923  \$ 'lex ', / ' 12=Well-cond., random complex ', 6x, ' ',
924  \$ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
925  \$ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
926  \$ ' complx ' )
927  9996 format( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
928  \$ 'with small random entries.', / ' 20=Matrix with large ran',
929  \$ 'dom entries. ', / )
930  9995 format( ' Tests performed with test threshold =', f8.2,
931  \$ / ' ( A denotes A on input and T denotes A on output)',
932  \$ / / ' 1 = 0 if T in Schur form (no sort), ',
933  \$ ' 1/ulp otherwise', /
934  \$ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
935  \$ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
936  \$ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
937  \$ ' 1/ulp otherwise', /
938  \$ ' 5 = 0 if T same no matter if VS computed (no sort),',
939  \$ ' 1/ulp otherwise', /
940  \$ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
941  \$ ', 1/ulp otherwise' )
942  9994 format( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
943  \$ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
944  \$ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
945  \$ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
946  \$ ' 1/ulp otherwise', /
947  \$ ' 11 = 0 if T same no matter if VS computed (sort),',
948  \$ ' 1/ulp otherwise', /
949  \$ ' 12 = 0 if WR, WI same no matter if VS computed (sort),',
950  \$ ' 1/ulp otherwise', /
951  \$ ' 13 = 0 if sorting succesful, 1/ulp otherwise', / )
952  9993 format( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
953  \$ ' type ', i2, ', test(', i2, ')=', g10.3 )
954  9992 format( ' DDRVES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
955  \$ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
956 *
957  return
958 *
959 * End of DDRVES
960 *
961  END