LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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 ..
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 ..
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, slabad, 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*
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 )
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 CALL slabad( unfl, ovfl )
587 ulp = slamch( 'Precision' )
588 ulpinv = one / ulp
589 rtulp = sqrt( ulp )
590 rtulpi = one / rtulp
591*
592* Loop over sizes, types
593*
594 nerrs = 0
595*
596 DO 140 jsize = 1, nsizes
597 n = nn( jsize )
598 IF( nsizes.NE.1 ) THEN
599 mtypes = min( maxtyp, ntypes )
600 ELSE
601 mtypes = min( maxtyp+1, ntypes )
602 END IF
603*
604 DO 130 jtype = 1, mtypes
605 IF( .NOT.dotype( jtype ) )
606 \$ GO TO 130
607*
608* Save ISEED in case of an error.
609*
610 DO 20 j = 1, 4
611 ioldsd( j ) = iseed( j )
612 20 CONTINUE
613*
614* Compute "A"
615*
616* Control parameters:
617*
618* KMAGN KCONDS KMODE KTYPE
619* =1 O(1) 1 clustered 1 zero
620* =2 large large clustered 2 identity
621* =3 small exponential Jordan
622* =4 arithmetic diagonal, (w/ eigenvalues)
623* =5 random log symmetric, w/ eigenvalues
624* =6 random general, w/ eigenvalues
625* =7 random diagonal
626* =8 random symmetric
627* =9 random general
628* =10 random triangular
629*
630 IF( mtypes.GT.maxtyp )
631 \$ GO TO 90
632*
633 itype = ktype( jtype )
634 imode = kmode( jtype )
635*
636* Compute norm
637*
638 GO TO ( 30, 40, 50 )kmagn( jtype )
639*
640 30 CONTINUE
641 anorm = one
642 GO TO 60
643*
644 40 CONTINUE
645 anorm = ovfl*ulp
646 GO TO 60
647*
648 50 CONTINUE
649 anorm = unfl*ulpinv
650 GO TO 60
651*
652 60 CONTINUE
653*
654 CALL slaset( 'Full', lda, n, zero, zero, a, lda )
655 iinfo = 0
656 cond = ulpinv
657*
658* Special Matrices -- Identity & Jordan block
659*
660* Zero
661*
662 IF( itype.EQ.1 ) THEN
663 iinfo = 0
664*
665 ELSE IF( itype.EQ.2 ) THEN
666*
667* Identity
668*
669 DO 70 jcol = 1, n
670 a( jcol, jcol ) = anorm
671 70 CONTINUE
672*
673 ELSE IF( itype.EQ.3 ) THEN
674*
675* Jordan Block
676*
677 DO 80 jcol = 1, n
678 a( jcol, jcol ) = anorm
679 IF( jcol.GT.1 )
680 \$ a( jcol, jcol-1 ) = one
681 80 CONTINUE
682*
683 ELSE IF( itype.EQ.4 ) THEN
684*
685* Diagonal Matrix, [Eigen]values Specified
686*
687 CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
688 \$ anorm, 0, 0, 'N', a, lda, work( n+1 ),
689 \$ iinfo )
690*
691 ELSE IF( itype.EQ.5 ) THEN
692*
693* Symmetric, eigenvalues specified
694*
695 CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
696 \$ anorm, n, n, 'N', a, lda, work( n+1 ),
697 \$ iinfo )
698*
699 ELSE IF( itype.EQ.6 ) THEN
700*
701* General, eigenvalues specified
702*
703 IF( kconds( jtype ).EQ.1 ) THEN
704 conds = one
705 ELSE IF( kconds( jtype ).EQ.2 ) THEN
706 conds = rtulpi
707 ELSE
708 conds = zero
709 END IF
710*
711 adumma( 1 ) = ' '
712 CALL slatme( n, 'S', iseed, work, imode, cond, one,
713 \$ adumma, 'T', 'T', 'T', work( n+1 ), 4,
714 \$ conds, n, n, anorm, a, lda, work( 2*n+1 ),
715 \$ iinfo )
716*
717 ELSE IF( itype.EQ.7 ) THEN
718*
719* Diagonal, random eigenvalues
720*
721 CALL slatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
722 \$ 'T', 'N', work( n+1 ), 1, one,
723 \$ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
724 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
725*
726 ELSE IF( itype.EQ.8 ) THEN
727*
728* Symmetric, random eigenvalues
729*
730 CALL slatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
731 \$ 'T', 'N', work( n+1 ), 1, one,
732 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
733 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
734*
735 ELSE IF( itype.EQ.9 ) THEN
736*
737* General, random eigenvalues
738*
739 CALL slatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
740 \$ 'T', 'N', work( n+1 ), 1, one,
741 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
742 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
743 IF( n.GE.4 ) THEN
744 CALL slaset( 'Full', 2, n, zero, zero, a, lda )
745 CALL slaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
746 \$ lda )
747 CALL slaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
748 \$ lda )
749 CALL slaset( 'Full', 1, n, zero, zero, a( n, 1 ),
750 \$ lda )
751 END IF
752*
753 ELSE IF( itype.EQ.10 ) THEN
754*
755* Triangular, random eigenvalues
756*
757 CALL slatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
758 \$ 'T', 'N', work( n+1 ), 1, one,
759 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
760 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
761*
762 ELSE
763*
764 iinfo = 1
765 END IF
766*
767 IF( iinfo.NE.0 ) THEN
768 WRITE( nounit, fmt = 9991 )'Generator', iinfo, n, jtype,
769 \$ ioldsd
770 info = abs( iinfo )
771 RETURN
772 END IF
773*
774 90 CONTINUE
775*
776* Test for minimal and generous workspace
777*
778 DO 120 iwk = 1, 2
779 IF( iwk.EQ.1 ) THEN
780 nnwork = 3*n
781 ELSE
782 nnwork = max( 3*n, 2*n*n )
783 END IF
784 nnwork = max( nnwork, 1 )
785*
786 CALL sget24( .false., jtype, thresh, ioldsd, nounit, n,
787 \$ a, lda, h, ht, wr, wi, wrt, wit, wrtmp,
788 \$ witmp, vs, ldvs, vs1, rcdein, rcdvin, nslct,
789 \$ islct, result, work, nnwork, iwork, bwork,
790 \$ info )
791*
792* Check for RESULT(j) > THRESH
793*
794 ntest = 0
795 nfail = 0
796 DO 100 j = 1, 15
797 IF( result( j ).GE.zero )
798 \$ ntest = ntest + 1
799 IF( result( j ).GE.thresh )
800 \$ nfail = nfail + 1
801 100 CONTINUE
802*
803 IF( nfail.GT.0 )
804 \$ ntestf = ntestf + 1
805 IF( ntestf.EQ.1 ) THEN
806 WRITE( nounit, fmt = 9999 )path
807 WRITE( nounit, fmt = 9998 )
808 WRITE( nounit, fmt = 9997 )
809 WRITE( nounit, fmt = 9996 )
810 WRITE( nounit, fmt = 9995 )thresh
811 WRITE( nounit, fmt = 9994 )
812 ntestf = 2
813 END IF
814*
815 DO 110 j = 1, 15
816 IF( result( j ).GE.thresh ) THEN
817 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
818 \$ j, result( j )
819 END IF
820 110 CONTINUE
821*
822 nerrs = nerrs + nfail
823 ntestt = ntestt + ntest
824*
825 120 CONTINUE
826 130 CONTINUE
827 140 CONTINUE
828*
829 150 CONTINUE
830*
831* Read in data from file to check accuracy of condition estimation
832* Read input data until N=0
833*
834 jtype = 0
835 160 CONTINUE
836 READ( niunit, fmt = *, END = 200 )N, nslct
837 IF( n.EQ.0 )
838 \$ GO TO 200
839 jtype = jtype + 1
840 iseed( 1 ) = jtype
841 IF( nslct.GT.0 )
842 \$ READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
843 DO 170 i = 1, n
844 READ( niunit, fmt = * )( a( i, j ), j = 1, n )
845 170 CONTINUE
846 READ( niunit, fmt = * )rcdein, rcdvin
847*
848 CALL sget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
849 \$ wr, wi, wrt, wit, wrtmp, witmp, vs, ldvs, vs1,
850 \$ rcdein, rcdvin, nslct, islct, result, work, lwork,
851 \$ iwork, bwork, info )
852*
853* Check for RESULT(j) > THRESH
854*
855 ntest = 0
856 nfail = 0
857 DO 180 j = 1, 17
858 IF( result( j ).GE.zero )
859 \$ ntest = ntest + 1
860 IF( result( j ).GE.thresh )
861 \$ nfail = nfail + 1
862 180 CONTINUE
863*
864 IF( nfail.GT.0 )
865 \$ ntestf = ntestf + 1
866 IF( ntestf.EQ.1 ) THEN
867 WRITE( nounit, fmt = 9999 )path
868 WRITE( nounit, fmt = 9998 )
869 WRITE( nounit, fmt = 9997 )
870 WRITE( nounit, fmt = 9996 )
871 WRITE( nounit, fmt = 9995 )thresh
872 WRITE( nounit, fmt = 9994 )
873 ntestf = 2
874 END IF
875 DO 190 j = 1, 17
876 IF( result( j ).GE.thresh ) THEN
877 WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
878 END IF
879 190 CONTINUE
880*
881 nerrs = nerrs + nfail
882 ntestt = ntestt + ntest
883 GO TO 160
884 200 CONTINUE
885*
886* Summary
887*
888 CALL slasum( path, nounit, nerrs, ntestt )
889*
890 9999 FORMAT( / 1x, a3, ' -- Real Schur Form Decomposition Expert ',
891 \$ 'Driver', / ' Matrix types (see SDRVSX for details):' )
892*
893 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
894 \$ ' ', ' 5=Diagonal: geometr. spaced entries.',
895 \$ / ' 2=Identity matrix. ', ' 6=Diagona',
896 \$ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
897 \$ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
898 \$ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
899 \$ 'mall, evenly spaced.' )
900 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
901 \$ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
902 \$ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
903 \$ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
904 \$ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
905 \$ 'lex ', / ' 12=Well-cond., random complex ', ' ',
906 \$ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
907 \$ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
908 \$ ' complx ' )
909 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
910 \$ 'with small random entries.', / ' 20=Matrix with large ran',
911 \$ 'dom entries. ', / )
912 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
913 \$ / ' ( A denotes A on input and T denotes A on output)',
914 \$ / / ' 1 = 0 if T in Schur form (no sort), ',
915 \$ ' 1/ulp otherwise', /
916 \$ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
917 \$ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
918 \$ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
919 \$ ' 1/ulp otherwise', /
920 \$ ' 5 = 0 if T same no matter if VS computed (no sort),',
921 \$ ' 1/ulp otherwise', /
922 \$ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
923 \$ ', 1/ulp otherwise' )
924 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
925 \$ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
926 \$ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
927 \$ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
928 \$ ' 1/ulp otherwise', /
929 \$ ' 11 = 0 if T same no matter what else computed (sort),',
930 \$ ' 1/ulp otherwise', /
931 \$ ' 12 = 0 if WR, WI same no matter what else computed ',
932 \$ '(sort), 1/ulp otherwise', /
933 \$ ' 13 = 0 if sorting successful, 1/ulp otherwise',
934 \$ / ' 14 = 0 if RCONDE same no matter what else computed,',
935 \$ ' 1/ulp otherwise', /
936 \$ ' 15 = 0 if RCONDv same no matter what else computed,',
937 \$ ' 1/ulp otherwise', /
938 \$ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
939 \$ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
940 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
941 \$ ' type ', i2, ', test(', i2, ')=', g10.3 )
942 9992 FORMAT( ' N=', i5, ', input example =', i3, ', test(', i2, ')=',
943 \$ g10.3 )
944 9991 FORMAT( ' SDRVSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
945 \$ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
946*
947 RETURN
948*
949* End of SDRVSX
950*
951 END
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:321
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 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 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 slasum(TYPE, IOUNIT, IE, NRUN)
SLASUM
Definition: slasum.f:41