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*> \date November 2011 400* 401*> \ingroup double_eig 402* 403* ===================================================================== 404 SUBROUTINE DDRVEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 405 $ NOUNIT, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, 406 $ VR, LDVR, LRE, LDLRE, RESULT, WORK, NWORK, 407 $ IWORK, INFO ) 408* 409* -- LAPACK test routine (version 3.4.0) -- 410* -- LAPACK is a software package provided by Univ. of Tennessee, -- 411* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 412* November 2011 413* 414* .. Scalar Arguments .. 415 INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NOUNIT, NSIZES, 416 $ NTYPES, NWORK 417 DOUBLE PRECISION THRESH 418* .. 419* .. Array Arguments .. 420 LOGICAL DOTYPE( * ) 421 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 422 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ), 423 $ RESULT( 7 ), VL( LDVL, * ), VR( LDVR, * ), 424 $ WI( * ), WI1( * ), WORK( * ), WR( * ), WR1( * ) 425* .. 426* 427* ===================================================================== 428* 429* .. Parameters .. 430 DOUBLE PRECISION ZERO, ONE 431 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 432 DOUBLE PRECISION TWO 433 PARAMETER ( TWO = 2.0D0 ) 434 INTEGER MAXTYP 435 PARAMETER ( MAXTYP = 21 ) 436* .. 437* .. Local Scalars .. 438 LOGICAL BADNN 439 CHARACTER*3 PATH 440 INTEGER IINFO, IMODE, ITYPE, IWK, J, JCOL, JJ, JSIZE, 441 $ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX, NNWORK, 442 $ NTEST, NTESTF, NTESTT 443 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TNRM, 444 $ ULP, ULPINV, UNFL, VMX, VRMX, VTST 445* .. 446* .. Local Arrays .. 447 CHARACTER ADUMMA( 1 ) 448 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), 449 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 450 $ KTYPE( MAXTYP ) 451 DOUBLE PRECISION DUM( 1 ), RES( 2 ) 452* .. 453* .. External Functions .. 454 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 455 EXTERNAL DLAMCH, DLAPY2, DNRM2 456* .. 457* .. External Subroutines .. 458 EXTERNAL DGEEV, DGET22, DLABAD, DLACPY, DLASET, DLASUM, 459 $ DLATME, DLATMR, DLATMS, XERBLA 460* .. 461* .. Intrinsic Functions .. 462 INTRINSIC ABS, MAX, MIN, SQRT 463* .. 464* .. Data statements .. 465 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 466 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 467 $ 3, 1, 2, 3 / 468 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 469 $ 1, 5, 5, 5, 4, 3, 1 / 470 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 471* .. 472* .. Executable Statements .. 473* 474 PATH( 1: 1 ) = 'Double precision' 475 PATH( 2: 3 ) = 'EV' 476* 477* Check for errors 478* 479 NTESTT = 0 480 NTESTF = 0 481 INFO = 0 482* 483* Important constants 484* 485 BADNN = .FALSE. 486 NMAX = 0 487 DO 10 J = 1, NSIZES 488 NMAX = MAX( NMAX, NN( J ) ) 489 IF( NN( J ).LT.0 ) 490 $ BADNN = .TRUE. 491 10 CONTINUE 492* 493* Check for errors 494* 495 IF( NSIZES.LT.0 ) THEN 496 INFO = -1 497 ELSE IF( BADNN ) THEN 498 INFO = -2 499 ELSE IF( NTYPES.LT.0 ) THEN 500 INFO = -3 501 ELSE IF( THRESH.LT.ZERO ) THEN 502 INFO = -6 503 ELSE IF( NOUNIT.LE.0 ) THEN 504 INFO = -7 505 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN 506 INFO = -9 507 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.NMAX ) THEN 508 INFO = -16 509 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.NMAX ) THEN 510 INFO = -18 511 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.NMAX ) THEN 512 INFO = -20 513 ELSE IF( 5*NMAX+2*NMAX**2.GT.NWORK ) THEN 514 INFO = -23 515 END IF 516* 517 IF( INFO.NE.0 ) THEN 518 CALL XERBLA( 'DDRVEV', -INFO ) 519 RETURN 520 END IF 521* 522* Quick return if nothing to do 523* 524 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 525 $ RETURN 526* 527* More Important constants 528* 529 UNFL = DLAMCH( 'Safe minimum' ) 530 OVFL = ONE / UNFL 531 CALL DLABAD( UNFL, OVFL ) 532 ULP = DLAMCH( 'Precision' ) 533 ULPINV = ONE / ULP 534 RTULP = SQRT( ULP ) 535 RTULPI = ONE / RTULP 536* 537* Loop over sizes, types 538* 539 NERRS = 0 540* 541 DO 270 JSIZE = 1, NSIZES 542 N = NN( JSIZE ) 543 IF( NSIZES.NE.1 ) THEN 544 MTYPES = MIN( MAXTYP, NTYPES ) 545 ELSE 546 MTYPES = MIN( MAXTYP+1, NTYPES ) 547 END IF 548* 549 DO 260 JTYPE = 1, MTYPES 550 IF( .NOT.DOTYPE( JTYPE ) ) 551 $ GO TO 260 552* 553* Save ISEED in case of an error. 554* 555 DO 20 J = 1, 4 556 IOLDSD( J ) = ISEED( J ) 557 20 CONTINUE 558* 559* Compute "A" 560* 561* Control parameters: 562* 563* KMAGN KCONDS KMODE KTYPE 564* =1 O(1) 1 clustered 1 zero 565* =2 large large clustered 2 identity 566* =3 small exponential Jordan 567* =4 arithmetic diagonal, (w/ eigenvalues) 568* =5 random log symmetric, w/ eigenvalues 569* =6 random general, w/ eigenvalues 570* =7 random diagonal 571* =8 random symmetric 572* =9 random general 573* =10 random triangular 574* 575 IF( MTYPES.GT.MAXTYP ) 576 $ GO TO 90 577* 578 ITYPE = KTYPE( JTYPE ) 579 IMODE = KMODE( JTYPE ) 580* 581* Compute norm 582* 583 GO TO ( 30, 40, 50 )KMAGN( JTYPE ) 584* 585 30 CONTINUE 586 ANORM = ONE 587 GO TO 60 588* 589 40 CONTINUE 590 ANORM = OVFL*ULP 591 GO TO 60 592* 593 50 CONTINUE 594 ANORM = UNFL*ULPINV 595 GO TO 60 596* 597 60 CONTINUE 598* 599 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 600 IINFO = 0 601 COND = ULPINV 602* 603* Special Matrices -- Identity & Jordan block 604* 605* Zero 606* 607 IF( ITYPE.EQ.1 ) THEN 608 IINFO = 0 609* 610 ELSE IF( ITYPE.EQ.2 ) THEN 611* 612* Identity 613* 614 DO 70 JCOL = 1, N 615 A( JCOL, JCOL ) = ANORM 616 70 CONTINUE 617* 618 ELSE IF( ITYPE.EQ.3 ) THEN 619* 620* Jordan Block 621* 622 DO 80 JCOL = 1, N 623 A( JCOL, JCOL ) = ANORM 624 IF( JCOL.GT.1 ) 625 $ A( JCOL, JCOL-1 ) = ONE 626 80 CONTINUE 627* 628 ELSE IF( ITYPE.EQ.4 ) THEN 629* 630* Diagonal Matrix, [Eigen]values Specified 631* 632 CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 633 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 634 $ IINFO ) 635* 636 ELSE IF( ITYPE.EQ.5 ) THEN 637* 638* Symmetric, eigenvalues specified 639* 640 CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 641 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 642 $ IINFO ) 643* 644 ELSE IF( ITYPE.EQ.6 ) THEN 645* 646* General, eigenvalues specified 647* 648 IF( KCONDS( JTYPE ).EQ.1 ) THEN 649 CONDS = ONE 650 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 651 CONDS = RTULPI 652 ELSE 653 CONDS = ZERO 654 END IF 655* 656 ADUMMA( 1 ) = ' ' 657 CALL DLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 658 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 659 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 660 $ IINFO ) 661* 662 ELSE IF( ITYPE.EQ.7 ) THEN 663* 664* Diagonal, random eigenvalues 665* 666 CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 667 $ 'T', 'N', WORK( N+1 ), 1, ONE, 668 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 669 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 670* 671 ELSE IF( ITYPE.EQ.8 ) THEN 672* 673* Symmetric, random eigenvalues 674* 675 CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 676 $ 'T', 'N', WORK( N+1 ), 1, ONE, 677 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 678 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 679* 680 ELSE IF( ITYPE.EQ.9 ) THEN 681* 682* General, random eigenvalues 683* 684 CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 685 $ 'T', 'N', WORK( N+1 ), 1, ONE, 686 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 687 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 688 IF( N.GE.4 ) THEN 689 CALL DLASET( 'Full', 2, N, ZERO, ZERO, A, LDA ) 690 CALL DLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ), 691 $ LDA ) 692 CALL DLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ), 693 $ LDA ) 694 CALL DLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ), 695 $ LDA ) 696 END IF 697* 698 ELSE IF( ITYPE.EQ.10 ) THEN 699* 700* Triangular, random eigenvalues 701* 702 CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 703 $ 'T', 'N', WORK( N+1 ), 1, ONE, 704 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 705 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 706* 707 ELSE 708* 709 IINFO = 1 710 END IF 711* 712 IF( IINFO.NE.0 ) THEN 713 WRITE( NOUNIT, FMT = 9993 )'Generator', IINFO, N, JTYPE, 714 $ IOLDSD 715 INFO = ABS( IINFO ) 716 RETURN 717 END IF 718* 719 90 CONTINUE 720* 721* Test for minimal and generous workspace 722* 723 DO 250 IWK = 1, 2 724 IF( IWK.EQ.1 ) THEN 725 NNWORK = 4*N 726 ELSE 727 NNWORK = 5*N + 2*N**2 728 END IF 729 NNWORK = MAX( NNWORK, 1 ) 730* 731* Initialize RESULT 732* 733 DO 100 J = 1, 7 734 RESULT( J ) = -ONE 735 100 CONTINUE 736* 737* Compute eigenvalues and eigenvectors, and test them 738* 739 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 740 CALL DGEEV( 'V', 'V', N, H, LDA, WR, WI, VL, LDVL, VR, 741 $ LDVR, WORK, NNWORK, IINFO ) 742 IF( IINFO.NE.0 ) THEN 743 RESULT( 1 ) = ULPINV 744 WRITE( NOUNIT, FMT = 9993 )'DGEEV1', IINFO, N, JTYPE, 745 $ IOLDSD 746 INFO = ABS( IINFO ) 747 GO TO 220 748 END IF 749* 750* Do Test (1) 751* 752 CALL DGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, WR, WI, 753 $ WORK, RES ) 754 RESULT( 1 ) = RES( 1 ) 755* 756* Do Test (2) 757* 758 CALL DGET22( 'T', 'N', 'T', N, A, LDA, VL, LDVL, WR, WI, 759 $ WORK, RES ) 760 RESULT( 2 ) = RES( 1 ) 761* 762* Do Test (3) 763* 764 DO 120 J = 1, N 765 TNRM = ONE 766 IF( WI( J ).EQ.ZERO ) THEN 767 TNRM = DNRM2( N, VR( 1, J ), 1 ) 768 ELSE IF( WI( J ).GT.ZERO ) THEN 769 TNRM = DLAPY2( DNRM2( N, VR( 1, J ), 1 ), 770 $ DNRM2( N, VR( 1, J+1 ), 1 ) ) 771 END IF 772 RESULT( 3 ) = MAX( RESULT( 3 ), 773 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 774 IF( WI( J ).GT.ZERO ) THEN 775 VMX = ZERO 776 VRMX = ZERO 777 DO 110 JJ = 1, N 778 VTST = DLAPY2( VR( JJ, J ), VR( JJ, J+1 ) ) 779 IF( VTST.GT.VMX ) 780 $ VMX = VTST 781 IF( VR( JJ, J+1 ).EQ.ZERO .AND. 782 $ ABS( VR( JJ, J ) ).GT.VRMX ) 783 $ VRMX = ABS( VR( JJ, J ) ) 784 110 CONTINUE 785 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 786 $ RESULT( 3 ) = ULPINV 787 END IF 788 120 CONTINUE 789* 790* Do Test (4) 791* 792 DO 140 J = 1, N 793 TNRM = ONE 794 IF( WI( J ).EQ.ZERO ) THEN 795 TNRM = DNRM2( N, VL( 1, J ), 1 ) 796 ELSE IF( WI( J ).GT.ZERO ) THEN 797 TNRM = DLAPY2( DNRM2( N, VL( 1, J ), 1 ), 798 $ DNRM2( N, VL( 1, J+1 ), 1 ) ) 799 END IF 800 RESULT( 4 ) = MAX( RESULT( 4 ), 801 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 802 IF( WI( J ).GT.ZERO ) THEN 803 VMX = ZERO 804 VRMX = ZERO 805 DO 130 JJ = 1, N 806 VTST = DLAPY2( VL( JJ, J ), VL( JJ, J+1 ) ) 807 IF( VTST.GT.VMX ) 808 $ VMX = VTST 809 IF( VL( JJ, J+1 ).EQ.ZERO .AND. 810 $ ABS( VL( JJ, J ) ).GT.VRMX ) 811 $ VRMX = ABS( VL( JJ, J ) ) 812 130 CONTINUE 813 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 814 $ RESULT( 4 ) = ULPINV 815 END IF 816 140 CONTINUE 817* 818* Compute eigenvalues only, and test them 819* 820 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 821 CALL DGEEV( 'N', 'N', N, H, LDA, WR1, WI1, DUM, 1, DUM, 822 $ 1, WORK, NNWORK, IINFO ) 823 IF( IINFO.NE.0 ) THEN 824 RESULT( 1 ) = ULPINV 825 WRITE( NOUNIT, FMT = 9993 )'DGEEV2', IINFO, N, JTYPE, 826 $ IOLDSD 827 INFO = ABS( IINFO ) 828 GO TO 220 829 END IF 830* 831* Do Test (5) 832* 833 DO 150 J = 1, N 834 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 835 $ RESULT( 5 ) = ULPINV 836 150 CONTINUE 837* 838* Compute eigenvalues and right eigenvectors, and test them 839* 840 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 841 CALL DGEEV( 'N', 'V', N, H, LDA, WR1, WI1, DUM, 1, LRE, 842 $ LDLRE, WORK, NNWORK, IINFO ) 843 IF( IINFO.NE.0 ) THEN 844 RESULT( 1 ) = ULPINV 845 WRITE( NOUNIT, FMT = 9993 )'DGEEV3', IINFO, N, JTYPE, 846 $ IOLDSD 847 INFO = ABS( IINFO ) 848 GO TO 220 849 END IF 850* 851* Do Test (5) again 852* 853 DO 160 J = 1, N 854 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 855 $ RESULT( 5 ) = ULPINV 856 160 CONTINUE 857* 858* Do Test (6) 859* 860 DO 180 J = 1, N 861 DO 170 JJ = 1, N 862 IF( VR( J, JJ ).NE.LRE( J, JJ ) ) 863 $ RESULT( 6 ) = ULPINV 864 170 CONTINUE 865 180 CONTINUE 866* 867* Compute eigenvalues and left eigenvectors, and test them 868* 869 CALL DLACPY( 'F', N, N, A, LDA, H, LDA ) 870 CALL DGEEV( 'V', 'N', N, H, LDA, WR1, WI1, LRE, LDLRE, 871 $ DUM, 1, WORK, NNWORK, IINFO ) 872 IF( IINFO.NE.0 ) THEN 873 RESULT( 1 ) = ULPINV 874 WRITE( NOUNIT, FMT = 9993 )'DGEEV4', IINFO, N, JTYPE, 875 $ IOLDSD 876 INFO = ABS( IINFO ) 877 GO TO 220 878 END IF 879* 880* Do Test (5) again 881* 882 DO 190 J = 1, N 883 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 884 $ RESULT( 5 ) = ULPINV 885 190 CONTINUE 886* 887* Do Test (7) 888* 889 DO 210 J = 1, N 890 DO 200 JJ = 1, N 891 IF( VL( J, JJ ).NE.LRE( J, JJ ) ) 892 $ RESULT( 7 ) = ULPINV 893 200 CONTINUE 894 210 CONTINUE 895* 896* End of Loop -- Check for RESULT(j) > THRESH 897* 898 220 CONTINUE 899* 900 NTEST = 0 901 NFAIL = 0 902 DO 230 J = 1, 7 903 IF( RESULT( J ).GE.ZERO ) 904 $ NTEST = NTEST + 1 905 IF( RESULT( J ).GE.THRESH ) 906 $ NFAIL = NFAIL + 1 907 230 CONTINUE 908* 909 IF( NFAIL.GT.0 ) 910 $ NTESTF = NTESTF + 1 911 IF( NTESTF.EQ.1 ) THEN 912 WRITE( NOUNIT, FMT = 9999 )PATH 913 WRITE( NOUNIT, FMT = 9998 ) 914 WRITE( NOUNIT, FMT = 9997 ) 915 WRITE( NOUNIT, FMT = 9996 ) 916 WRITE( NOUNIT, FMT = 9995 )THRESH 917 NTESTF = 2 918 END IF 919* 920 DO 240 J = 1, 7 921 IF( RESULT( J ).GE.THRESH ) THEN 922 WRITE( NOUNIT, FMT = 9994 )N, IWK, IOLDSD, JTYPE, 923 $ J, RESULT( J ) 924 END IF 925 240 CONTINUE 926* 927 NERRS = NERRS + NFAIL 928 NTESTT = NTESTT + NTEST 929* 930 250 CONTINUE 931 260 CONTINUE 932 270 CONTINUE 933* 934* Summary 935* 936 CALL DLASUM( PATH, NOUNIT, NERRS, NTESTT ) 937* 938 9999 FORMAT( / 1X, A3, ' -- Real Eigenvalue-Eigenvector Decomposition', 939 $ ' Driver', / ' Matrix types (see DDRVEV for details): ' ) 940* 941 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ', 942 $ ' ', ' 5=Diagonal: geometr. spaced entries.', 943 $ / ' 2=Identity matrix. ', ' 6=Diagona', 944 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ', 945 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ', 946 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s', 947 $ 'mall, evenly spaced.' ) 948 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev', 949 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e', 950 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ', 951 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond', 952 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp', 953 $ 'lex ', / ' 12=Well-cond., random complex ', 6X, ' ', 954 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi', 955 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.', 956 $ ' complx ' ) 957 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ', 958 $ 'with small random entries.', / ' 20=Matrix with large ran', 959 $ 'dom entries. ', / ) 960 9995 FORMAT( ' Tests performed with test threshold =', F8.2, 961 $ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ', 962 $ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ', 963 $ / ' 3 = | |VR(i)| - 1 | / ulp ', 964 $ / ' 4 = | |VL(i)| - 1 | / ulp ', 965 $ / ' 5 = 0 if W same no matter if VR or VL computed,', 966 $ ' 1/ulp otherwise', / 967 $ ' 6 = 0 if VR same no matter if VL computed,', 968 $ ' 1/ulp otherwise', / 969 $ ' 7 = 0 if VL same no matter if VR computed,', 970 $ ' 1/ulp otherwise', / ) 971 9994 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ), 972 $ ' type ', I2, ', test(', I2, ')=', G10.3 ) 973 9993 FORMAT( ' DDRVEV: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 974 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 975* 976 RETURN 977* 978* End of DDRVEV 979* 980 END 981