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