LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sdrves.f
Go to the documentation of this file.
1*> \brief \b SDRVES
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 SDRVES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12* NOUNIT, A, LDA, H, HT, WR, WI, WRT, WIT, VS,
13* LDVS, RESULT, WORK, NWORK, IWORK, BWORK, INFO )
14*
15* .. Scalar Arguments ..
16* INTEGER INFO, LDA, LDVS, NOUNIT, NSIZES, NTYPES, NWORK
17* REAL THRESH
18* ..
19* .. Array Arguments ..
20* LOGICAL BWORK( * ), DOTYPE( * )
21* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
22* REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ),
23* $ RESULT( 13 ), VS( LDVS, * ), WI( * ), WIT( * ),
24* $ WORK( * ), WR( * ), WRT( * )
25* ..
26*
27*
28*> \par Purpose:
29* =============
30*>
31*> \verbatim
32*>
33*> SDRVES checks the nonsymmetric eigenvalue (Schur form) problem
34*> driver SGEES.
35*>
36*> When SDRVES is called, a number of matrix "sizes" ("n's") and a
37*> number of matrix "types" are specified. For each size ("n")
38*> and each type of matrix, one matrix will be generated and used
39*> to test the nonsymmetric eigenroutines. For each matrix, 13
40*> tests will be performed:
41*>
42*> (1) 0 if T is in Schur form, 1/ulp otherwise
43*> (no sorting of eigenvalues)
44*>
45*> (2) | A - VS T VS' | / ( n |A| ulp )
46*>
47*> Here VS is the matrix of Schur eigenvectors, and T is in Schur
48*> form (no sorting of eigenvalues).
49*>
50*> (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
51*>
52*> (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T
53*> 1/ulp otherwise
54*> (no sorting of eigenvalues)
55*>
56*> (5) 0 if T(with VS) = T(without VS),
57*> 1/ulp otherwise
58*> (no sorting of eigenvalues)
59*>
60*> (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
61*> 1/ulp otherwise
62*> (no sorting of eigenvalues)
63*>
64*> (7) 0 if T is in Schur form, 1/ulp otherwise
65*> (with sorting of eigenvalues)
66*>
67*> (8) | A - VS T VS' | / ( n |A| ulp )
68*>
69*> Here VS is the matrix of Schur eigenvectors, and T is in Schur
70*> form (with sorting of eigenvalues).
71*>
72*> (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
73*>
74*> (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T
75*> 1/ulp otherwise
76*> (with sorting of eigenvalues)
77*>
78*> (11) 0 if T(with VS) = T(without VS),
79*> 1/ulp otherwise
80*> (with sorting of eigenvalues)
81*>
82*> (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
83*> 1/ulp otherwise
84*> (with sorting of eigenvalues)
85*>
86*> (13) if sorting worked and SDIM is the number of
87*> eigenvalues which were SELECTed
88*>
89*> The "sizes" are specified by an array NN(1:NSIZES); the value of
90*> each element NN(j) specifies one size.
91*> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
92*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
93*> Currently, the list of possible types is:
94*>
95*> (1) The zero matrix.
96*> (2) The identity matrix.
97*> (3) A (transposed) Jordan block, with 1's on the diagonal.
98*>
99*> (4) A diagonal matrix with evenly spaced entries
100*> 1, ..., ULP and random signs.
101*> (ULP = (first number larger than 1) - 1 )
102*> (5) A diagonal matrix with geometrically spaced entries
103*> 1, ..., ULP and random signs.
104*> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
105*> and random signs.
106*>
107*> (7) Same as (4), but multiplied by a constant near
108*> the overflow threshold
109*> (8) Same as (4), but multiplied by a constant near
110*> the underflow threshold
111*>
112*> (9) A matrix of the form U' T U, where U is orthogonal and
113*> T has evenly spaced entries 1, ..., ULP with random signs
114*> on the diagonal and random O(1) entries in the upper
115*> triangle.
116*>
117*> (10) A matrix of the form U' T U, where U is orthogonal and
118*> T has geometrically spaced entries 1, ..., ULP with random
119*> signs on the diagonal and random O(1) entries in the upper
120*> triangle.
121*>
122*> (11) A matrix of the form U' T U, where U is orthogonal and
123*> T has "clustered" entries 1, ULP,..., ULP with random
124*> signs on the diagonal and random O(1) entries in the upper
125*> triangle.
126*>
127*> (12) A matrix of the form U' T U, where U is orthogonal and
128*> T has real or complex conjugate paired eigenvalues randomly
129*> chosen from ( ULP, 1 ) and random O(1) entries in the upper
130*> triangle.
131*>
132*> (13) A matrix of the form X' T X, where X has condition
133*> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
134*> with random signs on the diagonal and random O(1) entries
135*> in the upper triangle.
136*>
137*> (14) A matrix of the form X' T X, where X has condition
138*> SQRT( ULP ) and T has geometrically spaced entries
139*> 1, ..., ULP with random signs on the diagonal and random
140*> O(1) entries in the upper triangle.
141*>
142*> (15) A matrix of the form X' T X, where X has condition
143*> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
144*> with random signs on the diagonal and random O(1) entries
145*> in the upper triangle.
146*>
147*> (16) A matrix of the form X' T X, where X has condition
148*> SQRT( ULP ) and T has real or complex conjugate paired
149*> eigenvalues randomly chosen from ( ULP, 1 ) and random
150*> O(1) entries in the upper triangle.
151*>
152*> (17) Same as (16), but multiplied by a constant
153*> near the overflow threshold
154*> (18) Same as (16), but multiplied by a constant
155*> near the underflow threshold
156*>
157*> (19) Nonsymmetric matrix with random entries chosen from (-1,1).
158*> If N is at least 4, all entries in first two rows and last
159*> row, and first column and last two columns are zero.
160*> (20) Same as (19), but multiplied by a constant
161*> near the overflow threshold
162*> (21) Same as (19), but multiplied by a constant
163*> near the underflow threshold
164*> \endverbatim
165*
166* Arguments:
167* ==========
168*
169*> \param[in] NSIZES
170*> \verbatim
171*> NSIZES is INTEGER
172*> The number of sizes of matrices to use. If it is zero,
173*> SDRVES does nothing. It must be at least zero.
174*> \endverbatim
175*>
176*> \param[in] NN
177*> \verbatim
178*> NN is INTEGER array, dimension (NSIZES)
179*> An array containing the sizes to be used for the matrices.
180*> Zero values will be skipped. The values must be at least
181*> zero.
182*> \endverbatim
183*>
184*> \param[in] NTYPES
185*> \verbatim
186*> NTYPES is INTEGER
187*> The number of elements in DOTYPE. If it is zero, SDRVES
188*> does nothing. It must be at least zero. If it is MAXTYP+1
189*> and NSIZES is 1, then an additional type, MAXTYP+1 is
190*> defined, which is to use whatever matrix is in A. This
191*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
192*> DOTYPE(MAXTYP+1) is .TRUE. .
193*> \endverbatim
194*>
195*> \param[in] DOTYPE
196*> \verbatim
197*> DOTYPE is LOGICAL array, dimension (NTYPES)
198*> If DOTYPE(j) is .TRUE., then for each size in NN a
199*> matrix of that size and of type j will be generated.
200*> If NTYPES is smaller than the maximum number of types
201*> defined (PARAMETER MAXTYP), then types NTYPES+1 through
202*> MAXTYP will not be generated. If NTYPES is larger
203*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
204*> will be ignored.
205*> \endverbatim
206*>
207*> \param[in,out] ISEED
208*> \verbatim
209*> ISEED is INTEGER array, dimension (4)
210*> On entry ISEED specifies the seed of the random number
211*> generator. The array elements should be between 0 and 4095;
212*> if not they will be reduced mod 4096. Also, ISEED(4) must
213*> be odd. The random number generator uses a linear
214*> congruential sequence limited to small integers, and so
215*> should produce machine independent random numbers. The
216*> values of ISEED are changed on exit, and can be used in the
217*> next call to SDRVES to continue the same random number
218*> sequence.
219*> \endverbatim
220*>
221*> \param[in] THRESH
222*> \verbatim
223*> THRESH is REAL
224*> A test will count as "failed" if the "error", computed as
225*> described above, exceeds THRESH. Note that the error
226*> is scaled to be O(1), so THRESH should be a reasonably
227*> small multiple of 1, e.g., 10 or 100. In particular,
228*> it should not depend on the precision (single vs. double)
229*> or the size of the matrix. It must be at least zero.
230*> \endverbatim
231*>
232*> \param[in] NOUNIT
233*> \verbatim
234*> NOUNIT is INTEGER
235*> The FORTRAN unit number for printing out error messages
236*> (e.g., if a routine returns INFO not equal to 0.)
237*> \endverbatim
238*>
239*> \param[out] A
240*> \verbatim
241*> A is REAL array, dimension (LDA, max(NN))
242*> Used to hold the matrix whose eigenvalues are to be
243*> computed. On exit, A contains the last matrix actually used.
244*> \endverbatim
245*>
246*> \param[in] LDA
247*> \verbatim
248*> LDA is INTEGER
249*> The leading dimension of A, and H. LDA must be at
250*> least 1 and at least max(NN).
251*> \endverbatim
252*>
253*> \param[out] H
254*> \verbatim
255*> H is REAL array, dimension (LDA, max(NN))
256*> Another copy of the test matrix A, modified by SGEES.
257*> \endverbatim
258*>
259*> \param[out] HT
260*> \verbatim
261*> HT is REAL array, dimension (LDA, max(NN))
262*> Yet another copy of the test matrix A, modified by SGEES.
263*> \endverbatim
264*>
265*> \param[out] WR
266*> \verbatim
267*> WR is REAL array, dimension (max(NN))
268*> \endverbatim
269*>
270*> \param[out] WI
271*> \verbatim
272*> WI is REAL array, dimension (max(NN))
273*>
274*> The real and imaginary parts of the eigenvalues of A.
275*> On exit, WR + WI*i are the eigenvalues of the matrix in A.
276*> \endverbatim
277*>
278*> \param[out] WRT
279*> \verbatim
280*> WRT is REAL array, dimension (max(NN))
281*> \endverbatim
282*>
283*> \param[out] WIT
284*> \verbatim
285*> WIT is REAL array, dimension (max(NN))
286*>
287*> Like WR, WI, these arrays contain the eigenvalues of A,
288*> but those computed when SGEES only computes a partial
289*> eigendecomposition, i.e. not Schur vectors
290*> \endverbatim
291*>
292*> \param[out] VS
293*> \verbatim
294*> VS is REAL array, dimension (LDVS, max(NN))
295*> VS holds the computed Schur vectors.
296*> \endverbatim
297*>
298*> \param[in] LDVS
299*> \verbatim
300*> LDVS is INTEGER
301*> Leading dimension of VS. Must be at least max(1,max(NN)).
302*> \endverbatim
303*>
304*> \param[out] RESULT
305*> \verbatim
306*> RESULT is REAL array, dimension (13)
307*> The values computed by the 13 tests described above.
308*> The values are currently limited to 1/ulp, to avoid overflow.
309*> \endverbatim
310*>
311*> \param[out] WORK
312*> \verbatim
313*> WORK is REAL array, dimension (NWORK)
314*> \endverbatim
315*>
316*> \param[in] NWORK
317*> \verbatim
318*> NWORK is INTEGER
319*> The number of entries in WORK. This must be at least
320*> 5*NN(j)+2*NN(j)**2 for all j.
321*> \endverbatim
322*>
323*> \param[out] IWORK
324*> \verbatim
325*> IWORK is INTEGER array, dimension (max(NN))
326*> \endverbatim
327*>
328*> \param[out] BWORK
329*> \verbatim
330*> BWORK is LOGICAL array, dimension (max(NN))
331*> \endverbatim
332*>
333*> \param[out] INFO
334*> \verbatim
335*> INFO is INTEGER
336*> If 0, then everything ran OK.
337*> -1: NSIZES < 0
338*> -2: Some NN(j) < 0
339*> -3: NTYPES < 0
340*> -6: THRESH < 0
341*> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
342*> -17: LDVS < 1 or LDVS < NMAX, where NMAX is max( NN(j) ).
343*> -20: NWORK too small.
344*> If SLATMR, SLATMS, SLATME or SGEES returns an error code,
345*> the absolute value of it is returned.
346*>
347*>-----------------------------------------------------------------------
348*>
349*> Some Local Variables and Parameters:
350*> ---- ----- --------- --- ----------
351*>
352*> ZERO, ONE Real 0 and 1.
353*> MAXTYP The number of types defined.
354*> NMAX Largest value in NN.
355*> NERRS The number of tests which have exceeded THRESH
356*> COND, CONDS,
357*> IMODE Values to be passed to the matrix generators.
358*> ANORM Norm of A; passed to matrix generators.
359*>
360*> OVFL, UNFL Overflow and underflow thresholds.
361*> ULP, ULPINV Finest relative precision and its inverse.
362*> RTULP, RTULPI Square roots of the previous 4 values.
363*>
364*> The following four arrays decode JTYPE:
365*> KTYPE(j) The general type (1-10) for type "j".
366*> KMODE(j) The MODE value to be passed to the matrix
367*> generator for type "j".
368*> KMAGN(j) The order of magnitude ( O(1),
369*> O(overflow^(1/2) ), O(underflow^(1/2) )
370*> KCONDS(j) Selectw whether CONDS is to be 1 or
371*> 1/sqrt(ulp). (0 means irrelevant.)
372*> \endverbatim
373*
374* Authors:
375* ========
376*
377*> \author Univ. of Tennessee
378*> \author Univ. of California Berkeley
379*> \author Univ. of Colorado Denver
380*> \author NAG Ltd.
381*
382*> \ingroup single_eig
383*
384* =====================================================================
385 SUBROUTINE sdrves( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
386 $ NOUNIT, A, LDA, H, HT, WR, WI, WRT, WIT, VS,
387 $ LDVS, RESULT, WORK, NWORK, IWORK, BWORK, INFO )
388*
389* -- LAPACK test routine --
390* -- LAPACK is a software package provided by Univ. of Tennessee, --
391* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
392*
393* .. Scalar Arguments ..
394 INTEGER INFO, LDA, LDVS, NOUNIT, NSIZES, NTYPES, NWORK
395 REAL THRESH
396* ..
397* .. Array Arguments ..
398 LOGICAL BWORK( * ), DOTYPE( * )
399 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
400 REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ),
401 $ result( 13 ), vs( ldvs, * ), wi( * ), wit( * ),
402 $ work( * ), wr( * ), wrt( * )
403* ..
404*
405* =====================================================================
406*
407* .. Parameters ..
408 REAL ZERO, ONE
409 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
410 INTEGER MAXTYP
411 parameter( maxtyp = 21 )
412* ..
413* .. Local Scalars ..
414 LOGICAL BADNN
415 CHARACTER SORT
416 CHARACTER*3 PATH
417 INTEGER I, IINFO, IMODE, ISORT, ITYPE, IWK, J, JCOL,
418 $ jsize, jtype, knteig, lwork, mtypes, n,
419 $ nerrs, nfail, nmax, nnwork, ntest, ntestf,
420 $ ntestt, rsub, sdim
421 REAL ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TMP,
422 $ ULP, ULPINV, UNFL
423* ..
424* .. Local Arrays ..
425 CHARACTER ADUMMA( 1 )
426 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
427 $ kmagn( maxtyp ), kmode( maxtyp ),
428 $ ktype( maxtyp )
429 REAL RES( 2 )
430* ..
431* .. Arrays in Common ..
432 LOGICAL SELVAL( 20 )
433 REAL SELWI( 20 ), SELWR( 20 )
434* ..
435* .. Scalars in Common ..
436 INTEGER SELDIM, SELOPT
437* ..
438* .. Common blocks ..
439 COMMON / sslct / selopt, seldim, selval, selwr, selwi
440* ..
441* .. External Functions ..
442 LOGICAL SSLECT
443 REAL SLAMCH
444 EXTERNAL sslect, slamch
445* ..
446* .. External Subroutines ..
447 EXTERNAL sgees, shst01, slabad, slacpy, slasum, slatme,
449* ..
450* .. Intrinsic Functions ..
451 INTRINSIC abs, max, sign, sqrt
452* ..
453* .. Data statements ..
454 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
455 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
456 $ 3, 1, 2, 3 /
457 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
458 $ 1, 5, 5, 5, 4, 3, 1 /
459 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
460* ..
461* .. Executable Statements ..
462*
463 path( 1: 1 ) = 'Single precision'
464 path( 2: 3 ) = 'ES'
465*
466* Check for errors
467*
468 ntestt = 0
469 ntestf = 0
470 info = 0
471 selopt = 0
472*
473* Important constants
474*
475 badnn = .false.
476 nmax = 0
477 DO 10 j = 1, nsizes
478 nmax = max( nmax, nn( j ) )
479 IF( nn( j ).LT.0 )
480 $ badnn = .true.
481 10 CONTINUE
482*
483* Check for errors
484*
485 IF( nsizes.LT.0 ) THEN
486 info = -1
487 ELSE IF( badnn ) THEN
488 info = -2
489 ELSE IF( ntypes.LT.0 ) THEN
490 info = -3
491 ELSE IF( thresh.LT.zero ) THEN
492 info = -6
493 ELSE IF( nounit.LE.0 ) THEN
494 info = -7
495 ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
496 info = -9
497 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax ) THEN
498 info = -17
499 ELSE IF( 5*nmax+2*nmax**2.GT.nwork ) THEN
500 info = -20
501 END IF
502*
503 IF( info.NE.0 ) THEN
504 CALL xerbla( 'SDRVES', -info )
505 RETURN
506 END IF
507*
508* Quick return if nothing to do
509*
510 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
511 $ RETURN
512*
513* More Important constants
514*
515 unfl = slamch( 'Safe minimum' )
516 ovfl = one / unfl
517 CALL slabad( unfl, ovfl )
518 ulp = slamch( 'Precision' )
519 ulpinv = one / ulp
520 rtulp = sqrt( ulp )
521 rtulpi = one / rtulp
522*
523* Loop over sizes, types
524*
525 nerrs = 0
526*
527 DO 270 jsize = 1, nsizes
528 n = nn( jsize )
529 mtypes = maxtyp
530 IF( nsizes.EQ.1 .AND. ntypes.EQ.maxtyp+1 )
531 $ mtypes = mtypes + 1
532*
533 DO 260 jtype = 1, mtypes
534 IF( .NOT.dotype( jtype ) )
535 $ GO TO 260
536*
537* Save ISEED in case of an error.
538*
539 DO 20 j = 1, 4
540 ioldsd( j ) = iseed( j )
541 20 CONTINUE
542*
543* Compute "A"
544*
545* Control parameters:
546*
547* KMAGN KCONDS KMODE KTYPE
548* =1 O(1) 1 clustered 1 zero
549* =2 large large clustered 2 identity
550* =3 small exponential Jordan
551* =4 arithmetic diagonal, (w/ eigenvalues)
552* =5 random log symmetric, w/ eigenvalues
553* =6 random general, w/ eigenvalues
554* =7 random diagonal
555* =8 random symmetric
556* =9 random general
557* =10 random triangular
558*
559 IF( mtypes.GT.maxtyp )
560 $ GO TO 90
561*
562 itype = ktype( jtype )
563 imode = kmode( jtype )
564*
565* Compute norm
566*
567 GO TO ( 30, 40, 50 )kmagn( jtype )
568*
569 30 CONTINUE
570 anorm = one
571 GO TO 60
572*
573 40 CONTINUE
574 anorm = ovfl*ulp
575 GO TO 60
576*
577 50 CONTINUE
578 anorm = unfl*ulpinv
579 GO TO 60
580*
581 60 CONTINUE
582*
583 CALL slaset( 'Full', lda, n, zero, zero, a, lda )
584 iinfo = 0
585 cond = ulpinv
586*
587* Special Matrices -- Identity & Jordan block
588*
589* Zero
590*
591 IF( itype.EQ.1 ) THEN
592 iinfo = 0
593*
594 ELSE IF( itype.EQ.2 ) THEN
595*
596* Identity
597*
598 DO 70 jcol = 1, n
599 a( jcol, jcol ) = anorm
600 70 CONTINUE
601*
602 ELSE IF( itype.EQ.3 ) THEN
603*
604* Jordan Block
605*
606 DO 80 jcol = 1, n
607 a( jcol, jcol ) = anorm
608 IF( jcol.GT.1 )
609 $ a( jcol, jcol-1 ) = one
610 80 CONTINUE
611*
612 ELSE IF( itype.EQ.4 ) THEN
613*
614* Diagonal Matrix, [Eigen]values Specified
615*
616 CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
617 $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
618 $ iinfo )
619*
620 ELSE IF( itype.EQ.5 ) THEN
621*
622* Symmetric, eigenvalues specified
623*
624 CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
625 $ anorm, n, n, 'N', a, lda, work( n+1 ),
626 $ iinfo )
627*
628 ELSE IF( itype.EQ.6 ) THEN
629*
630* General, eigenvalues specified
631*
632 IF( kconds( jtype ).EQ.1 ) THEN
633 conds = one
634 ELSE IF( kconds( jtype ).EQ.2 ) THEN
635 conds = rtulpi
636 ELSE
637 conds = zero
638 END IF
639*
640 adumma( 1 ) = ' '
641 CALL slatme( n, 'S', iseed, work, imode, cond, one,
642 $ adumma, 'T', 'T', 'T', work( n+1 ), 4,
643 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
644 $ iinfo )
645*
646 ELSE IF( itype.EQ.7 ) THEN
647*
648* Diagonal, random eigenvalues
649*
650 CALL slatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
651 $ 'T', 'N', work( n+1 ), 1, one,
652 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
653 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
654*
655 ELSE IF( itype.EQ.8 ) THEN
656*
657* Symmetric, random eigenvalues
658*
659 CALL slatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
660 $ 'T', 'N', work( n+1 ), 1, one,
661 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
662 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
663*
664 ELSE IF( itype.EQ.9 ) THEN
665*
666* General, random eigenvalues
667*
668 CALL slatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
669 $ 'T', 'N', work( n+1 ), 1, one,
670 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
671 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
672 IF( n.GE.4 ) THEN
673 CALL slaset( 'Full', 2, n, zero, zero, a, lda )
674 CALL slaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
675 $ lda )
676 CALL slaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
677 $ lda )
678 CALL slaset( 'Full', 1, n, zero, zero, a( n, 1 ),
679 $ lda )
680 END IF
681*
682 ELSE IF( itype.EQ.10 ) THEN
683*
684* Triangular, random eigenvalues
685*
686 CALL slatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
687 $ 'T', 'N', work( n+1 ), 1, one,
688 $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
689 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
690*
691 ELSE
692*
693 iinfo = 1
694 END IF
695*
696 IF( iinfo.NE.0 ) THEN
697 WRITE( nounit, fmt = 9992 )'Generator', iinfo, n, jtype,
698 $ ioldsd
699 info = abs( iinfo )
700 RETURN
701 END IF
702*
703 90 CONTINUE
704*
705* Test for minimal and generous workspace
706*
707 DO 250 iwk = 1, 2
708 IF( iwk.EQ.1 ) THEN
709 nnwork = 3*n
710 ELSE
711 nnwork = 5*n + 2*n**2
712 END IF
713 nnwork = max( nnwork, 1 )
714*
715* Initialize RESULT
716*
717 DO 100 j = 1, 13
718 result( j ) = -one
719 100 CONTINUE
720*
721* Test with and without sorting of eigenvalues
722*
723 DO 210 isort = 0, 1
724 IF( isort.EQ.0 ) THEN
725 sort = 'N'
726 rsub = 0
727 ELSE
728 sort = 'S'
729 rsub = 6
730 END IF
731*
732* Compute Schur form and Schur vectors, and test them
733*
734 CALL slacpy( 'F', n, n, a, lda, h, lda )
735 CALL sgees( 'V', sort, sslect, n, h, lda, sdim, wr,
736 $ wi, vs, ldvs, work, nnwork, bwork, iinfo )
737 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
738 result( 1+rsub ) = ulpinv
739 WRITE( nounit, fmt = 9992 )'SGEES1', iinfo, n,
740 $ jtype, ioldsd
741 info = abs( iinfo )
742 GO TO 220
743 END IF
744*
745* Do Test (1) or Test (7)
746*
747 result( 1+rsub ) = zero
748 DO 120 j = 1, n - 2
749 DO 110 i = j + 2, n
750 IF( h( i, j ).NE.zero )
751 $ result( 1+rsub ) = ulpinv
752 110 CONTINUE
753 120 CONTINUE
754 DO 130 i = 1, n - 2
755 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.
756 $ zero )result( 1+rsub ) = ulpinv
757 130 CONTINUE
758 DO 140 i = 1, n - 1
759 IF( h( i+1, i ).NE.zero ) THEN
760 IF( h( i, i ).NE.h( i+1, i+1 ) .OR.
761 $ h( i, i+1 ).EQ.zero .OR.
762 $ sign( one, h( i+1, i ) ).EQ.
763 $ sign( one, h( i, i+1 ) ) )result( 1+rsub )
764 $ = ulpinv
765 END IF
766 140 CONTINUE
767*
768* Do Tests (2) and (3) or Tests (8) and (9)
769*
770 lwork = max( 1, 2*n*n )
771 CALL shst01( n, 1, n, a, lda, h, lda, vs, ldvs, work,
772 $ lwork, res )
773 result( 2+rsub ) = res( 1 )
774 result( 3+rsub ) = res( 2 )
775*
776* Do Test (4) or Test (10)
777*
778 result( 4+rsub ) = zero
779 DO 150 i = 1, n
780 IF( h( i, i ).NE.wr( i ) )
781 $ result( 4+rsub ) = ulpinv
782 150 CONTINUE
783 IF( n.GT.1 ) THEN
784 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
785 $ result( 4+rsub ) = ulpinv
786 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
787 $ result( 4+rsub ) = ulpinv
788 END IF
789 DO 160 i = 1, n - 1
790 IF( h( i+1, i ).NE.zero ) THEN
791 tmp = sqrt( abs( h( i+1, i ) ) )*
792 $ sqrt( abs( h( i, i+1 ) ) )
793 result( 4+rsub ) = max( result( 4+rsub ),
794 $ abs( wi( i )-tmp ) /
795 $ max( ulp*tmp, unfl ) )
796 result( 4+rsub ) = max( result( 4+rsub ),
797 $ abs( wi( i+1 )+tmp ) /
798 $ max( ulp*tmp, unfl ) )
799 ELSE IF( i.GT.1 ) THEN
800 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.
801 $ zero .AND. wi( i ).NE.zero )result( 4+rsub )
802 $ = ulpinv
803 END IF
804 160 CONTINUE
805*
806* Do Test (5) or Test (11)
807*
808 CALL slacpy( 'F', n, n, a, lda, ht, lda )
809 CALL sgees( 'N', sort, sslect, n, ht, lda, sdim, wrt,
810 $ wit, vs, ldvs, work, nnwork, bwork,
811 $ iinfo )
812 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
813 result( 5+rsub ) = ulpinv
814 WRITE( nounit, fmt = 9992 )'SGEES2', iinfo, n,
815 $ jtype, ioldsd
816 info = abs( iinfo )
817 GO TO 220
818 END IF
819*
820 result( 5+rsub ) = zero
821 DO 180 j = 1, n
822 DO 170 i = 1, n
823 IF( h( i, j ).NE.ht( i, j ) )
824 $ result( 5+rsub ) = ulpinv
825 170 CONTINUE
826 180 CONTINUE
827*
828* Do Test (6) or Test (12)
829*
830 result( 6+rsub ) = zero
831 DO 190 i = 1, n
832 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
833 $ result( 6+rsub ) = ulpinv
834 190 CONTINUE
835*
836* Do Test (13)
837*
838 IF( isort.EQ.1 ) THEN
839 result( 13 ) = zero
840 knteig = 0
841 DO 200 i = 1, n
842 IF( sslect( wr( i ), wi( i ) ) .OR.
843 $ sslect( wr( i ), -wi( i ) ) )
844 $ knteig = knteig + 1
845 IF( i.LT.n ) THEN
846 IF( ( sslect( wr( i+1 ),
847 $ wi( i+1 ) ) .OR. sslect( wr( i+1 ),
848 $ -wi( i+1 ) ) ) .AND.
849 $ ( .NOT.( sslect( wr( i ),
850 $ wi( i ) ) .OR. sslect( wr( i ),
851 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )
852 $ result( 13 ) = ulpinv
853 END IF
854 200 CONTINUE
855 IF( sdim.NE.knteig ) THEN
856 result( 13 ) = ulpinv
857 END IF
858 END IF
859*
860 210 CONTINUE
861*
862* End of Loop -- Check for RESULT(j) > THRESH
863*
864 220 CONTINUE
865*
866 ntest = 0
867 nfail = 0
868 DO 230 j = 1, 13
869 IF( result( j ).GE.zero )
870 $ ntest = ntest + 1
871 IF( result( j ).GE.thresh )
872 $ nfail = nfail + 1
873 230 CONTINUE
874*
875 IF( nfail.GT.0 )
876 $ ntestf = ntestf + 1
877 IF( ntestf.EQ.1 ) THEN
878 WRITE( nounit, fmt = 9999 )path
879 WRITE( nounit, fmt = 9998 )
880 WRITE( nounit, fmt = 9997 )
881 WRITE( nounit, fmt = 9996 )
882 WRITE( nounit, fmt = 9995 )thresh
883 WRITE( nounit, fmt = 9994 )
884 ntestf = 2
885 END IF
886*
887 DO 240 j = 1, 13
888 IF( result( j ).GE.thresh ) THEN
889 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
890 $ j, result( j )
891 END IF
892 240 CONTINUE
893*
894 nerrs = nerrs + nfail
895 ntestt = ntestt + ntest
896*
897 250 CONTINUE
898 260 CONTINUE
899 270 CONTINUE
900*
901* Summary
902*
903 CALL slasum( path, nounit, nerrs, ntestt )
904*
905 9999 FORMAT( / 1x, a3, ' -- Real Schur Form Decomposition Driver',
906 $ / ' Matrix types (see SDRVES for details): ' )
907*
908 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
909 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
910 $ / ' 2=Identity matrix. ', ' 6=Diagona',
911 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
912 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
913 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
914 $ 'mall, evenly spaced.' )
915 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
916 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
917 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
918 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
919 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
920 $ 'lex ', / ' 12=Well-cond., random complex ', 6x, ' ',
921 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
922 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
923 $ ' complx ' )
924 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
925 $ 'with small random entries.', / ' 20=Matrix with large ran',
926 $ 'dom entries. ', / )
927 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
928 $ / ' ( A denotes A on input and T denotes A on output)',
929 $ / / ' 1 = 0 if T in Schur form (no sort), ',
930 $ ' 1/ulp otherwise', /
931 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
932 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
933 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
934 $ ' 1/ulp otherwise', /
935 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
936 $ ' 1/ulp otherwise', /
937 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
938 $ ', 1/ulp otherwise' )
939 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
940 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
941 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
942 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
943 $ ' 1/ulp otherwise', /
944 $ ' 11 = 0 if T same no matter if VS computed (sort),',
945 $ ' 1/ulp otherwise', /
946 $ ' 12 = 0 if WR, WI same no matter if VS computed (sort),',
947 $ ' 1/ulp otherwise', /
948 $ ' 13 = 0 if sorting successful, 1/ulp otherwise', / )
949 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
950 $ ' type ', i2, ', test(', i2, ')=', g10.3 )
951 9992 FORMAT( ' SDRVES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
952 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
953*
954 RETURN
955*
956* End of SDRVES
957*
958 END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 sgees(JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI, VS, LDVS, WORK, LWORK, BWORK, INFO)
SGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
Definition: sgees.f:216
subroutine sdrves(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, H, HT, WR, WI, WRT, WIT, VS, LDVS, RESULT, WORK, NWORK, IWORK, BWORK, INFO)
SDRVES
Definition: sdrves.f:388
subroutine shst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
SHST01
Definition: shst01.f:134
subroutine slasum(TYPE, IOUNIT, IE, NRUN)
SLASUM
Definition: slasum.f:41