LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ddrves.f
Go to the documentation of this file.
1*> \brief \b DDRVES
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 DDRVES( 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* DOUBLE PRECISION THRESH
18* ..
19* .. Array Arguments ..
20* LOGICAL BWORK( * ), DOTYPE( * )
21* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
22* DOUBLE PRECISION 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*> DDRVES checks the nonsymmetric eigenvalue (Schur form) problem
34*> driver DGEES.
35*>
36*> When DDRVES 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*> DDRVES 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, DDRVES
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 DDRVES to continue the same random number
218*> sequence.
219*> \endverbatim
220*>
221*> \param[in] THRESH
222*> \verbatim
223*> THRESH is DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDA, max(NN))
256*> Another copy of the test matrix A, modified by DGEES.
257*> \endverbatim
258*>
259*> \param[out] HT
260*> \verbatim
261*> HT is DOUBLE PRECISION array, dimension (LDA, max(NN))
262*> Yet another copy of the test matrix A, modified by DGEES.
263*> \endverbatim
264*>
265*> \param[out] WR
266*> \verbatim
267*> WR is DOUBLE PRECISION array, dimension (max(NN))
268*> \endverbatim
269*>
270*> \param[out] WI
271*> \verbatim
272*> WI is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (max(NN))
281*> \endverbatim
282*>
283*> \param[out] WIT
284*> \verbatim
285*> WIT is DOUBLE PRECISION array, dimension (max(NN))
286*>
287*> Like WR, WI, these arrays contain the eigenvalues of A,
288*> but those computed when DGEES only computes a partial
289*> eigendecomposition, i.e. not Schur vectors
290*> \endverbatim
291*>
292*> \param[out] VS
293*> \verbatim
294*> VS is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLATMR, SLATMS, SLATME or DGEES 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 double_eig
383*
384* =====================================================================
385 SUBROUTINE ddrves( 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 DOUBLE PRECISION THRESH
396* ..
397* .. Array Arguments ..
398 LOGICAL BWORK( * ), DOTYPE( * )
399 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
400 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
401 $ result( 13 ), vs( ldvs, * ), wi( * ), wit( * ),
402 $ work( * ), wr( * ), wrt( * )
403* ..
404*
405* =====================================================================
406*
407* .. Parameters ..
408 DOUBLE PRECISION ZERO, ONE
409 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
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, nerrs,
419 $ nfail, nmax, nnwork, ntest, ntestf, ntestt,
420 $ rsub, sdim
421 DOUBLE PRECISION 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 DOUBLE PRECISION RES( 2 )
430* ..
431* .. Arrays in Common ..
432 LOGICAL SELVAL( 20 )
433 DOUBLE PRECISION 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 DSLECT
443 DOUBLE PRECISION DLAMCH
444 EXTERNAL dslect, dlamch
445* ..
446* .. External Subroutines ..
447 EXTERNAL dgees, dhst01, dlacpy, dlaset, dlasum, dlatme,
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 ) = 'Double 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( 'DDRVES', -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 = dlamch( 'Safe minimum' )
516 ovfl = one / unfl
517 ulp = dlamch( 'Precision' )
518 ulpinv = one / ulp
519 rtulp = sqrt( ulp )
520 rtulpi = one / rtulp
521*
522* Loop over sizes, types
523*
524 nerrs = 0
525*
526 DO 270 jsize = 1, nsizes
527 n = nn( jsize )
528 mtypes = maxtyp
529 IF( nsizes.EQ.1 .AND. ntypes.EQ.maxtyp+1 )
530 $ mtypes = mtypes + 1
531*
532 DO 260 jtype = 1, mtypes
533 IF( .NOT.dotype( jtype ) )
534 $ GO TO 260
535*
536* Save ISEED in case of an error.
537*
538 DO 20 j = 1, 4
539 ioldsd( j ) = iseed( j )
540 20 CONTINUE
541*
542* Compute "A"
543*
544* Control parameters:
545*
546* KMAGN KCONDS KMODE KTYPE
547* =1 O(1) 1 clustered 1 zero
548* =2 large large clustered 2 identity
549* =3 small exponential Jordan
550* =4 arithmetic diagonal, (w/ eigenvalues)
551* =5 random log symmetric, w/ eigenvalues
552* =6 random general, w/ eigenvalues
553* =7 random diagonal
554* =8 random symmetric
555* =9 random general
556* =10 random triangular
557*
558 IF( mtypes.GT.maxtyp )
559 $ GO TO 90
560*
561 itype = ktype( jtype )
562 imode = kmode( jtype )
563*
564* Compute norm
565*
566 GO TO ( 30, 40, 50 )kmagn( jtype )
567*
568 30 CONTINUE
569 anorm = one
570 GO TO 60
571*
572 40 CONTINUE
573 anorm = ovfl*ulp
574 GO TO 60
575*
576 50 CONTINUE
577 anorm = unfl*ulpinv
578 GO TO 60
579*
580 60 CONTINUE
581*
582 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
583 iinfo = 0
584 cond = ulpinv
585*
586* Special Matrices -- Identity & Jordan block
587*
588* Zero
589*
590 IF( itype.EQ.1 ) THEN
591 iinfo = 0
592*
593 ELSE IF( itype.EQ.2 ) THEN
594*
595* Identity
596*
597 DO 70 jcol = 1, n
598 a( jcol, jcol ) = anorm
599 70 CONTINUE
600*
601 ELSE IF( itype.EQ.3 ) THEN
602*
603* Jordan Block
604*
605 DO 80 jcol = 1, n
606 a( jcol, jcol ) = anorm
607 IF( jcol.GT.1 )
608 $ a( jcol, jcol-1 ) = one
609 80 CONTINUE
610*
611 ELSE IF( itype.EQ.4 ) THEN
612*
613* Diagonal Matrix, [Eigen]values Specified
614*
615 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
616 $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
617 $ iinfo )
618*
619 ELSE IF( itype.EQ.5 ) THEN
620*
621* Symmetric, eigenvalues specified
622*
623 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
624 $ anorm, n, n, 'N', a, lda, work( n+1 ),
625 $ iinfo )
626*
627 ELSE IF( itype.EQ.6 ) THEN
628*
629* General, eigenvalues specified
630*
631 IF( kconds( jtype ).EQ.1 ) THEN
632 conds = one
633 ELSE IF( kconds( jtype ).EQ.2 ) THEN
634 conds = rtulpi
635 ELSE
636 conds = zero
637 END IF
638*
639 adumma( 1 ) = ' '
640 CALL dlatme( n, 'S', iseed, work, imode, cond, one,
641 $ adumma, 'T', 'T', 'T', work( n+1 ), 4,
642 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
643 $ iinfo )
644*
645 ELSE IF( itype.EQ.7 ) THEN
646*
647* Diagonal, random eigenvalues
648*
649 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
650 $ 'T', 'N', work( n+1 ), 1, one,
651 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
652 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
653*
654 ELSE IF( itype.EQ.8 ) THEN
655*
656* Symmetric, random eigenvalues
657*
658 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
659 $ 'T', 'N', work( n+1 ), 1, one,
660 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
661 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
662*
663 ELSE IF( itype.EQ.9 ) THEN
664*
665* General, random eigenvalues
666*
667 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
668 $ 'T', 'N', work( n+1 ), 1, one,
669 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
670 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
671 IF( n.GE.4 ) THEN
672 CALL dlaset( 'Full', 2, n, zero, zero, a, lda )
673 CALL dlaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
674 $ lda )
675 CALL dlaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
676 $ lda )
677 CALL dlaset( 'Full', 1, n, zero, zero, a( n, 1 ),
678 $ lda )
679 END IF
680*
681 ELSE IF( itype.EQ.10 ) THEN
682*
683* Triangular, random eigenvalues
684*
685 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
686 $ 'T', 'N', work( n+1 ), 1, one,
687 $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
688 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
689*
690 ELSE
691*
692 iinfo = 1
693 END IF
694*
695 IF( iinfo.NE.0 ) THEN
696 WRITE( nounit, fmt = 9992 )'Generator', iinfo, n, jtype,
697 $ ioldsd
698 info = abs( iinfo )
699 RETURN
700 END IF
701*
702 90 CONTINUE
703*
704* Test for minimal and generous workspace
705*
706 DO 250 iwk = 1, 2
707 IF( iwk.EQ.1 ) THEN
708 nnwork = 3*n
709 ELSE
710 nnwork = 5*n + 2*n**2
711 END IF
712 nnwork = max( nnwork, 1 )
713*
714* Initialize RESULT
715*
716 DO 100 j = 1, 13
717 result( j ) = -one
718 100 CONTINUE
719*
720* Test with and without sorting of eigenvalues
721*
722 DO 210 isort = 0, 1
723 IF( isort.EQ.0 ) THEN
724 sort = 'N'
725 rsub = 0
726 ELSE
727 sort = 'S'
728 rsub = 6
729 END IF
730*
731* Compute Schur form and Schur vectors, and test them
732*
733 CALL dlacpy( 'F', n, n, a, lda, h, lda )
734 CALL dgees( 'V', sort, dslect, n, h, lda, sdim, wr,
735 $ wi, vs, ldvs, work, nnwork, bwork, iinfo )
736 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
737 result( 1+rsub ) = ulpinv
738 WRITE( nounit, fmt = 9992 )'DGEES1', iinfo, n,
739 $ jtype, ioldsd
740 info = abs( iinfo )
741 GO TO 220
742 END IF
743*
744* Do Test (1) or Test (7)
745*
746 result( 1+rsub ) = zero
747 DO 120 j = 1, n - 2
748 DO 110 i = j + 2, n
749 IF( h( i, j ).NE.zero )
750 $ result( 1+rsub ) = ulpinv
751 110 CONTINUE
752 120 CONTINUE
753 DO 130 i = 1, n - 2
754 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.
755 $ zero )result( 1+rsub ) = ulpinv
756 130 CONTINUE
757 DO 140 i = 1, n - 1
758 IF( h( i+1, i ).NE.zero ) THEN
759 IF( h( i, i ).NE.h( i+1, i+1 ) .OR.
760 $ h( i, i+1 ).EQ.zero .OR.
761 $ sign( one, h( i+1, i ) ).EQ.
762 $ sign( one, h( i, i+1 ) ) )result( 1+rsub )
763 $ = ulpinv
764 END IF
765 140 CONTINUE
766*
767* Do Tests (2) and (3) or Tests (8) and (9)
768*
769 lwork = max( 1, 2*n*n )
770 CALL dhst01( n, 1, n, a, lda, h, lda, vs, ldvs, work,
771 $ lwork, res )
772 result( 2+rsub ) = res( 1 )
773 result( 3+rsub ) = res( 2 )
774*
775* Do Test (4) or Test (10)
776*
777 result( 4+rsub ) = zero
778 DO 150 i = 1, n
779 IF( h( i, i ).NE.wr( i ) )
780 $ result( 4+rsub ) = ulpinv
781 150 CONTINUE
782 IF( n.GT.1 ) THEN
783 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
784 $ result( 4+rsub ) = ulpinv
785 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
786 $ result( 4+rsub ) = ulpinv
787 END IF
788 DO 160 i = 1, n - 1
789 IF( h( i+1, i ).NE.zero ) THEN
790 tmp = sqrt( abs( h( i+1, i ) ) )*
791 $ sqrt( abs( h( i, i+1 ) ) )
792 result( 4+rsub ) = max( result( 4+rsub ),
793 $ abs( wi( i )-tmp ) /
794 $ max( ulp*tmp, unfl ) )
795 result( 4+rsub ) = max( result( 4+rsub ),
796 $ abs( wi( i+1 )+tmp ) /
797 $ max( ulp*tmp, unfl ) )
798 ELSE IF( i.GT.1 ) THEN
799 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.
800 $ zero .AND. wi( i ).NE.zero )result( 4+rsub )
801 $ = ulpinv
802 END IF
803 160 CONTINUE
804*
805* Do Test (5) or Test (11)
806*
807 CALL dlacpy( 'F', n, n, a, lda, ht, lda )
808 CALL dgees( 'N', sort, dslect, n, ht, lda, sdim, wrt,
809 $ wit, vs, ldvs, work, nnwork, bwork,
810 $ iinfo )
811 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
812 result( 5+rsub ) = ulpinv
813 WRITE( nounit, fmt = 9992 )'DGEES2', iinfo, n,
814 $ jtype, ioldsd
815 info = abs( iinfo )
816 GO TO 220
817 END IF
818*
819 result( 5+rsub ) = zero
820 DO 180 j = 1, n
821 DO 170 i = 1, n
822 IF( h( i, j ).NE.ht( i, j ) )
823 $ result( 5+rsub ) = ulpinv
824 170 CONTINUE
825 180 CONTINUE
826*
827* Do Test (6) or Test (12)
828*
829 result( 6+rsub ) = zero
830 DO 190 i = 1, n
831 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
832 $ result( 6+rsub ) = ulpinv
833 190 CONTINUE
834*
835* Do Test (13)
836*
837 IF( isort.EQ.1 ) THEN
838 result( 13 ) = zero
839 knteig = 0
840 DO 200 i = 1, n
841 IF( dslect( wr( i ), wi( i ) ) .OR.
842 $ dslect( wr( i ), -wi( i ) ) )
843 $ knteig = knteig + 1
844 IF( i.LT.n ) THEN
845 IF( ( dslect( wr( i+1 ),
846 $ wi( i+1 ) ) .OR. dslect( wr( i+1 ),
847 $ -wi( i+1 ) ) ) .AND.
848 $ ( .NOT.( dslect( wr( i ),
849 $ wi( i ) ) .OR. dslect( wr( i ),
850 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )
851 $ result( 13 ) = ulpinv
852 END IF
853 200 CONTINUE
854 IF( sdim.NE.knteig ) THEN
855 result( 13 ) = ulpinv
856 END IF
857 END IF
858*
859 210 CONTINUE
860*
861* End of Loop -- Check for RESULT(j) > THRESH
862*
863 220 CONTINUE
864*
865 ntest = 0
866 nfail = 0
867 DO 230 j = 1, 13
868 IF( result( j ).GE.zero )
869 $ ntest = ntest + 1
870 IF( result( j ).GE.thresh )
871 $ nfail = nfail + 1
872 230 CONTINUE
873*
874 IF( nfail.GT.0 )
875 $ ntestf = ntestf + 1
876 IF( ntestf.EQ.1 ) THEN
877 WRITE( nounit, fmt = 9999 )path
878 WRITE( nounit, fmt = 9998 )
879 WRITE( nounit, fmt = 9997 )
880 WRITE( nounit, fmt = 9996 )
881 WRITE( nounit, fmt = 9995 )thresh
882 WRITE( nounit, fmt = 9994 )
883 ntestf = 2
884 END IF
885*
886 DO 240 j = 1, 13
887 IF( result( j ).GE.thresh ) THEN
888 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
889 $ j, result( j )
890 END IF
891 240 CONTINUE
892*
893 nerrs = nerrs + nfail
894 ntestt = ntestt + ntest
895*
896 250 CONTINUE
897 260 CONTINUE
898 270 CONTINUE
899*
900* Summary
901*
902 CALL dlasum( path, nounit, nerrs, ntestt )
903*
904 9999 FORMAT( / 1x, a3, ' -- Real Schur Form Decomposition Driver',
905 $ / ' Matrix types (see DDRVES for details): ' )
906*
907 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
908 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
909 $ / ' 2=Identity matrix. ', ' 6=Diagona',
910 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
911 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
912 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
913 $ 'mall, evenly spaced.' )
914 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
915 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
916 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
917 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
918 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
919 $ 'lex ', / ' 12=Well-cond., random complex ', 6x, ' ',
920 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
921 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
922 $ ' complx ' )
923 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
924 $ 'with small random entries.', / ' 20=Matrix with large ran',
925 $ 'dom entries. ', / )
926 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
927 $ / ' ( A denotes A on input and T denotes A on output)',
928 $ / / ' 1 = 0 if T in Schur form (no sort), ',
929 $ ' 1/ulp otherwise', /
930 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
931 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
932 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
933 $ ' 1/ulp otherwise', /
934 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
935 $ ' 1/ulp otherwise', /
936 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
937 $ ', 1/ulp otherwise' )
938 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
939 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
940 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
941 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
942 $ ' 1/ulp otherwise', /
943 $ ' 11 = 0 if T same no matter if VS computed (sort),',
944 $ ' 1/ulp otherwise', /
945 $ ' 12 = 0 if WR, WI same no matter if VS computed (sort),',
946 $ ' 1/ulp otherwise', /
947 $ ' 13 = 0 if sorting successful, 1/ulp otherwise', / )
948 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
949 $ ' type ', i2, ', test(', i2, ')=', g10.3 )
950 9992 FORMAT( ' DDRVES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
951 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
952*
953 RETURN
954*
955* End of DDRVES
956*
957 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ddrves(nsizes, nn, ntypes, dotype, iseed, thresh, nounit, a, lda, h, ht, wr, wi, wrt, wit, vs, ldvs, result, work, nwork, iwork, bwork, info)
DDRVES
Definition ddrves.f:388
subroutine dhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
DHST01
Definition dhst01.f:134
subroutine dlasum(type, iounit, ie, nrun)
DLASUM
Definition dlasum.f:43
subroutine dlatme(n, dist, iseed, d, mode, cond, dmax, ei, rsign, upper, sim, ds, modes, conds, kl, ku, anorm, a, lda, work, info)
DLATME
Definition dlatme.f:332
subroutine dlatmr(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)
DLATMR
Definition dlatmr.f:471
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
Definition dlatms.f:321
subroutine dgees(jobvs, sort, select, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info)
DGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
Definition dgees.f:216
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110