LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zdrvsx.f
Go to the documentation of this file.
1*> \brief \b ZDRVSX
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 ZDRVSX( 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* DOUBLE PRECISION THRESH
20* ..
21* .. Array Arguments ..
22* LOGICAL BWORK( * ), DOTYPE( * )
23* INTEGER ISEED( 4 ), NN( * )
24* DOUBLE PRECISION RESULT( 17 ), RWORK( * )
25* COMPLEX*16 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*> ZDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
37*> expert driver ZGEESX.
38*>
39*> ZDRVSX 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 ZDRVSX 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 ZGEESX 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 ZGEESX 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 ZDRVSX 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 COMPLEX*16 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*16 array, dimension (LDA, max(NN))
311*> Another copy of the test matrix A, modified by ZGEESX.
312*> \endverbatim
313*>
314*> \param[out] HT
315*> \verbatim
316*> HT is COMPLEX*16 array, dimension (LDA, max(NN))
317*> Yet another copy of the test matrix A, modified by ZGEESX.
318*> \endverbatim
319*>
320*> \param[out] W
321*> \verbatim
322*> W is COMPLEX*16 array, dimension (max(NN))
323*> The computed eigenvalues of A.
324*> \endverbatim
325*>
326*> \param[out] WT
327*> \verbatim
328*> WT is COMPLEX*16 array, dimension (max(NN))
329*> Like W, this array contains the eigenvalues of A,
330*> but those computed when ZGEESX only computes a partial
331*> eigendecomposition, i.e. not Schur vectors
332*> \endverbatim
333*>
334*> \param[out] WTMP
335*> \verbatim
336*> WTMP is COMPLEX*16 array, dimension (max(NN))
337*> More temporary storage for eigenvalues.
338*> \endverbatim
339*>
340*> \param[out] VS
341*> \verbatim
342*> VS is COMPLEX*16 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*16 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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, ZLATMR, CLATMS, CLATME or ZGET24 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 complex16_eig
429*
430* =====================================================================
431 SUBROUTINE zdrvsx( 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 DOUBLE PRECISION THRESH
444* ..
445* .. Array Arguments ..
446 LOGICAL BWORK( * ), DOTYPE( * )
447 INTEGER ISEED( 4 ), NN( * )
448 DOUBLE PRECISION RESULT( 17 ), RWORK( * )
449 COMPLEX*16 A( LDA, * ), H( LDA, * ), HT( LDA, * ),
450 \$ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
451 \$ work( * ), wt( * ), wtmp( * )
452* ..
453*
454* =====================================================================
455*
456* .. Parameters ..
457 COMPLEX*16 CZERO
458 PARAMETER ( CZERO = ( 0.0d+0, 0.0d+0 ) )
459 COMPLEX*16 CONE
460 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
461 DOUBLE PRECISION ZERO, ONE
462 parameter( zero = 0.0d+0, one = 1.0d+0 )
463 INTEGER MAXTYP
464 parameter( maxtyp = 21 )
465* ..
466* .. Local Scalars ..
468 CHARACTER*3 PATH
469 INTEGER I, IINFO, IMODE, ISRT, ITYPE, IWK, J, JCOL,
470 \$ jsize, jtype, mtypes, n, nerrs, nfail, nmax,
471 \$ nnwork, nslct, ntest, ntestf, ntestt
472 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH
492 EXTERNAL DLAMCH
493* ..
494* .. External Subroutines ..
495 EXTERNAL dlabad, dlasum, xerbla, zget24, zlaset, zlatme,
496 \$ zlatmr, zlatms
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 ) = 'Zomplex 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( 'ZDRVSX', -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 = dlamch( 'Safe minimum' )
569 ovfl = one / unfl
570 CALL dlabad( unfl, ovfl )
571 ulp = dlamch( 'Precision' )
572 ulpinv = one / ulp
573 rtulp = sqrt( ulp )
574 rtulpi = one / rtulp
575*
576* Loop over sizes, types
577*
578 nerrs = 0
579*
580 DO 140 jsize = 1, nsizes
581 n = nn( jsize )
582 IF( nsizes.NE.1 ) THEN
583 mtypes = min( maxtyp, ntypes )
584 ELSE
585 mtypes = min( maxtyp+1, ntypes )
586 END IF
587*
588 DO 130 jtype = 1, mtypes
589 IF( .NOT.dotype( jtype ) )
590 \$ GO TO 130
591*
592* Save ISEED in case of an error.
593*
594 DO 20 j = 1, 4
595 ioldsd( j ) = iseed( j )
596 20 CONTINUE
597*
598* Compute "A"
599*
600* Control parameters:
601*
602* KMAGN KCONDS KMODE KTYPE
603* =1 O(1) 1 clustered 1 zero
604* =2 large large clustered 2 identity
605* =3 small exponential Jordan
606* =4 arithmetic diagonal, (w/ eigenvalues)
607* =5 random log symmetric, w/ eigenvalues
608* =6 random general, w/ eigenvalues
609* =7 random diagonal
610* =8 random symmetric
611* =9 random general
612* =10 random triangular
613*
614 IF( mtypes.GT.maxtyp )
615 \$ GO TO 90
616*
617 itype = ktype( jtype )
618 imode = kmode( jtype )
619*
620* Compute norm
621*
622 GO TO ( 30, 40, 50 )kmagn( jtype )
623*
624 30 CONTINUE
625 anorm = one
626 GO TO 60
627*
628 40 CONTINUE
629 anorm = ovfl*ulp
630 GO TO 60
631*
632 50 CONTINUE
633 anorm = unfl*ulpinv
634 GO TO 60
635*
636 60 CONTINUE
637*
638 CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
639 iinfo = 0
640 cond = ulpinv
641*
642* Special Matrices -- Identity & Jordan block
643*
644 IF( itype.EQ.1 ) THEN
645*
646* Zero
647*
648 iinfo = 0
649*
650 ELSE IF( itype.EQ.2 ) THEN
651*
652* Identity
653*
654 DO 70 jcol = 1, n
655 a( jcol, jcol ) = anorm
656 70 CONTINUE
657*
658 ELSE IF( itype.EQ.3 ) THEN
659*
660* Jordan Block
661*
662 DO 80 jcol = 1, n
663 a( jcol, jcol ) = anorm
664 IF( jcol.GT.1 )
665 \$ a( jcol, jcol-1 ) = cone
666 80 CONTINUE
667*
668 ELSE IF( itype.EQ.4 ) THEN
669*
670* Diagonal Matrix, [Eigen]values Specified
671*
672 CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
673 \$ anorm, 0, 0, 'N', a, lda, work( n+1 ),
674 \$ iinfo )
675*
676 ELSE IF( itype.EQ.5 ) THEN
677*
678* Symmetric, eigenvalues specified
679*
680 CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode, cond,
681 \$ anorm, n, n, 'N', a, lda, work( n+1 ),
682 \$ iinfo )
683*
684 ELSE IF( itype.EQ.6 ) THEN
685*
686* General, eigenvalues specified
687*
688 IF( kconds( jtype ).EQ.1 ) THEN
689 conds = one
690 ELSE IF( kconds( jtype ).EQ.2 ) THEN
691 conds = rtulpi
692 ELSE
693 conds = zero
694 END IF
695*
696 CALL zlatme( n, 'D', iseed, work, imode, cond, cone,
697 \$ 'T', 'T', 'T', rwork, 4, conds, n, n, anorm,
698 \$ a, lda, work( 2*n+1 ), iinfo )
699*
700 ELSE IF( itype.EQ.7 ) THEN
701*
702* Diagonal, random eigenvalues
703*
704 CALL zlatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
705 \$ 'T', 'N', work( n+1 ), 1, one,
706 \$ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
707 \$ zero, anorm, 'NO', a, lda, idumma, iinfo )
708*
709 ELSE IF( itype.EQ.8 ) THEN
710*
711* Symmetric, random eigenvalues
712*
713 CALL zlatmr( n, n, 'D', iseed, 'H', work, 6, one, cone,
714 \$ 'T', 'N', work( n+1 ), 1, one,
715 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
716 \$ zero, anorm, 'NO', a, lda, idumma, iinfo )
717*
718 ELSE IF( itype.EQ.9 ) THEN
719*
720* General, random eigenvalues
721*
722 CALL zlatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
723 \$ 'T', 'N', work( n+1 ), 1, one,
724 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
725 \$ zero, anorm, 'NO', a, lda, idumma, iinfo )
726 IF( n.GE.4 ) THEN
727 CALL zlaset( 'Full', 2, n, czero, czero, a, lda )
728 CALL zlaset( 'Full', n-3, 1, czero, czero, a( 3, 1 ),
729 \$ lda )
730 CALL zlaset( 'Full', n-3, 2, czero, czero,
731 \$ a( 3, n-1 ), lda )
732 CALL zlaset( 'Full', 1, n, czero, czero, a( n, 1 ),
733 \$ lda )
734 END IF
735*
736 ELSE IF( itype.EQ.10 ) THEN
737*
738* Triangular, random eigenvalues
739*
740 CALL zlatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
741 \$ 'T', 'N', work( n+1 ), 1, one,
742 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
743 \$ zero, anorm, 'NO', a, lda, idumma, iinfo )
744*
745 ELSE
746*
747 iinfo = 1
748 END IF
749*
750 IF( iinfo.NE.0 ) THEN
751 WRITE( nounit, fmt = 9991 )'Generator', iinfo, n, jtype,
752 \$ ioldsd
753 info = abs( iinfo )
754 RETURN
755 END IF
756*
757 90 CONTINUE
758*
759* Test for minimal and generous workspace
760*
761 DO 120 iwk = 1, 2
762 IF( iwk.EQ.1 ) THEN
763 nnwork = 2*n
764 ELSE
765 nnwork = max( 2*n, n*( n+1 ) / 2 )
766 END IF
767 nnwork = max( nnwork, 1 )
768*
769 CALL zget24( .false., jtype, thresh, ioldsd, nounit, n,
770 \$ a, lda, h, ht, w, wt, wtmp, vs, ldvs, vs1,
771 \$ rcdein, rcdvin, nslct, islct, 0, result,
772 \$ work, nnwork, rwork, bwork, info )
773*
774* Check for RESULT(j) > THRESH
775*
776 ntest = 0
777 nfail = 0
778 DO 100 j = 1, 15
779 IF( result( j ).GE.zero )
780 \$ ntest = ntest + 1
781 IF( result( j ).GE.thresh )
782 \$ nfail = nfail + 1
783 100 CONTINUE
784*
785 IF( nfail.GT.0 )
786 \$ ntestf = ntestf + 1
787 IF( ntestf.EQ.1 ) THEN
788 WRITE( nounit, fmt = 9999 )path
789 WRITE( nounit, fmt = 9998 )
790 WRITE( nounit, fmt = 9997 )
791 WRITE( nounit, fmt = 9996 )
792 WRITE( nounit, fmt = 9995 )thresh
793 WRITE( nounit, fmt = 9994 )
794 ntestf = 2
795 END IF
796*
797 DO 110 j = 1, 15
798 IF( result( j ).GE.thresh ) THEN
799 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
800 \$ j, result( j )
801 END IF
802 110 CONTINUE
803*
804 nerrs = nerrs + nfail
805 ntestt = ntestt + ntest
806*
807 120 CONTINUE
808 130 CONTINUE
809 140 CONTINUE
810*
811 150 CONTINUE
812*
813* Read in data from file to check accuracy of condition estimation
814* Read input data until N=0
815*
816 jtype = 0
817 160 CONTINUE
818 READ( niunit, fmt = *, END = 200 )N, NSLCT, isrt
819 IF( n.EQ.0 )
820 \$ GO TO 200
821 jtype = jtype + 1
822 iseed( 1 ) = jtype
823 READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
824 DO 170 i = 1, n
825 READ( niunit, fmt = * )( a( i, j ), j = 1, n )
826 170 CONTINUE
827 READ( niunit, fmt = * )rcdein, rcdvin
828*
829 CALL zget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
830 \$ w, wt, wtmp, vs, ldvs, vs1, rcdein, rcdvin, nslct,
831 \$ islct, isrt, result, work, lwork, rwork, bwork,
832 \$ info )
833*
834* Check for RESULT(j) > THRESH
835*
836 ntest = 0
837 nfail = 0
838 DO 180 j = 1, 17
839 IF( result( j ).GE.zero )
840 \$ ntest = ntest + 1
841 IF( result( j ).GE.thresh )
842 \$ nfail = nfail + 1
843 180 CONTINUE
844*
845 IF( nfail.GT.0 )
846 \$ ntestf = ntestf + 1
847 IF( ntestf.EQ.1 ) THEN
848 WRITE( nounit, fmt = 9999 )path
849 WRITE( nounit, fmt = 9998 )
850 WRITE( nounit, fmt = 9997 )
851 WRITE( nounit, fmt = 9996 )
852 WRITE( nounit, fmt = 9995 )thresh
853 WRITE( nounit, fmt = 9994 )
854 ntestf = 2
855 END IF
856 DO 190 j = 1, 17
857 IF( result( j ).GE.thresh ) THEN
858 WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
859 END IF
860 190 CONTINUE
861*
862 nerrs = nerrs + nfail
863 ntestt = ntestt + ntest
864 GO TO 160
865 200 CONTINUE
866*
867* Summary
868*
869 CALL dlasum( path, nounit, nerrs, ntestt )
870*
871 9999 FORMAT( / 1x, a3, ' -- Complex Schur Form Decomposition Expert ',
872 \$ 'Driver', / ' Matrix types (see ZDRVSX for details): ' )
873*
874 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
875 \$ ' ', ' 5=Diagonal: geometr. spaced entries.',
876 \$ / ' 2=Identity matrix. ', ' 6=Diagona',
877 \$ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
878 \$ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
879 \$ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
880 \$ 'mall, evenly spaced.' )
881 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
882 \$ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
883 \$ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
884 \$ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
885 \$ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
886 \$ 'lex ', / ' 12=Well-cond., random complex ', ' ',
887 \$ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
888 \$ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
889 \$ ' complx ' )
890 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
891 \$ 'with small random entries.', / ' 20=Matrix with large ran',
892 \$ 'dom entries. ', / )
893 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
894 \$ / ' ( A denotes A on input and T denotes A on output)',
895 \$ / / ' 1 = 0 if T in Schur form (no sort), ',
896 \$ ' 1/ulp otherwise', /
897 \$ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
898 \$ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
899 \$ / ' 4 = 0 if W are eigenvalues of T (no sort),',
900 \$ ' 1/ulp otherwise', /
901 \$ ' 5 = 0 if T same no matter if VS computed (no sort),',
902 \$ ' 1/ulp otherwise', /
903 \$ ' 6 = 0 if W same no matter if VS computed (no sort)',
904 \$ ', 1/ulp otherwise' )
905 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
906 \$ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
907 \$ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
908 \$ / ' 10 = 0 if W are eigenvalues of T (sort),',
909 \$ ' 1/ulp otherwise', /
910 \$ ' 11 = 0 if T same no matter what else computed (sort),',
911 \$ ' 1/ulp otherwise', /
912 \$ ' 12 = 0 if W same no matter what else computed ',
913 \$ '(sort), 1/ulp otherwise', /
914 \$ ' 13 = 0 if sorting successful, 1/ulp otherwise',
915 \$ / ' 14 = 0 if RCONDE same no matter what else computed,',
916 \$ ' 1/ulp otherwise', /
917 \$ ' 15 = 0 if RCONDv same no matter what else computed,',
918 \$ ' 1/ulp otherwise', /
919 \$ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
920 \$ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
921 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
922 \$ ' type ', i2, ', test(', i2, ')=', g10.3 )
923 9992 FORMAT( ' N=', i5, ', input example =', i3, ', test(', i2, ')=',
924 \$ g10.3 )
925 9991 FORMAT( ' ZDRVSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
926 \$ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
927*
928 RETURN
929*
930* End of ZDRVSX
931*
932 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zget24(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)
ZGET24
Definition: zget24.f:335
subroutine zdrvsx(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS, LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK, INFO)
ZDRVSX
Definition: zdrvsx.f:435
subroutine zlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
ZLATMS
Definition: zlatms.f:332
subroutine zlatmr(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)
ZLATMR
Definition: zlatmr.f:490
subroutine zlatme(N, DIST, ISEED, D, MODE, COND, DMAX, RSIGN, UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, LDA, WORK, INFO)
ZLATME
Definition: zlatme.f:301
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:43