LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cdrvsx.f
Go to the documentation of this file.
1 *> \brief \b CDRVSX
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 CDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS,
13 * LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK,
14 * INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
18 * $ NTYPES
19 * REAL THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL BWORK( * ), DOTYPE( * )
23 * INTEGER ISEED( 4 ), NN( * )
24 * REAL RESULT( 17 ), RWORK( * )
25 * COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
26 * $ VS( LDVS, * ), VS1( LDVS, * ), W( * ),
27 * $ WORK( * ), WT( * ), WTMP( * )
28 * ..
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> CDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
37 *> expert driver CGEESX.
38 *>
39 *> CDRVSX 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 CDRVSX 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 W 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 W are eigenvalues of T
83 *> 1/ulp otherwise
84 *> If workspace sufficient, also compare W 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 complex angles.
121 *> (ULP = (first number larger than 1) - 1 )
122 *> (5) A diagonal matrix with geometrically spaced entries
123 *> 1, ..., ULP and random complex angles.
124 *> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
125 *> and random complex angles.
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 unitary and
133 *> T has evenly spaced entries 1, ..., ULP with random
134 *> complex angles on the diagonal and random O(1) entries in
135 *> the upper triangle.
136 *>
137 *> (10) A matrix of the form U' T U, where U is unitary and
138 *> T has geometrically spaced entries 1, ..., ULP with random
139 *> complex angles on the diagonal and random O(1) entries in
140 *> the upper 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 *> complex angles on the diagonal and random O(1) entries in
145 *> the upper triangle.
146 *>
147 *> (12) A matrix of the form U' T U, where U is unitary and
148 *> T has complex eigenvalues randomly chosen from
149 *> ULP < |z| < 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 complex angles on the diagonal and random O(1)
155 *> entries 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 complex angles on the diagonal
160 *> and random 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 complex angles on the diagonal and random O(1)
165 *> entries 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 complex eigenvalues randomly chosen
169 *> from ULP < |z| < 1 and random O(1) entries in the upper
170 *> 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 CGEESX 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 CGEESX 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 CDRVSX to continue the same random number
266 *> sequence.
267 *> \endverbatim
268 *>
269 *> \param[in] THRESH
270 *> \verbatim
271 *> THRESH is REAL
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 COMPLEX 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 COMPLEX array, dimension (LDA, max(NN))
311 *> Another copy of the test matrix A, modified by CGEESX.
312 *> \endverbatim
313 *>
314 *> \param[out] HT
315 *> \verbatim
316 *> HT is COMPLEX array, dimension (LDA, max(NN))
317 *> Yet another copy of the test matrix A, modified by CGEESX.
318 *> \endverbatim
319 *>
320 *> \param[out] W
321 *> \verbatim
322 *> W is COMPLEX array, dimension (max(NN))
323 *> The computed eigenvalues of A.
324 *> \endverbatim
325 *>
326 *> \param[out] WT
327 *> \verbatim
328 *> WT is COMPLEX array, dimension (max(NN))
329 *> Like W, this array contains the eigenvalues of A,
330 *> but those computed when CGEESX only computes a partial
331 *> eigendecomposition, i.e. not Schur vectors
332 *> \endverbatim
333 *>
334 *> \param[out] WTMP
335 *> \verbatim
336 *> WTMP is COMPLEX array, dimension (max(NN))
337 *> More temporary storage for eigenvalues.
338 *> \endverbatim
339 *>
340 *> \param[out] VS
341 *> \verbatim
342 *> VS is COMPLEX array, dimension (LDVS, max(NN))
343 *> VS holds the computed Schur vectors.
344 *> \endverbatim
345 *>
346 *> \param[in] LDVS
347 *> \verbatim
348 *> LDVS is INTEGER
349 *> Leading dimension of VS. Must be at least max(1,max(NN)).
350 *> \endverbatim
351 *>
352 *> \param[out] VS1
353 *> \verbatim
354 *> VS1 is COMPLEX array, dimension (LDVS, max(NN))
355 *> VS1 holds another copy of the computed Schur vectors.
356 *> \endverbatim
357 *>
358 *> \param[out] RESULT
359 *> \verbatim
360 *> RESULT is REAL array, dimension (17)
361 *> The values computed by the 17 tests described above.
362 *> The values are currently limited to 1/ulp, to avoid overflow.
363 *> \endverbatim
364 *>
365 *> \param[out] WORK
366 *> \verbatim
367 *> WORK is COMPLEX array, dimension (LWORK)
368 *> \endverbatim
369 *>
370 *> \param[in] LWORK
371 *> \verbatim
372 *> LWORK is INTEGER
373 *> The number of entries in WORK. This must be at least
374 *> max(1,2*NN(j)**2) for all j.
375 *> \endverbatim
376 *>
377 *> \param[out] RWORK
378 *> \verbatim
379 *> RWORK is REAL array, dimension (max(NN))
380 *> \endverbatim
381 *>
382 *> \param[out] BWORK
383 *> \verbatim
384 *> BWORK is LOGICAL array, dimension (max(NN))
385 *> \endverbatim
386 *>
387 *> \param[out] INFO
388 *> \verbatim
389 *> INFO is INTEGER
390 *> If 0, successful exit.
391 *> <0, input parameter -INFO is incorrect
392 *> >0, CLATMR, CLATMS, CLATME or CGET24 returned an error
393 *> code and INFO is its absolute value
394 *>
395 *>-----------------------------------------------------------------------
396 *>
397 *> Some Local Variables and Parameters:
398 *> ---- ----- --------- --- ----------
399 *> ZERO, ONE Real 0 and 1.
400 *> MAXTYP The number of types defined.
401 *> NMAX Largest value in NN.
402 *> NERRS The number of tests which have exceeded THRESH
403 *> COND, CONDS,
404 *> IMODE Values to be passed to the matrix generators.
405 *> ANORM Norm of A; passed to matrix generators.
406 *>
407 *> OVFL, UNFL Overflow and underflow thresholds.
408 *> ULP, ULPINV Finest relative precision and its inverse.
409 *> RTULP, RTULPI Square roots of the previous 4 values.
410 *> The following four arrays decode JTYPE:
411 *> KTYPE(j) The general type (1-10) for type "j".
412 *> KMODE(j) The MODE value to be passed to the matrix
413 *> generator for type "j".
414 *> KMAGN(j) The order of magnitude ( O(1),
415 *> O(overflow^(1/2) ), O(underflow^(1/2) )
416 *> KCONDS(j) Selectw whether CONDS is to be 1 or
417 *> 1/sqrt(ulp). (0 means irrelevant.)
418 *> \endverbatim
419 *
420 * Authors:
421 * ========
422 *
423 *> \author Univ. of Tennessee
424 *> \author Univ. of California Berkeley
425 *> \author Univ. of Colorado Denver
426 *> \author NAG Ltd.
427 *
428 *> \date November 2011
429 *
430 *> \ingroup complex_eig
431 *
432 * =====================================================================
433  SUBROUTINE cdrvsx( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
434  $ niunit, nounit, a, lda, h, ht, w, wt, wtmp, vs,
435  $ ldvs, vs1, result, work, lwork, rwork, bwork,
436  $ info )
437 *
438 * -- LAPACK test routine (version 3.4.0) --
439 * -- LAPACK is a software package provided by Univ. of Tennessee, --
440 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
441 * November 2011
442 *
443 * .. Scalar Arguments ..
444  INTEGER info, lda, ldvs, lwork, niunit, nounit, nsizes,
445  $ ntypes
446  REAL thresh
447 * ..
448 * .. Array Arguments ..
449  LOGICAL bwork( * ), dotype( * )
450  INTEGER iseed( 4 ), nn( * )
451  REAL result( 17 ), rwork( * )
452  COMPLEX a( lda, * ), h( lda, * ), ht( lda, * ),
453  $ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
454  $ work( * ), wt( * ), wtmp( * )
455 * ..
456 *
457 * =====================================================================
458 *
459 * .. Parameters ..
460  COMPLEX czero
461  parameter( czero = ( 0.0e+0, 0.0e+0 ) )
462  COMPLEX cone
463  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
464  REAL zero, one
465  parameter( zero = 0.0e+0, one = 1.0e+0 )
466  INTEGER maxtyp
467  parameter( maxtyp = 21 )
468 * ..
469 * .. Local Scalars ..
470  LOGICAL badnn
471  CHARACTER*3 path
472  INTEGER i, iinfo, imode, isrt, itype, iwk, j, jcol,
473  $ jsize, jtype, mtypes, n, nerrs, nfail,
474  $ nmax, nnwork, nslct, ntest, ntestf, ntestt
475  REAL anorm, cond, conds, ovfl, rcdein, rcdvin,
476  $ rtulp, rtulpi, ulp, ulpinv, unfl
477 * ..
478 * .. Local Arrays ..
479  INTEGER idumma( 1 ), ioldsd( 4 ), islct( 20 ),
480  $ kconds( maxtyp ), kmagn( maxtyp ),
481  $ kmode( maxtyp ), ktype( maxtyp )
482 * ..
483 * .. Arrays in Common ..
484  LOGICAL selval( 20 )
485  REAL selwi( 20 ), selwr( 20 )
486 * ..
487 * .. Scalars in Common ..
488  INTEGER seldim, selopt
489 * ..
490 * .. Common blocks ..
491  common / sslct / selopt, seldim, selval, selwr, selwi
492 * ..
493 * .. External Functions ..
494  REAL slamch
495  EXTERNAL slamch
496 * ..
497 * .. External Subroutines ..
498  EXTERNAL cget24, clatme, clatmr, clatms, claset, slabad,
499  $ slasum, xerbla
500 * ..
501 * .. Intrinsic Functions ..
502  INTRINSIC abs, max, min, sqrt
503 * ..
504 * .. Data statements ..
505  DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
506  DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
507  $ 3, 1, 2, 3 /
508  DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
509  $ 1, 5, 5, 5, 4, 3, 1 /
510  DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
511 * ..
512 * .. Executable Statements ..
513 *
514  path( 1: 1 ) = 'Complex precision'
515  path( 2: 3 ) = 'SX'
516 *
517 * Check for errors
518 *
519  ntestt = 0
520  ntestf = 0
521  info = 0
522 *
523 * Important constants
524 *
525  badnn = .false.
526 *
527 * 8 is the largest dimension in the input file of precomputed
528 * problems
529 *
530  nmax = 8
531  DO 10 j = 1, nsizes
532  nmax = max( nmax, nn( j ) )
533  IF( nn( j ).LT.0 )
534  $ badnn = .true.
535  10 continue
536 *
537 * Check for errors
538 *
539  IF( nsizes.LT.0 ) THEN
540  info = -1
541  ELSE IF( badnn ) THEN
542  info = -2
543  ELSE IF( ntypes.LT.0 ) THEN
544  info = -3
545  ELSE IF( thresh.LT.zero ) THEN
546  info = -6
547  ELSE IF( niunit.LE.0 ) THEN
548  info = -7
549  ELSE IF( nounit.LE.0 ) THEN
550  info = -8
551  ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
552  info = -10
553  ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax ) THEN
554  info = -20
555  ELSE IF( max( 3*nmax, 2*nmax**2 ).GT.lwork ) THEN
556  info = -24
557  END IF
558 *
559  IF( info.NE.0 ) THEN
560  CALL xerbla( 'CDRVSX', -info )
561  return
562  END IF
563 *
564 * If nothing to do check on NIUNIT
565 *
566  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
567  $ go to 150
568 *
569 * More Important constants
570 *
571  unfl = slamch( 'Safe minimum' )
572  ovfl = one / unfl
573  CALL slabad( unfl, ovfl )
574  ulp = slamch( 'Precision' )
575  ulpinv = one / ulp
576  rtulp = sqrt( ulp )
577  rtulpi = one / rtulp
578 *
579 * Loop over sizes, types
580 *
581  nerrs = 0
582 *
583  DO 140 jsize = 1, nsizes
584  n = nn( jsize )
585  IF( nsizes.NE.1 ) THEN
586  mtypes = min( maxtyp, ntypes )
587  ELSE
588  mtypes = min( maxtyp+1, ntypes )
589  END IF
590 *
591  DO 130 jtype = 1, mtypes
592  IF( .NOT.dotype( jtype ) )
593  $ go to 130
594 *
595 * Save ISEED in case of an error.
596 *
597  DO 20 j = 1, 4
598  ioldsd( j ) = iseed( j )
599  20 continue
600 *
601 * Compute "A"
602 *
603 * Control parameters:
604 *
605 * KMAGN KCONDS KMODE KTYPE
606 * =1 O(1) 1 clustered 1 zero
607 * =2 large large clustered 2 identity
608 * =3 small exponential Jordan
609 * =4 arithmetic diagonal, (w/ eigenvalues)
610 * =5 random log symmetric, w/ eigenvalues
611 * =6 random general, w/ eigenvalues
612 * =7 random diagonal
613 * =8 random symmetric
614 * =9 random general
615 * =10 random triangular
616 *
617  IF( mtypes.GT.maxtyp )
618  $ go to 90
619 *
620  itype = ktype( jtype )
621  imode = kmode( jtype )
622 *
623 * Compute norm
624 *
625  go to( 30, 40, 50 )kmagn( jtype )
626 *
627  30 continue
628  anorm = one
629  go to 60
630 *
631  40 continue
632  anorm = ovfl*ulp
633  go to 60
634 *
635  50 continue
636  anorm = unfl*ulpinv
637  go to 60
638 *
639  60 continue
640 *
641  CALL claset( 'Full', lda, n, czero, czero, a, lda )
642  iinfo = 0
643  cond = ulpinv
644 *
645 * Special Matrices -- Identity & Jordan block
646 *
647  IF( itype.EQ.1 ) THEN
648 *
649 * Zero
650 *
651  iinfo = 0
652 *
653  ELSE IF( itype.EQ.2 ) THEN
654 *
655 * Identity
656 *
657  DO 70 jcol = 1, n
658  a( jcol, jcol ) = anorm
659  70 continue
660 *
661  ELSE IF( itype.EQ.3 ) THEN
662 *
663 * Jordan Block
664 *
665  DO 80 jcol = 1, n
666  a( jcol, jcol ) = anorm
667  IF( jcol.GT.1 )
668  $ a( jcol, jcol-1 ) = cone
669  80 continue
670 *
671  ELSE IF( itype.EQ.4 ) THEN
672 *
673 * Diagonal Matrix, [Eigen]values Specified
674 *
675  CALL clatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
676  $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
677  $ iinfo )
678 *
679  ELSE IF( itype.EQ.5 ) THEN
680 *
681 * Symmetric, eigenvalues specified
682 *
683  CALL clatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
684  $ anorm, n, n, 'N', a, lda, work( n+1 ),
685  $ iinfo )
686 *
687  ELSE IF( itype.EQ.6 ) THEN
688 *
689 * General, eigenvalues specified
690 *
691  IF( kconds( jtype ).EQ.1 ) THEN
692  conds = one
693  ELSE IF( kconds( jtype ).EQ.2 ) THEN
694  conds = rtulpi
695  ELSE
696  conds = zero
697  END IF
698 *
699  CALL clatme( n, 'D', iseed, work, imode, cond, cone,
700  $ 'T', 'T', 'T', rwork, 4, conds, n, n, anorm,
701  $ a, lda, work( 2*n+1 ), iinfo )
702 *
703  ELSE IF( itype.EQ.7 ) THEN
704 *
705 * Diagonal, random eigenvalues
706 *
707  CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
708  $ 'T', 'N', work( n+1 ), 1, one,
709  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
710  $ zero, anorm, 'NO', a, lda, idumma, iinfo )
711 *
712  ELSE IF( itype.EQ.8 ) THEN
713 *
714 * Symmetric, random eigenvalues
715 *
716  CALL clatmr( n, n, 'D', iseed, 'H', work, 6, one, cone,
717  $ 'T', 'N', work( n+1 ), 1, one,
718  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
719  $ zero, anorm, 'NO', a, lda, idumma, iinfo )
720 *
721  ELSE IF( itype.EQ.9 ) THEN
722 *
723 * General, random eigenvalues
724 *
725  CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
726  $ 'T', 'N', work( n+1 ), 1, one,
727  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
728  $ zero, anorm, 'NO', a, lda, idumma, iinfo )
729  IF( n.GE.4 ) THEN
730  CALL claset( 'Full', 2, n, czero, czero, a, lda )
731  CALL claset( 'Full', n-3, 1, czero, czero, a( 3, 1 ),
732  $ lda )
733  CALL claset( 'Full', n-3, 2, czero, czero,
734  $ a( 3, n-1 ), lda )
735  CALL claset( 'Full', 1, n, czero, czero, a( n, 1 ),
736  $ lda )
737  END IF
738 *
739  ELSE IF( itype.EQ.10 ) THEN
740 *
741 * Triangular, random eigenvalues
742 *
743  CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
744  $ 'T', 'N', work( n+1 ), 1, one,
745  $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
746  $ zero, anorm, 'NO', a, lda, idumma, iinfo )
747 *
748  ELSE
749 *
750  iinfo = 1
751  END IF
752 *
753  IF( iinfo.NE.0 ) THEN
754  WRITE( nounit, fmt = 9991 )'Generator', iinfo, n, jtype,
755  $ ioldsd
756  info = abs( iinfo )
757  return
758  END IF
759 *
760  90 continue
761 *
762 * Test for minimal and generous workspace
763 *
764  DO 120 iwk = 1, 2
765  IF( iwk.EQ.1 ) THEN
766  nnwork = 2*n
767  ELSE
768  nnwork = max( 2*n, n*( n+1 ) / 2 )
769  END IF
770  nnwork = max( nnwork, 1 )
771 *
772  CALL cget24( .false., jtype, thresh, ioldsd, nounit, n,
773  $ a, lda, h, ht, w, wt, wtmp, vs, ldvs, vs1,
774  $ rcdein, rcdvin, nslct, islct, 0, result,
775  $ work, nnwork, rwork, bwork, info )
776 *
777 * Check for RESULT(j) > THRESH
778 *
779  ntest = 0
780  nfail = 0
781  DO 100 j = 1, 15
782  IF( result( j ).GE.zero )
783  $ ntest = ntest + 1
784  IF( result( j ).GE.thresh )
785  $ nfail = nfail + 1
786  100 continue
787 *
788  IF( nfail.GT.0 )
789  $ ntestf = ntestf + 1
790  IF( ntestf.EQ.1 ) THEN
791  WRITE( nounit, fmt = 9999 )path
792  WRITE( nounit, fmt = 9998 )
793  WRITE( nounit, fmt = 9997 )
794  WRITE( nounit, fmt = 9996 )
795  WRITE( nounit, fmt = 9995 )thresh
796  WRITE( nounit, fmt = 9994 )
797  ntestf = 2
798  END IF
799 *
800  DO 110 j = 1, 15
801  IF( result( j ).GE.thresh ) THEN
802  WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
803  $ j, result( j )
804  END IF
805  110 continue
806 *
807  nerrs = nerrs + nfail
808  ntestt = ntestt + ntest
809 *
810  120 continue
811  130 continue
812  140 continue
813 *
814  150 continue
815 *
816 * Read in data from file to check accuracy of condition estimation
817 * Read input data until N=0
818 *
819  jtype = 0
820  160 continue
821  READ( niunit, fmt = *, END = 200 )n, nslct, isrt
822  IF( n.EQ.0 )
823  $ go to 200
824  jtype = jtype + 1
825  iseed( 1 ) = jtype
826  READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
827  DO 170 i = 1, n
828  READ( niunit, fmt = * )( a( i, j ), j = 1, n )
829  170 continue
830  READ( niunit, fmt = * )rcdein, rcdvin
831 *
832  CALL cget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
833  $ w, wt, wtmp, vs, ldvs, vs1, rcdein, rcdvin, nslct,
834  $ islct, isrt, result, work, lwork, rwork, bwork,
835  $ info )
836 *
837 * Check for RESULT(j) > THRESH
838 *
839  ntest = 0
840  nfail = 0
841  DO 180 j = 1, 17
842  IF( result( j ).GE.zero )
843  $ ntest = ntest + 1
844  IF( result( j ).GE.thresh )
845  $ nfail = nfail + 1
846  180 continue
847 *
848  IF( nfail.GT.0 )
849  $ ntestf = ntestf + 1
850  IF( ntestf.EQ.1 ) THEN
851  WRITE( nounit, fmt = 9999 )path
852  WRITE( nounit, fmt = 9998 )
853  WRITE( nounit, fmt = 9997 )
854  WRITE( nounit, fmt = 9996 )
855  WRITE( nounit, fmt = 9995 )thresh
856  WRITE( nounit, fmt = 9994 )
857  ntestf = 2
858  END IF
859  DO 190 j = 1, 17
860  IF( result( j ).GE.thresh ) THEN
861  WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
862  END IF
863  190 continue
864 *
865  nerrs = nerrs + nfail
866  ntestt = ntestt + ntest
867  go to 160
868  200 continue
869 *
870 * Summary
871 *
872  CALL slasum( path, nounit, nerrs, ntestt )
873 *
874  9999 format( / 1x, a3, ' -- Complex Schur Form Decomposition Expert ',
875  $ 'Driver', / ' Matrix types (see CDRVSX for details): ' )
876 *
877  9998 format( / ' Special Matrices:', / ' 1=Zero matrix. ',
878  $ ' ', ' 5=Diagonal: geometr. spaced entries.',
879  $ / ' 2=Identity matrix. ', ' 6=Diagona',
880  $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
881  $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
882  $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
883  $ 'mall, evenly spaced.' )
884  9997 format( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
885  $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
886  $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
887  $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
888  $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
889  $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
890  $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
891  $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
892  $ ' complx ' )
893  9996 format( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
894  $ 'with small random entries.', / ' 20=Matrix with large ran',
895  $ 'dom entries. ', / )
896  9995 format( ' Tests performed with test threshold =', f8.2,
897  $ / ' ( A denotes A on input and T denotes A on output)',
898  $ / / ' 1 = 0 if T in Schur form (no sort), ',
899  $ ' 1/ulp otherwise', /
900  $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
901  $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
902  $ / ' 4 = 0 if W are eigenvalues of T (no sort),',
903  $ ' 1/ulp otherwise', /
904  $ ' 5 = 0 if T same no matter if VS computed (no sort),',
905  $ ' 1/ulp otherwise', /
906  $ ' 6 = 0 if W same no matter if VS computed (no sort)',
907  $ ', 1/ulp otherwise' )
908  9994 format( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
909  $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
910  $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
911  $ / ' 10 = 0 if W are eigenvalues of T (sort),',
912  $ ' 1/ulp otherwise', /
913  $ ' 11 = 0 if T same no matter what else computed (sort),',
914  $ ' 1/ulp otherwise', /
915  $ ' 12 = 0 if W same no matter what else computed ',
916  $ '(sort), 1/ulp otherwise', /
917  $ ' 13 = 0 if sorting succesful, 1/ulp otherwise',
918  $ / ' 14 = 0 if RCONDE same no matter what else computed,',
919  $ ' 1/ulp otherwise', /
920  $ ' 15 = 0 if RCONDv same no matter what else computed,',
921  $ ' 1/ulp otherwise', /
922  $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
923  $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
924  9993 format( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
925  $ ' type ', i2, ', test(', i2, ')=', g10.3 )
926  9992 format( ' N=', i5, ', input example =', i3, ', test(', i2, ')=',
927  $ g10.3 )
928  9991 format( ' CDRVSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
929  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
930 *
931  return
932 *
933 * End of CDRVSX
934 *
935  END