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, DLABAD, DLACPY, DLASET, DLASUM, 456 $ DLATME, DLATMR, DLATMS, XERBLA 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 CALL DLABAD( UNFL, OVFL ) 529 ULP = DLAMCH( 'Precision' ) 530 ULPINV = ONE / ULP 531 RTULP = SQRT( ULP ) 532 RTULPI = ONE / RTULP 533* 534* Loop over sizes, types 535* 536 NERRS = 0 537* 538 DO 270 JSIZE = 1, NSIZES 539 N = NN( JSIZE ) 540 IF( NSIZES.NE.1 ) THEN 541 MTYPES = MIN( MAXTYP, NTYPES ) 542 ELSE 543 MTYPES = MIN( MAXTYP+1, NTYPES ) 544 END IF 545* 546 DO 260 JTYPE = 1, MTYPES 547 IF( .NOT.DOTYPE( JTYPE ) ) 548 $ GO TO 260 549* 550* Save ISEED in case of an error. 551* 552 DO 20 J = 1, 4 553 IOLDSD( J ) = ISEED( J ) 554 20 CONTINUE 555* 556* Compute "A" 557* 558* Control parameters: 559* 560* KMAGN KCONDS KMODE KTYPE 561* =1 O(1) 1 clustered 1 zero 562* =2 large large clustered 2 identity 563* =3 small exponential Jordan 564* =4 arithmetic diagonal, (w/ eigenvalues) 565* =5 random log symmetric, w/ eigenvalues 566* =6 random general, w/ eigenvalues 567* =7 random diagonal 568* =8 random symmetric 569* =9 random general 570* =10 random triangular 571* 572 IF( MTYPES.GT.MAXTYP ) 573 $ GO TO 90 574* 575 ITYPE = KTYPE( JTYPE ) 576 IMODE = KMODE( JTYPE ) 577* 578* Compute norm 579* 580 GO TO ( 30, 40, 50 )KMAGN( JTYPE ) 581* 582 30 CONTINUE 583 ANORM = ONE 584 GO TO 60 585* 586 40 CONTINUE 587 ANORM = OVFL*ULP 588 GO TO 60 589* 590 50 CONTINUE 591 ANORM = UNFL*ULPINV 592 GO TO 60 593* 594 60 CONTINUE 595* 596 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 597 IINFO = 0 598 COND = ULPINV 599* 600* Special Matrices -- Identity & Jordan block 601* 602* Zero 603* 604 IF( ITYPE.EQ.1 ) THEN 605 IINFO = 0 606* 607 ELSE IF( ITYPE.EQ.2 ) THEN 608* 609* Identity 610* 611 DO 70 JCOL = 1, N 612 A( JCOL, JCOL ) = ANORM 613 70 CONTINUE 614* 615 ELSE IF( ITYPE.EQ.3 ) THEN 616* 617* Jordan Block 618* 619 DO 80 JCOL = 1, N 620 A( JCOL, JCOL ) = ANORM 621 IF( JCOL.GT.1 ) 622 $ A( JCOL, JCOL-1 ) = ONE 623 80 CONTINUE 624* 625 ELSE IF( ITYPE.EQ.4 ) THEN 626* 627* Diagonal Matrix, [Eigen]values Specified 628* 629 CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 630 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 631 $ IINFO ) 632* 633 ELSE IF( ITYPE.EQ.5 ) THEN 634* 635* Symmetric, eigenvalues specified 636* 637 CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 638 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 639 $ IINFO ) 640* 641 ELSE IF( ITYPE.EQ.6 ) THEN 642* 643* General, eigenvalues specified 644* 645 IF( KCONDS( JTYPE ).EQ.1 ) THEN 646 CONDS = ONE 647 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 648 CONDS = RTULPI 649 ELSE 650 CONDS = ZERO 651 END IF 652* 653 ADUMMA( 1 ) = ' ' 654 CALL DLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 655 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 656 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 657 $ IINFO ) 658* 659 ELSE IF( ITYPE.EQ.7 ) THEN 660* 661* Diagonal, random eigenvalues 662* 663 CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 664 $ 'T', 'N', WORK( N+1 ), 1, ONE, 665 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 666 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 667* 668 ELSE IF( ITYPE.EQ.8 ) THEN 669* 670* Symmetric, random eigenvalues 671* 672 CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 673 $ 'T', 'N', WORK( N+1 ), 1, ONE, 674 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 675 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 676* 677 ELSE IF( ITYPE.EQ.9 ) THEN 678* 679* General, random eigenvalues 680* 681 CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 682 $ 'T', 'N', WORK( N+1 ), 1, ONE, 683 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 684 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 685 IF( N.GE.4 ) THEN 686 CALL DLASET( 'Full', 2, N, ZERO, ZERO, A, LDA ) 687 CALL DLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ), 688 $ LDA ) 689 CALL DLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ), 690 $ LDA ) 691 CALL DLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ), 692 $ LDA ) 693 END IF 694* 695 ELSE IF( ITYPE.EQ.10 ) THEN 696* 697* Triangular, random eigenvalues 698* 699 CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 700 $ 'T', 'N', WORK( N+1 ), 1, ONE, 701 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 702 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 703* 704 ELSE 705* 706 IINFO = 1 707 END IF 708* 709 IF( IINFO.NE.0 ) THEN 710 WRITE( NOUNIT, FMT = 9993 )'Generator', IINFO, N, JTYPE, 711 $ IOLDSD 712 INFO = ABS( IINFO ) 713 RETURN 714 END IF 715* 716 90 CONTINUE 717* 718* Test for minimal and generous workspace 719* 720 DO 250 IWK = 1, 2 721 IF( IWK.EQ.1 ) THEN 722 NNWORK = 4*N 723 ELSE 724 NNWORK = 5*N + 2*N**2 725 END IF 726 NNWORK = MAX( NNWORK, 1 ) 727* 728* Initialize RESULT 729* 730 DO 100 J = 1, 7 731 RESULT( J ) = -ONE 732 100 CONTINUE 733* 734* Compute eigenvalues and eigenvectors, and test them 735* 736 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 737 CALL DGEEV( 'V', 'V', N, H, LDA, WR, WI, VL, LDVL, VR, 738 $ LDVR, WORK, NNWORK, IINFO ) 739 IF( IINFO.NE.0 ) THEN 740 RESULT( 1 ) = ULPINV 741 WRITE( NOUNIT, FMT = 9993 )'DGEEV1', IINFO, N, JTYPE, 742 $ IOLDSD 743 INFO = ABS( IINFO ) 744 GO TO 220 745 END IF 746* 747* Do Test (1) 748* 749 CALL DGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, WR, WI, 750 $ WORK, RES ) 751 RESULT( 1 ) = RES( 1 ) 752* 753* Do Test (2) 754* 755 CALL DGET22( 'T', 'N', 'T', N, A, LDA, VL, LDVL, WR, WI, 756 $ WORK, RES ) 757 RESULT( 2 ) = RES( 1 ) 758* 759* Do Test (3) 760* 761 DO 120 J = 1, N 762 TNRM = ONE 763 IF( WI( J ).EQ.ZERO ) THEN 764 TNRM = DNRM2( N, VR( 1, J ), 1 ) 765 ELSE IF( WI( J ).GT.ZERO ) THEN 766 TNRM = DLAPY2( DNRM2( N, VR( 1, J ), 1 ), 767 $ DNRM2( N, VR( 1, J+1 ), 1 ) ) 768 END IF 769 RESULT( 3 ) = MAX( RESULT( 3 ), 770 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 771 IF( WI( J ).GT.ZERO ) THEN 772 VMX = ZERO 773 VRMX = ZERO 774 DO 110 JJ = 1, N 775 VTST = DLAPY2( VR( JJ, J ), VR( JJ, J+1 ) ) 776 IF( VTST.GT.VMX ) 777 $ VMX = VTST 778 IF( VR( JJ, J+1 ).EQ.ZERO .AND. 779 $ ABS( VR( JJ, J ) ).GT.VRMX ) 780 $ VRMX = ABS( VR( JJ, J ) ) 781 110 CONTINUE 782 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 783 $ RESULT( 3 ) = ULPINV 784 END IF 785 120 CONTINUE 786* 787* Do Test (4) 788* 789 DO 140 J = 1, N 790 TNRM = ONE 791 IF( WI( J ).EQ.ZERO ) THEN 792 TNRM = DNRM2( N, VL( 1, J ), 1 ) 793 ELSE IF( WI( J ).GT.ZERO ) THEN 794 TNRM = DLAPY2( DNRM2( N, VL( 1, J ), 1 ), 795 $ DNRM2( N, VL( 1, J+1 ), 1 ) ) 796 END IF 797 RESULT( 4 ) = MAX( RESULT( 4 ), 798 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 799 IF( WI( J ).GT.ZERO ) THEN 800 VMX = ZERO 801 VRMX = ZERO 802 DO 130 JJ = 1, N 803 VTST = DLAPY2( VL( JJ, J ), VL( JJ, J+1 ) ) 804 IF( VTST.GT.VMX ) 805 $ VMX = VTST 806 IF( VL( JJ, J+1 ).EQ.ZERO .AND. 807 $ ABS( VL( JJ, J ) ).GT.VRMX ) 808 $ VRMX = ABS( VL( JJ, J ) ) 809 130 CONTINUE 810 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 811 $ RESULT( 4 ) = ULPINV 812 END IF 813 140 CONTINUE 814* 815* Compute eigenvalues only, and test them 816* 817 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 818 CALL DGEEV( 'N', 'N', N, H, LDA, WR1, WI1, DUM, 1, DUM, 819 $ 1, WORK, NNWORK, IINFO ) 820 IF( IINFO.NE.0 ) THEN 821 RESULT( 1 ) = ULPINV 822 WRITE( NOUNIT, FMT = 9993 )'DGEEV2', IINFO, N, JTYPE, 823 $ IOLDSD 824 INFO = ABS( IINFO ) 825 GO TO 220 826 END IF 827* 828* Do Test (5) 829* 830 DO 150 J = 1, N 831 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 832 $ RESULT( 5 ) = ULPINV 833 150 CONTINUE 834* 835* Compute eigenvalues and right eigenvectors, and test them 836* 837 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 838 CALL DGEEV( 'N', 'V', N, H, LDA, WR1, WI1, DUM, 1, LRE, 839 $ LDLRE, WORK, NNWORK, IINFO ) 840 IF( IINFO.NE.0 ) THEN 841 RESULT( 1 ) = ULPINV 842 WRITE( NOUNIT, FMT = 9993 )'DGEEV3', IINFO, N, JTYPE, 843 $ IOLDSD 844 INFO = ABS( IINFO ) 845 GO TO 220 846 END IF 847* 848* Do Test (5) again 849* 850 DO 160 J = 1, N 851 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 852 $ RESULT( 5 ) = ULPINV 853 160 CONTINUE 854* 855* Do Test (6) 856* 857 DO 180 J = 1, N 858 DO 170 JJ = 1, N 859 IF( VR( J, JJ ).NE.LRE( J, JJ ) ) 860 $ RESULT( 6 ) = ULPINV 861 170 CONTINUE 862 180 CONTINUE 863* 864* Compute eigenvalues and left eigenvectors, and test them 865* 866 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 867 CALL DGEEV( 'V', 'N', N, H, LDA, WR1, WI1, LRE, LDLRE, 868 $ DUM, 1, WORK, NNWORK, IINFO ) 869 IF( IINFO.NE.0 ) THEN 870 RESULT( 1 ) = ULPINV 871 WRITE( NOUNIT, FMT = 9993 )'DGEEV4', IINFO, N, JTYPE, 872 $ IOLDSD 873 INFO = ABS( IINFO ) 874 GO TO 220 875 END IF 876* 877* Do Test (5) again 878* 879 DO 190 J = 1, N 880 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 881 $ RESULT( 5 ) = ULPINV 882 190 CONTINUE 883* 884* Do Test (7) 885* 886 DO 210 J = 1, N 887 DO 200 JJ = 1, N 888 IF( VL( J, JJ ).NE.LRE( J, JJ ) ) 889 $ RESULT( 7 ) = ULPINV 890 200 CONTINUE 891 210 CONTINUE 892* 893* End of Loop -- Check for RESULT(j) > THRESH 894* 895 220 CONTINUE 896* 897 NTEST = 0 898 NFAIL = 0 899 DO 230 J = 1, 7 900 IF( RESULT( J ).GE.ZERO ) 901 $ NTEST = NTEST + 1 902 IF( RESULT( J ).GE.THRESH ) 903 $ NFAIL = NFAIL + 1 904 230 CONTINUE 905* 906 IF( NFAIL.GT.0 ) 907 $ NTESTF = NTESTF + 1 908 IF( NTESTF.EQ.1 ) THEN 909 WRITE( NOUNIT, FMT = 9999 )PATH 910 WRITE( NOUNIT, FMT = 9998 ) 911 WRITE( NOUNIT, FMT = 9997 ) 912 WRITE( NOUNIT, FMT = 9996 ) 913 WRITE( NOUNIT, FMT = 9995 )THRESH 914 NTESTF = 2 915 END IF 916* 917 DO 240 J = 1, 7 918 IF( RESULT( J ).GE.THRESH ) THEN 919 WRITE( NOUNIT, FMT = 9994 )N, IWK, IOLDSD, JTYPE, 920 $ J, RESULT( J ) 921 END IF 922 240 CONTINUE 923* 924 NERRS = NERRS + NFAIL 925 NTESTT = NTESTT + NTEST 926* 927 250 CONTINUE 928 260 CONTINUE 929 270 CONTINUE 930* 931* Summary 932* 933 CALL DLASUM( PATH, NOUNIT, NERRS, NTESTT ) 934* 935 9999 FORMAT( / 1X, A3, ' -- Real Eigenvalue-Eigenvector Decomposition', 936 $ ' Driver', / ' Matrix types (see DDRVEV for details): ' ) 937* 938 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ', 939 $ ' ', ' 5=Diagonal: geometr. spaced entries.', 940 $ / ' 2=Identity matrix. ', ' 6=Diagona', 941 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ', 942 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ', 943 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s', 944 $ 'mall, evenly spaced.' ) 945 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev', 946 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e', 947 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ', 948 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond', 949 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp', 950 $ 'lex ', / ' 12=Well-cond., random complex ', 6X, ' ', 951 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi', 952 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.', 953 $ ' complx ' ) 954 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ', 955 $ 'with small random entries.', / ' 20=Matrix with large ran', 956 $ 'dom entries. ', / ) 957 9995 FORMAT( ' Tests performed with test threshold =', F8.2, 958 $ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ', 959 $ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ', 960 $ / ' 3 = | |VR(i)| - 1 | / ulp ', 961 $ / ' 4 = | |VL(i)| - 1 | / ulp ', 962 $ / ' 5 = 0 if W same no matter if VR or VL computed,', 963 $ ' 1/ulp otherwise', / 964 $ ' 6 = 0 if VR same no matter if VL computed,', 965 $ ' 1/ulp otherwise', / 966 $ ' 7 = 0 if VL same no matter if VR computed,', 967 $ ' 1/ulp otherwise', / ) 968 9994 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ), 969 $ ' type ', I2, ', test(', I2, ')=', G10.3 ) 970 9993 FORMAT( ' DDRVEV: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 971 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 972* 973 RETURN 974* 975* End of DDRVEV 976* 977 END 978