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