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