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