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