LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \ingroup complex_eig
429*
430* =====================================================================
431 SUBROUTINE cdrvsx( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
432 $ NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS,
433 $ LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK,
434 $ INFO )
435*
436* -- LAPACK test routine --
437* -- LAPACK is a software package provided by Univ. of Tennessee, --
438* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
439*
440* .. Scalar Arguments ..
441 INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
442 $ NTYPES
443 REAL THRESH
444* ..
445* .. Array Arguments ..
446 LOGICAL BWORK( * ), DOTYPE( * )
447 INTEGER ISEED( 4 ), NN( * )
448 REAL RESULT( 17 ), RWORK( * )
449 COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
450 $ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
451 $ work( * ), wt( * ), wtmp( * )
452* ..
453*
454* =====================================================================
455*
456* .. Parameters ..
457 COMPLEX CZERO
458 PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ) )
459 COMPLEX CONE
460 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
461 REAL ZERO, ONE
462 parameter( zero = 0.0e+0, one = 1.0e+0 )
463 INTEGER MAXTYP
464 parameter( maxtyp = 21 )
465* ..
466* .. Local Scalars ..
467 LOGICAL BADNN
468 CHARACTER*3 PATH
469 INTEGER I, IINFO, IMODE, ISRT, ITYPE, IWK, J, JCOL,
470 $ jsize, jtype, mtypes, n, nerrs, nfail,
471 $ nmax, nnwork, nslct, ntest, ntestf, ntestt
472 REAL ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
473 $ RTULP, RTULPI, ULP, ULPINV, UNFL
474* ..
475* .. Local Arrays ..
476 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
477 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ),
478 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
479* ..
480* .. Arrays in Common ..
481 LOGICAL SELVAL( 20 )
482 REAL SELWI( 20 ), SELWR( 20 )
483* ..
484* .. Scalars in Common ..
485 INTEGER SELDIM, SELOPT
486* ..
487* .. Common blocks ..
488 COMMON / sslct / selopt, seldim, selval, selwr, selwi
489* ..
490* .. External Functions ..
491 REAL SLAMCH
492 EXTERNAL SLAMCH
493* ..
494* .. External Subroutines ..
495 EXTERNAL cget24, clatme, clatmr, clatms, claset,
496 $ slasum, xerbla
497* ..
498* .. Intrinsic Functions ..
499 INTRINSIC abs, max, min, sqrt
500* ..
501* .. Data statements ..
502 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
503 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
504 $ 3, 1, 2, 3 /
505 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
506 $ 1, 5, 5, 5, 4, 3, 1 /
507 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
508* ..
509* .. Executable Statements ..
510*
511 path( 1: 1 ) = 'Complex precision'
512 path( 2: 3 ) = 'SX'
513*
514* Check for errors
515*
516 ntestt = 0
517 ntestf = 0
518 info = 0
519*
520* Important constants
521*
522 badnn = .false.
523*
524* 8 is the largest dimension in the input file of precomputed
525* problems
526*
527 nmax = 8
528 DO 10 j = 1, nsizes
529 nmax = max( nmax, nn( j ) )
530 IF( nn( j ).LT.0 )
531 $ badnn = .true.
532 10 CONTINUE
533*
534* Check for errors
535*
536 IF( nsizes.LT.0 ) THEN
537 info = -1
538 ELSE IF( badnn ) THEN
539 info = -2
540 ELSE IF( ntypes.LT.0 ) THEN
541 info = -3
542 ELSE IF( thresh.LT.zero ) THEN
543 info = -6
544 ELSE IF( niunit.LE.0 ) THEN
545 info = -7
546 ELSE IF( nounit.LE.0 ) THEN
547 info = -8
548 ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
549 info = -10
550 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax ) THEN
551 info = -20
552 ELSE IF( max( 3*nmax, 2*nmax**2 ).GT.lwork ) THEN
553 info = -24
554 END IF
555*
556 IF( info.NE.0 ) THEN
557 CALL xerbla( 'CDRVSX', -info )
558 RETURN
559 END IF
560*
561* If nothing to do check on NIUNIT
562*
563 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
564 $ GO TO 150
565*
566* More Important constants
567*
568 unfl = slamch( 'Safe minimum' )
569 ovfl = one / unfl
570 ulp = slamch( 'Precision' )
571 ulpinv = one / ulp
572 rtulp = sqrt( ulp )
573 rtulpi = one / rtulp
574*
575* Loop over sizes, types
576*
577 nerrs = 0
578*
579 DO 140 jsize = 1, nsizes
580 n = nn( jsize )
581 IF( nsizes.NE.1 ) THEN
582 mtypes = min( maxtyp, ntypes )
583 ELSE
584 mtypes = min( maxtyp+1, ntypes )
585 END IF
586*
587 DO 130 jtype = 1, mtypes
588 IF( .NOT.dotype( jtype ) )
589 $ GO TO 130
590*
591* Save ISEED in case of an error.
592*
593 DO 20 j = 1, 4
594 ioldsd( j ) = iseed( j )
595 20 CONTINUE
596*
597* Compute "A"
598*
599* Control parameters:
600*
601* KMAGN KCONDS KMODE KTYPE
602* =1 O(1) 1 clustered 1 zero
603* =2 large large clustered 2 identity
604* =3 small exponential Jordan
605* =4 arithmetic diagonal, (w/ eigenvalues)
606* =5 random log symmetric, w/ eigenvalues
607* =6 random general, w/ eigenvalues
608* =7 random diagonal
609* =8 random symmetric
610* =9 random general
611* =10 random triangular
612*
613 IF( mtypes.GT.maxtyp )
614 $ GO TO 90
615*
616 itype = ktype( jtype )
617 imode = kmode( jtype )
618*
619* Compute norm
620*
621 GO TO ( 30, 40, 50 )kmagn( jtype )
622*
623 30 CONTINUE
624 anorm = one
625 GO TO 60
626*
627 40 CONTINUE
628 anorm = ovfl*ulp
629 GO TO 60
630*
631 50 CONTINUE
632 anorm = unfl*ulpinv
633 GO TO 60
634*
635 60 CONTINUE
636*
637 CALL claset( 'Full', lda, n, czero, czero, a, lda )
638 iinfo = 0
639 cond = ulpinv
640*
641* Special Matrices -- Identity & Jordan block
642*
643 IF( itype.EQ.1 ) THEN
644*
645* Zero
646*
647 iinfo = 0
648*
649 ELSE IF( itype.EQ.2 ) THEN
650*
651* Identity
652*
653 DO 70 jcol = 1, n
654 a( jcol, jcol ) = anorm
655 70 CONTINUE
656*
657 ELSE IF( itype.EQ.3 ) THEN
658*
659* Jordan Block
660*
661 DO 80 jcol = 1, n
662 a( jcol, jcol ) = anorm
663 IF( jcol.GT.1 )
664 $ a( jcol, jcol-1 ) = cone
665 80 CONTINUE
666*
667 ELSE IF( itype.EQ.4 ) THEN
668*
669* Diagonal Matrix, [Eigen]values Specified
670*
671 CALL clatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
672 $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
673 $ iinfo )
674*
675 ELSE IF( itype.EQ.5 ) THEN
676*
677* Symmetric, eigenvalues specified
678*
679 CALL clatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
680 $ anorm, n, n, 'N', a, lda, work( n+1 ),
681 $ iinfo )
682*
683 ELSE IF( itype.EQ.6 ) THEN
684*
685* General, eigenvalues specified
686*
687 IF( kconds( jtype ).EQ.1 ) THEN
688 conds = one
689 ELSE IF( kconds( jtype ).EQ.2 ) THEN
690 conds = rtulpi
691 ELSE
692 conds = zero
693 END IF
694*
695 CALL clatme( n, 'D', iseed, work, imode, cond, cone,
696 $ 'T', 'T', 'T', rwork, 4, conds, n, n, anorm,
697 $ a, lda, work( 2*n+1 ), iinfo )
698*
699 ELSE IF( itype.EQ.7 ) THEN
700*
701* Diagonal, random eigenvalues
702*
703 CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
704 $ 'T', 'N', work( n+1 ), 1, one,
705 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
706 $ zero, anorm, 'NO', a, lda, idumma, iinfo )
707*
708 ELSE IF( itype.EQ.8 ) THEN
709*
710* Symmetric, random eigenvalues
711*
712 CALL clatmr( n, n, 'D', iseed, 'H', work, 6, one, cone,
713 $ 'T', 'N', work( n+1 ), 1, one,
714 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
715 $ zero, anorm, 'NO', a, lda, idumma, iinfo )
716*
717 ELSE IF( itype.EQ.9 ) THEN
718*
719* General, random eigenvalues
720*
721 CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
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, idumma, iinfo )
725 IF( n.GE.4 ) THEN
726 CALL claset( 'Full', 2, n, czero, czero, a, lda )
727 CALL claset( 'Full', n-3, 1, czero, czero, a( 3, 1 ),
728 $ lda )
729 CALL claset( 'Full', n-3, 2, czero, czero,
730 $ a( 3, n-1 ), lda )
731 CALL claset( 'Full', 1, n, czero, czero, a( n, 1 ),
732 $ lda )
733 END IF
734*
735 ELSE IF( itype.EQ.10 ) THEN
736*
737* Triangular, random eigenvalues
738*
739 CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
740 $ 'T', 'N', work( n+1 ), 1, one,
741 $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
742 $ zero, anorm, 'NO', a, lda, idumma, iinfo )
743*
744 ELSE
745*
746 iinfo = 1
747 END IF
748*
749 IF( iinfo.NE.0 ) THEN
750 WRITE( nounit, fmt = 9991 )'Generator', iinfo, n, jtype,
751 $ ioldsd
752 info = abs( iinfo )
753 RETURN
754 END IF
755*
756 90 CONTINUE
757*
758* Test for minimal and generous workspace
759*
760 DO 120 iwk = 1, 2
761 IF( iwk.EQ.1 ) THEN
762 nnwork = 2*n
763 ELSE
764 nnwork = max( 2*n, n*( n+1 ) / 2 )
765 END IF
766 nnwork = max( nnwork, 1 )
767*
768 CALL cget24( .false., jtype, thresh, ioldsd, nounit, n,
769 $ a, lda, h, ht, w, wt, wtmp, vs, ldvs, vs1,
770 $ rcdein, rcdvin, nslct, islct, 0, result,
771 $ work, nnwork, rwork, bwork, info )
772*
773* Check for RESULT(j) > THRESH
774*
775 ntest = 0
776 nfail = 0
777 DO 100 j = 1, 15
778 IF( result( j ).GE.zero )
779 $ ntest = ntest + 1
780 IF( result( j ).GE.thresh )
781 $ nfail = nfail + 1
782 100 CONTINUE
783*
784 IF( nfail.GT.0 )
785 $ ntestf = ntestf + 1
786 IF( ntestf.EQ.1 ) THEN
787 WRITE( nounit, fmt = 9999 )path
788 WRITE( nounit, fmt = 9998 )
789 WRITE( nounit, fmt = 9997 )
790 WRITE( nounit, fmt = 9996 )
791 WRITE( nounit, fmt = 9995 )thresh
792 WRITE( nounit, fmt = 9994 )
793 ntestf = 2
794 END IF
795*
796 DO 110 j = 1, 15
797 IF( result( j ).GE.thresh ) THEN
798 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
799 $ j, result( j )
800 END IF
801 110 CONTINUE
802*
803 nerrs = nerrs + nfail
804 ntestt = ntestt + ntest
805*
806 120 CONTINUE
807 130 CONTINUE
808 140 CONTINUE
809*
810 150 CONTINUE
811*
812* Read in data from file to check accuracy of condition estimation
813* Read input data until N=0
814*
815 jtype = 0
816 160 CONTINUE
817 READ( niunit, fmt = *, END = 200 )N, NSLCT, isrt
818 IF( n.EQ.0 )
819 $ GO TO 200
820 jtype = jtype + 1
821 iseed( 1 ) = jtype
822 READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
823 DO 170 i = 1, n
824 READ( niunit, fmt = * )( a( i, j ), j = 1, n )
825 170 CONTINUE
826 READ( niunit, fmt = * )rcdein, rcdvin
827*
828 CALL cget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
829 $ w, wt, wtmp, vs, ldvs, vs1, rcdein, rcdvin, nslct,
830 $ islct, isrt, result, work, lwork, rwork, bwork,
831 $ info )
832*
833* Check for RESULT(j) > THRESH
834*
835 ntest = 0
836 nfail = 0
837 DO 180 j = 1, 17
838 IF( result( j ).GE.zero )
839 $ ntest = ntest + 1
840 IF( result( j ).GE.thresh )
841 $ nfail = nfail + 1
842 180 CONTINUE
843*
844 IF( nfail.GT.0 )
845 $ ntestf = ntestf + 1
846 IF( ntestf.EQ.1 ) THEN
847 WRITE( nounit, fmt = 9999 )path
848 WRITE( nounit, fmt = 9998 )
849 WRITE( nounit, fmt = 9997 )
850 WRITE( nounit, fmt = 9996 )
851 WRITE( nounit, fmt = 9995 )thresh
852 WRITE( nounit, fmt = 9994 )
853 ntestf = 2
854 END IF
855 DO 190 j = 1, 17
856 IF( result( j ).GE.thresh ) THEN
857 WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
858 END IF
859 190 CONTINUE
860*
861 nerrs = nerrs + nfail
862 ntestt = ntestt + ntest
863 GO TO 160
864 200 CONTINUE
865*
866* Summary
867*
868 CALL slasum( path, nounit, nerrs, ntestt )
869*
870 9999 FORMAT( / 1x, a3, ' -- Complex Schur Form Decomposition Expert ',
871 $ 'Driver', / ' Matrix types (see CDRVSX for details): ' )
872*
873 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
874 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
875 $ / ' 2=Identity matrix. ', ' 6=Diagona',
876 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
877 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
878 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
879 $ 'mall, evenly spaced.' )
880 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
881 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
882 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
883 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
884 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
885 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
886 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
887 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
888 $ ' complx ' )
889 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
890 $ 'with small random entries.', / ' 20=Matrix with large ran',
891 $ 'dom entries. ', / )
892 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
893 $ / ' ( A denotes A on input and T denotes A on output)',
894 $ / / ' 1 = 0 if T in Schur form (no sort), ',
895 $ ' 1/ulp otherwise', /
896 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
897 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
898 $ / ' 4 = 0 if W are eigenvalues of T (no sort),',
899 $ ' 1/ulp otherwise', /
900 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
901 $ ' 1/ulp otherwise', /
902 $ ' 6 = 0 if W same no matter if VS computed (no sort)',
903 $ ', 1/ulp otherwise' )
904 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
905 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
906 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
907 $ / ' 10 = 0 if W are eigenvalues of T (sort),',
908 $ ' 1/ulp otherwise', /
909 $ ' 11 = 0 if T same no matter what else computed (sort),',
910 $ ' 1/ulp otherwise', /
911 $ ' 12 = 0 if W same no matter what else computed ',
912 $ '(sort), 1/ulp otherwise', /
913 $ ' 13 = 0 if sorting successful, 1/ulp otherwise',
914 $ / ' 14 = 0 if RCONDE same no matter what else computed,',
915 $ ' 1/ulp otherwise', /
916 $ ' 15 = 0 if RCONDv same no matter what else computed,',
917 $ ' 1/ulp otherwise', /
918 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
919 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
920 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
921 $ ' type ', i2, ', test(', i2, ')=', g10.3 )
922 9992 FORMAT( ' N=', i5, ', input example =', i3, ', test(', i2, ')=',
923 $ g10.3 )
924 9991 FORMAT( ' CDRVSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
925 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
926*
927 RETURN
928*
929* End of CDRVSX
930*
931 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cdrvsx(nsizes, nn, ntypes, dotype, iseed, thresh, niunit, nounit, a, lda, h, ht, w, wt, wtmp, vs, ldvs, vs1, result, work, lwork, rwork, bwork, info)
CDRVSX
Definition cdrvsx.f:435
subroutine cget24(comp, jtype, thresh, iseed, nounit, n, a, lda, h, ht, w, wt, wtmp, vs, ldvs, vs1, rcdein, rcdvin, nslct, islct, isrt, result, work, lwork, rwork, bwork, info)
CGET24
Definition cget24.f:335
subroutine clatme(n, dist, iseed, d, mode, cond, dmax, rsign, upper, sim, ds, modes, conds, kl, ku, anorm, a, lda, work, info)
CLATME
Definition clatme.f:301
subroutine clatmr(m, n, dist, iseed, sym, d, mode, cond, dmax, rsign, grade, dl, model, condl, dr, moder, condr, pivtng, ipivot, kl, ku, sparse, anorm, pack, a, lda, iwork, info)
CLATMR
Definition clatmr.f:490
subroutine clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
Definition clatms.f:332
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
subroutine slasum(type, iounit, ie, nrun)
SLASUM
Definition slasum.f:41