1*> \brief \b SCHKHS 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 SCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 12* NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, WR1, 13* WI1, WR2, WI2, WR3, WI3, EVECTL, EVECTR, EVECTY, 14* EVECTX, UU, TAU, WORK, NWORK, IWORK, SELECT, 15* RESULT, INFO ) 16* 17* .. Scalar Arguments .. 18* INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK 19* REAL THRESH 20* .. 21* .. Array Arguments .. 22* LOGICAL DOTYPE( * ), SELECT( * ) 23* INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 24* REAL A( LDA, * ), EVECTL( LDU, * ), 25* $ EVECTR( LDU, * ), EVECTX( LDU, * ), 26* $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), 27* $ T1( LDA, * ), T2( LDA, * ), TAU( * ), 28* $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), 29* $ WI1( * ), WI2( * ), WI3( * ), WORK( * ), 30* $ WR1( * ), WR2( * ), WR3( * ), Z( LDU, * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SCHKHS checks the nonsymmetric eigenvalue problem routines. 40*> 41*> SGEHRD factors A as U H U' , where ' means transpose, 42*> H is hessenberg, and U is an orthogonal matrix. 43*> 44*> SORGHR generates the orthogonal matrix U. 45*> 46*> SORMHR multiplies a matrix by the orthogonal matrix U. 47*> 48*> SHSEQR factors H as Z T Z' , where Z is orthogonal and 49*> T is "quasi-triangular", and the eigenvalue vector W. 50*> 51*> STREVC computes the left and right eigenvector matrices 52*> L and R for T. 53*> 54*> SHSEIN computes the left and right eigenvector matrices 55*> Y and X for H, using inverse iteration. 56*> 57*> When SCHKHS is called, a number of matrix "sizes" ("n's") and a 58*> number of matrix "types" are specified. For each size ("n") 59*> and each type of matrix, one matrix will be generated and used 60*> to test the nonsymmetric eigenroutines. For each matrix, 14 61*> tests will be performed: 62*> 63*> (1) | A - U H U**T | / ( |A| n ulp ) 64*> 65*> (2) | I - UU**T | / ( n ulp ) 66*> 67*> (3) | H - Z T Z**T | / ( |H| n ulp ) 68*> 69*> (4) | I - ZZ**T | / ( n ulp ) 70*> 71*> (5) | A - UZ H (UZ)**T | / ( |A| n ulp ) 72*> 73*> (6) | I - UZ (UZ)**T | / ( n ulp ) 74*> 75*> (7) | T(Z computed) - T(Z not computed) | / ( |T| ulp ) 76*> 77*> (8) | W(Z computed) - W(Z not computed) | / ( |W| ulp ) 78*> 79*> (9) | TR - RW | / ( |T| |R| ulp ) 80*> 81*> (10) | L**H T - W**H L | / ( |T| |L| ulp ) 82*> 83*> (11) | HX - XW | / ( |H| |X| ulp ) 84*> 85*> (12) | Y**H H - W**H Y | / ( |H| |Y| ulp ) 86*> 87*> (13) | AX - XW | / ( |A| |X| ulp ) 88*> 89*> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp ) 90*> 91*> The "sizes" are specified by an array NN(1:NSIZES); the value of 92*> each element NN(j) specifies one size. 93*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 94*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 95*> Currently, the list of possible types is: 96*> 97*> (1) The zero matrix. 98*> (2) The identity matrix. 99*> (3) A (transposed) Jordan block, with 1's on the diagonal. 100*> 101*> (4) A diagonal matrix with evenly spaced entries 102*> 1, ..., ULP and random signs. 103*> (ULP = (first number larger than 1) - 1 ) 104*> (5) A diagonal matrix with geometrically spaced entries 105*> 1, ..., ULP and random signs. 106*> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 107*> and random signs. 108*> 109*> (7) Same as (4), but multiplied by SQRT( overflow threshold ) 110*> (8) Same as (4), but multiplied by SQRT( underflow threshold ) 111*> 112*> (9) A matrix of the form U' T U, where U is orthogonal and 113*> T has evenly spaced entries 1, ..., ULP with random signs 114*> on the diagonal and random O(1) entries in the upper 115*> triangle. 116*> 117*> (10) A matrix of the form U' T U, where U is orthogonal and 118*> T has geometrically spaced entries 1, ..., ULP with random 119*> signs on the diagonal and random O(1) entries in the upper 120*> triangle. 121*> 122*> (11) A matrix of the form U' T U, where U is orthogonal and 123*> T has "clustered" entries 1, ULP,..., ULP with random 124*> signs on the diagonal and random O(1) entries in the upper 125*> triangle. 126*> 127*> (12) A matrix of the form U' T U, where U is orthogonal and 128*> T has real or complex conjugate paired eigenvalues randomly 129*> chosen from ( ULP, 1 ) and random O(1) entries in the upper 130*> triangle. 131*> 132*> (13) A matrix of the form X' T X, where X has condition 133*> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP 134*> with random signs on the diagonal and random O(1) entries 135*> in the upper triangle. 136*> 137*> (14) A matrix of the form X' T X, where X has condition 138*> SQRT( ULP ) and T has geometrically spaced entries 139*> 1, ..., ULP with random signs on the diagonal and random 140*> O(1) entries in the upper triangle. 141*> 142*> (15) A matrix of the form X' T X, where X has condition 143*> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP 144*> with random signs on the diagonal and random O(1) entries 145*> in the upper triangle. 146*> 147*> (16) A matrix of the form X' T X, where X has condition 148*> SQRT( ULP ) and T has real or complex conjugate paired 149*> eigenvalues randomly chosen from ( ULP, 1 ) and random 150*> O(1) entries in the upper triangle. 151*> 152*> (17) Same as (16), but multiplied by SQRT( overflow threshold ) 153*> (18) Same as (16), but multiplied by SQRT( underflow threshold ) 154*> 155*> (19) Nonsymmetric matrix with random entries chosen from (-1,1). 156*> (20) Same as (19), but multiplied by SQRT( overflow threshold ) 157*> (21) Same as (19), but multiplied by SQRT( underflow threshold ) 158*> \endverbatim 159* 160* Arguments: 161* ========== 162* 163*> \verbatim 164*> NSIZES - INTEGER 165*> The number of sizes of matrices to use. If it is zero, 166*> SCHKHS does nothing. It must be at least zero. 167*> Not modified. 168*> 169*> NN - INTEGER array, dimension (NSIZES) 170*> An array containing the sizes to be used for the matrices. 171*> Zero values will be skipped. The values must be at least 172*> zero. 173*> Not modified. 174*> 175*> NTYPES - INTEGER 176*> The number of elements in DOTYPE. If it is zero, SCHKHS 177*> does nothing. It must be at least zero. If it is MAXTYP+1 178*> and NSIZES is 1, then an additional type, MAXTYP+1 is 179*> defined, which is to use whatever matrix is in A. This 180*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 181*> DOTYPE(MAXTYP+1) is .TRUE. . 182*> Not modified. 183*> 184*> DOTYPE - LOGICAL array, dimension (NTYPES) 185*> If DOTYPE(j) is .TRUE., then for each size in NN a 186*> matrix of that size and of type j will be generated. 187*> If NTYPES is smaller than the maximum number of types 188*> defined (PARAMETER MAXTYP), then types NTYPES+1 through 189*> MAXTYP will not be generated. If NTYPES is larger 190*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 191*> will be ignored. 192*> Not modified. 193*> 194*> ISEED - INTEGER array, dimension (4) 195*> On entry ISEED specifies the seed of the random number 196*> generator. The array elements should be between 0 and 4095; 197*> if not they will be reduced mod 4096. Also, ISEED(4) must 198*> be odd. The random number generator uses a linear 199*> congruential sequence limited to small integers, and so 200*> should produce machine independent random numbers. The 201*> values of ISEED are changed on exit, and can be used in the 202*> next call to SCHKHS to continue the same random number 203*> sequence. 204*> Modified. 205*> 206*> THRESH - REAL 207*> A test will count as "failed" if the "error", computed as 208*> described above, exceeds THRESH. Note that the error 209*> is scaled to be O(1), so THRESH should be a reasonably 210*> small multiple of 1, e.g., 10 or 100. In particular, 211*> it should not depend on the precision (single vs. double) 212*> or the size of the matrix. It must be at least zero. 213*> Not modified. 214*> 215*> NOUNIT - INTEGER 216*> The FORTRAN unit number for printing out error messages 217*> (e.g., if a routine returns IINFO not equal to 0.) 218*> Not modified. 219*> 220*> A - REAL array, dimension (LDA,max(NN)) 221*> Used to hold the matrix whose eigenvalues are to be 222*> computed. On exit, A contains the last matrix actually 223*> used. 224*> Modified. 225*> 226*> LDA - INTEGER 227*> The leading dimension of A, H, T1 and T2. It must be at 228*> least 1 and at least max( NN ). 229*> Not modified. 230*> 231*> H - REAL array, dimension (LDA,max(NN)) 232*> The upper hessenberg matrix computed by SGEHRD. On exit, 233*> H contains the Hessenberg form of the matrix in A. 234*> Modified. 235*> 236*> T1 - REAL array, dimension (LDA,max(NN)) 237*> The Schur (="quasi-triangular") matrix computed by SHSEQR 238*> if Z is computed. On exit, T1 contains the Schur form of 239*> the matrix in A. 240*> Modified. 241*> 242*> T2 - REAL array, dimension (LDA,max(NN)) 243*> The Schur matrix computed by SHSEQR when Z is not computed. 244*> This should be identical to T1. 245*> Modified. 246*> 247*> LDU - INTEGER 248*> The leading dimension of U, Z, UZ and UU. It must be at 249*> least 1 and at least max( NN ). 250*> Not modified. 251*> 252*> U - REAL array, dimension (LDU,max(NN)) 253*> The orthogonal matrix computed by SGEHRD. 254*> Modified. 255*> 256*> Z - REAL array, dimension (LDU,max(NN)) 257*> The orthogonal matrix computed by SHSEQR. 258*> Modified. 259*> 260*> UZ - REAL array, dimension (LDU,max(NN)) 261*> The product of U times Z. 262*> Modified. 263*> 264*> WR1 - REAL array, dimension (max(NN)) 265*> WI1 - REAL array, dimension (max(NN)) 266*> The real and imaginary parts of the eigenvalues of A, 267*> as computed when Z is computed. 268*> On exit, WR1 + WI1*i are the eigenvalues of the matrix in A. 269*> Modified. 270*> 271*> WR2 - REAL array, dimension (max(NN)) 272*> WI2 - REAL array, dimension (max(NN)) 273*> The real and imaginary parts of the eigenvalues of A, 274*> as computed when T is computed but not Z. 275*> On exit, WR2 + WI2*i are the eigenvalues of the matrix in A. 276*> Modified. 277*> 278*> WR3 - REAL array, dimension (max(NN)) 279*> WI3 - REAL array, dimension (max(NN)) 280*> Like WR1, WI1, these arrays contain the eigenvalues of A, 281*> but those computed when SHSEQR only computes the 282*> eigenvalues, i.e., not the Schur vectors and no more of the 283*> Schur form than is necessary for computing the 284*> eigenvalues. 285*> Modified. 286*> 287*> EVECTL - REAL array, dimension (LDU,max(NN)) 288*> The (upper triangular) left eigenvector matrix for the 289*> matrix in T1. For complex conjugate pairs, the real part 290*> is stored in one row and the imaginary part in the next. 291*> Modified. 292*> 293*> EVECTR - REAL array, dimension (LDU,max(NN)) 294*> The (upper triangular) right eigenvector matrix for the 295*> matrix in T1. For complex conjugate pairs, the real part 296*> is stored in one column and the imaginary part in the next. 297*> Modified. 298*> 299*> EVECTY - REAL array, dimension (LDU,max(NN)) 300*> The left eigenvector matrix for the 301*> matrix in H. For complex conjugate pairs, the real part 302*> is stored in one row and the imaginary part in the next. 303*> Modified. 304*> 305*> EVECTX - REAL array, dimension (LDU,max(NN)) 306*> The right eigenvector matrix for the 307*> matrix in H. For complex conjugate pairs, the real part 308*> is stored in one column and the imaginary part in the next. 309*> Modified. 310*> 311*> UU - REAL array, dimension (LDU,max(NN)) 312*> Details of the orthogonal matrix computed by SGEHRD. 313*> Modified. 314*> 315*> TAU - REAL array, dimension(max(NN)) 316*> Further details of the orthogonal matrix computed by SGEHRD. 317*> Modified. 318*> 319*> WORK - REAL array, dimension (NWORK) 320*> Workspace. 321*> Modified. 322*> 323*> NWORK - INTEGER 324*> The number of entries in WORK. NWORK >= 4*NN(j)*NN(j) + 2. 325*> 326*> IWORK - INTEGER array, dimension (max(NN)) 327*> Workspace. 328*> Modified. 329*> 330*> SELECT - LOGICAL array, dimension (max(NN)) 331*> Workspace. 332*> Modified. 333*> 334*> RESULT - REAL array, dimension (14) 335*> The values computed by the fourteen tests described above. 336*> The values are currently limited to 1/ulp, to avoid 337*> overflow. 338*> Modified. 339*> 340*> INFO - INTEGER 341*> If 0, then everything ran OK. 342*> -1: NSIZES < 0 343*> -2: Some NN(j) < 0 344*> -3: NTYPES < 0 345*> -6: THRESH < 0 346*> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). 347*> -14: LDU < 1 or LDU < NMAX. 348*> -28: NWORK too small. 349*> If SLATMR, SLATMS, or SLATME returns an error code, the 350*> absolute value of it is returned. 351*> If 1, then SHSEQR could not find all the shifts. 352*> If 2, then the EISPACK code (for small blocks) failed. 353*> If >2, then 30*N iterations were not enough to find an 354*> eigenvalue or to decompose the problem. 355*> Modified. 356*> 357*>----------------------------------------------------------------------- 358*> 359*> Some Local Variables and Parameters: 360*> ---- ----- --------- --- ---------- 361*> 362*> ZERO, ONE Real 0 and 1. 363*> MAXTYP The number of types defined. 364*> MTEST The number of tests defined: care must be taken 365*> that (1) the size of RESULT, (2) the number of 366*> tests actually performed, and (3) MTEST agree. 367*> NTEST The number of tests performed on this matrix 368*> so far. This should be less than MTEST, and 369*> equal to it by the last test. It will be less 370*> if any of the routines being tested indicates 371*> that it could not compute the matrices that 372*> would be tested. 373*> NMAX Largest value in NN. 374*> NMATS The number of matrices generated so far. 375*> NERRS The number of tests which have exceeded THRESH 376*> so far (computed by SLAFTS). 377*> COND, CONDS, 378*> IMODE Values to be passed to the matrix generators. 379*> ANORM Norm of A; passed to matrix generators. 380*> 381*> OVFL, UNFL Overflow and underflow thresholds. 382*> ULP, ULPINV Finest relative precision and its inverse. 383*> RTOVFL, RTUNFL, 384*> RTULP, RTULPI Square roots of the previous 4 values. 385*> 386*> The following four arrays decode JTYPE: 387*> KTYPE(j) The general type (1-10) for type "j". 388*> KMODE(j) The MODE value to be passed to the matrix 389*> generator for type "j". 390*> KMAGN(j) The order of magnitude ( O(1), 391*> O(overflow^(1/2) ), O(underflow^(1/2) ) 392*> KCONDS(j) Selects whether CONDS is to be 1 or 393*> 1/sqrt(ulp). (0 means irrelevant.) 394*> \endverbatim 395* 396* Authors: 397* ======== 398* 399*> \author Univ. of Tennessee 400*> \author Univ. of California Berkeley 401*> \author Univ. of Colorado Denver 402*> \author NAG Ltd. 403* 404*> \ingroup single_eig 405* 406* ===================================================================== 407 SUBROUTINE SCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 408 $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, WR1, 409 $ WI1, WR2, WI2, WR3, WI3, EVECTL, EVECTR, 410 $ EVECTY, EVECTX, UU, TAU, WORK, NWORK, IWORK, 411 $ SELECT, RESULT, INFO ) 412* 413* -- LAPACK test routine -- 414* -- LAPACK is a software package provided by Univ. of Tennessee, -- 415* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 416* 417* .. Scalar Arguments .. 418 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK 419 REAL THRESH 420* .. 421* .. Array Arguments .. 422 LOGICAL DOTYPE( * ), SELECT( * ) 423 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 424 REAL A( LDA, * ), EVECTL( LDU, * ), 425 $ EVECTR( LDU, * ), EVECTX( LDU, * ), 426 $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), 427 $ T1( LDA, * ), T2( LDA, * ), TAU( * ), 428 $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), 429 $ WI1( * ), WI2( * ), WI3( * ), WORK( * ), 430 $ WR1( * ), WR2( * ), WR3( * ), Z( LDU, * ) 431* .. 432* 433* ===================================================================== 434* 435* .. Parameters .. 436 REAL ZERO, ONE 437 PARAMETER ( ZERO = 0.0, ONE = 1.0 ) 438 INTEGER MAXTYP 439 PARAMETER ( MAXTYP = 21 ) 440* .. 441* .. Local Scalars .. 442 LOGICAL BADNN, MATCH 443 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL, 444 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS, 445 $ NMATS, NMAX, NSELC, NSELR, NTEST, NTESTT 446 REAL ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP, 447 $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL 448* .. 449* .. Local Arrays .. 450 CHARACTER ADUMMA( 1 ) 451 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), 452 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 453 $ KTYPE( MAXTYP ) 454 REAL DUMMA( 6 ) 455* .. 456* .. External Functions .. 457 REAL SLAMCH 458 EXTERNAL SLAMCH 459* .. 460* .. External Subroutines .. 461 EXTERNAL SCOPY, SGEHRD, SGEMM, SGET10, SGET22, SHSEIN, 462 $ SHSEQR, SHST01, SLABAD, SLACPY, SLAFTS, SLASET, 463 $ SLASUM, SLATME, SLATMR, SLATMS, SORGHR, SORMHR, 464 $ STREVC, XERBLA 465* .. 466* .. Intrinsic Functions .. 467 INTRINSIC ABS, MAX, MIN, REAL, SQRT 468* .. 469* .. Data statements .. 470 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 471 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 472 $ 3, 1, 2, 3 / 473 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 474 $ 1, 5, 5, 5, 4, 3, 1 / 475 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 476* .. 477* .. Executable Statements .. 478* 479* Check for errors 480* 481 NTESTT = 0 482 INFO = 0 483* 484 BADNN = .FALSE. 485 NMAX = 0 486 DO 10 J = 1, NSIZES 487 NMAX = MAX( NMAX, NN( J ) ) 488 IF( NN( J ).LT.0 ) 489 $ BADNN = .TRUE. 490 10 CONTINUE 491* 492* Check for errors 493* 494 IF( NSIZES.LT.0 ) THEN 495 INFO = -1 496 ELSE IF( BADNN ) THEN 497 INFO = -2 498 ELSE IF( NTYPES.LT.0 ) THEN 499 INFO = -3 500 ELSE IF( THRESH.LT.ZERO ) THEN 501 INFO = -6 502 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 503 INFO = -9 504 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN 505 INFO = -14 506 ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN 507 INFO = -28 508 END IF 509* 510 IF( INFO.NE.0 ) THEN 511 CALL XERBLA( 'SCHKHS', -INFO ) 512 RETURN 513 END IF 514* 515* Quick return if possible 516* 517 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 518 $ RETURN 519* 520* More important constants 521* 522 UNFL = SLAMCH( 'Safe minimum' ) 523 OVFL = SLAMCH( 'Overflow' ) 524 CALL SLABAD( UNFL, OVFL ) 525 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 526 ULPINV = ONE / ULP 527 RTUNFL = SQRT( UNFL ) 528 RTOVFL = SQRT( OVFL ) 529 RTULP = SQRT( ULP ) 530 RTULPI = ONE / RTULP 531* 532* Loop over sizes, types 533* 534 NERRS = 0 535 NMATS = 0 536* 537 DO 270 JSIZE = 1, NSIZES 538 N = NN( JSIZE ) 539 IF( N.EQ.0 ) 540 $ GO TO 270 541 N1 = MAX( 1, N ) 542 ANINV = ONE / REAL( N1 ) 543* 544 IF( NSIZES.NE.1 ) THEN 545 MTYPES = MIN( MAXTYP, NTYPES ) 546 ELSE 547 MTYPES = MIN( MAXTYP+1, NTYPES ) 548 END IF 549* 550 DO 260 JTYPE = 1, MTYPES 551 IF( .NOT.DOTYPE( JTYPE ) ) 552 $ GO TO 260 553 NMATS = NMATS + 1 554 NTEST = 0 555* 556* Save ISEED in case of an error. 557* 558 DO 20 J = 1, 4 559 IOLDSD( J ) = ISEED( J ) 560 20 CONTINUE 561* 562* Initialize RESULT 563* 564 DO 30 J = 1, 14 565 RESULT( J ) = ZERO 566 30 CONTINUE 567* 568* Compute "A" 569* 570* Control parameters: 571* 572* KMAGN KCONDS KMODE KTYPE 573* =1 O(1) 1 clustered 1 zero 574* =2 large large clustered 2 identity 575* =3 small exponential Jordan 576* =4 arithmetic diagonal, (w/ eigenvalues) 577* =5 random log symmetric, w/ eigenvalues 578* =6 random general, w/ eigenvalues 579* =7 random diagonal 580* =8 random symmetric 581* =9 random general 582* =10 random triangular 583* 584 IF( MTYPES.GT.MAXTYP ) 585 $ GO TO 100 586* 587 ITYPE = KTYPE( JTYPE ) 588 IMODE = KMODE( JTYPE ) 589* 590* Compute norm 591* 592 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 593* 594 40 CONTINUE 595 ANORM = ONE 596 GO TO 70 597* 598 50 CONTINUE 599 ANORM = ( RTOVFL*ULP )*ANINV 600 GO TO 70 601* 602 60 CONTINUE 603 ANORM = RTUNFL*N*ULPINV 604 GO TO 70 605* 606 70 CONTINUE 607* 608 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 609 IINFO = 0 610 COND = ULPINV 611* 612* Special Matrices 613* 614 IF( ITYPE.EQ.1 ) THEN 615* 616* Zero 617* 618 IINFO = 0 619* 620 ELSE IF( ITYPE.EQ.2 ) THEN 621* 622* Identity 623* 624 DO 80 JCOL = 1, N 625 A( JCOL, JCOL ) = ANORM 626 80 CONTINUE 627* 628 ELSE IF( ITYPE.EQ.3 ) THEN 629* 630* Jordan Block 631* 632 DO 90 JCOL = 1, N 633 A( JCOL, JCOL ) = ANORM 634 IF( JCOL.GT.1 ) 635 $ A( JCOL, JCOL-1 ) = ONE 636 90 CONTINUE 637* 638 ELSE IF( ITYPE.EQ.4 ) THEN 639* 640* Diagonal Matrix, [Eigen]values Specified 641* 642 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 643 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 644 $ IINFO ) 645* 646 ELSE IF( ITYPE.EQ.5 ) THEN 647* 648* Symmetric, eigenvalues specified 649* 650 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 651 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 652 $ IINFO ) 653* 654 ELSE IF( ITYPE.EQ.6 ) THEN 655* 656* General, eigenvalues specified 657* 658 IF( KCONDS( JTYPE ).EQ.1 ) THEN 659 CONDS = ONE 660 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 661 CONDS = RTULPI 662 ELSE 663 CONDS = ZERO 664 END IF 665* 666 ADUMMA( 1 ) = ' ' 667 CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 668 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 669 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 670 $ IINFO ) 671* 672 ELSE IF( ITYPE.EQ.7 ) THEN 673* 674* Diagonal, random eigenvalues 675* 676 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 677 $ 'T', 'N', WORK( N+1 ), 1, ONE, 678 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 679 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 680* 681 ELSE IF( ITYPE.EQ.8 ) THEN 682* 683* Symmetric, random eigenvalues 684* 685 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 686 $ 'T', 'N', WORK( N+1 ), 1, ONE, 687 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 688 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 689* 690 ELSE IF( ITYPE.EQ.9 ) THEN 691* 692* General, random eigenvalues 693* 694 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 695 $ 'T', 'N', WORK( N+1 ), 1, ONE, 696 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 697 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 698* 699 ELSE IF( ITYPE.EQ.10 ) THEN 700* 701* Triangular, random eigenvalues 702* 703 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 704 $ 'T', 'N', WORK( N+1 ), 1, ONE, 705 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 706 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 707* 708 ELSE 709* 710 IINFO = 1 711 END IF 712* 713 IF( IINFO.NE.0 ) THEN 714 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 715 $ IOLDSD 716 INFO = ABS( IINFO ) 717 RETURN 718 END IF 719* 720 100 CONTINUE 721* 722* Call SGEHRD to compute H and U, do tests. 723* 724 CALL SLACPY( ' ', N, N, A, LDA, H, LDA ) 725* 726 NTEST = 1 727* 728 ILO = 1 729 IHI = N 730* 731 CALL SGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ), 732 $ NWORK-N, IINFO ) 733* 734 IF( IINFO.NE.0 ) THEN 735 RESULT( 1 ) = ULPINV 736 WRITE( NOUNIT, FMT = 9999 )'SGEHRD', IINFO, N, JTYPE, 737 $ IOLDSD 738 INFO = ABS( IINFO ) 739 GO TO 250 740 END IF 741* 742 DO 120 J = 1, N - 1 743 UU( J+1, J ) = ZERO 744 DO 110 I = J + 2, N 745 U( I, J ) = H( I, J ) 746 UU( I, J ) = H( I, J ) 747 H( I, J ) = ZERO 748 110 CONTINUE 749 120 CONTINUE 750 CALL SCOPY( N-1, WORK, 1, TAU, 1 ) 751 CALL SORGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ), 752 $ NWORK-N, IINFO ) 753 NTEST = 2 754* 755 CALL SHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK, 756 $ NWORK, RESULT( 1 ) ) 757* 758* Call SHSEQR to compute T1, T2 and Z, do tests. 759* 760* Eigenvalues only (WR3,WI3) 761* 762 CALL SLACPY( ' ', N, N, H, LDA, T2, LDA ) 763 NTEST = 3 764 RESULT( 3 ) = ULPINV 765* 766 CALL SHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, WR3, WI3, UZ, 767 $ LDU, WORK, NWORK, IINFO ) 768 IF( IINFO.NE.0 ) THEN 769 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(E)', IINFO, N, JTYPE, 770 $ IOLDSD 771 IF( IINFO.LE.N+2 ) THEN 772 INFO = ABS( IINFO ) 773 GO TO 250 774 END IF 775 END IF 776* 777* Eigenvalues (WR2,WI2) and Full Schur Form (T2) 778* 779 CALL SLACPY( ' ', N, N, H, LDA, T2, LDA ) 780* 781 CALL SHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, WR2, WI2, UZ, 782 $ LDU, WORK, NWORK, IINFO ) 783 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 784 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(S)', IINFO, N, JTYPE, 785 $ IOLDSD 786 INFO = ABS( IINFO ) 787 GO TO 250 788 END IF 789* 790* Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors 791* (UZ) 792* 793 CALL SLACPY( ' ', N, N, H, LDA, T1, LDA ) 794 CALL SLACPY( ' ', N, N, U, LDU, UZ, LDU ) 795* 796 CALL SHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, WR1, WI1, UZ, 797 $ LDU, WORK, NWORK, IINFO ) 798 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 799 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(V)', IINFO, N, JTYPE, 800 $ IOLDSD 801 INFO = ABS( IINFO ) 802 GO TO 250 803 END IF 804* 805* Compute Z = U' UZ 806* 807 CALL SGEMM( 'T', 'N', N, N, N, ONE, U, LDU, UZ, LDU, ZERO, 808 $ Z, LDU ) 809 NTEST = 8 810* 811* Do Tests 3: | H - Z T Z' | / ( |H| n ulp ) 812* and 4: | I - Z Z' | / ( n ulp ) 813* 814 CALL SHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK, 815 $ NWORK, RESULT( 3 ) ) 816* 817* Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp ) 818* and 6: | I - UZ (UZ)' | / ( n ulp ) 819* 820 CALL SHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK, 821 $ NWORK, RESULT( 5 ) ) 822* 823* Do Test 7: | T2 - T1 | / ( |T| n ulp ) 824* 825 CALL SGET10( N, N, T2, LDA, T1, LDA, WORK, RESULT( 7 ) ) 826* 827* Do Test 8: | W2 - W1 | / ( max(|W1|,|W2|) ulp ) 828* 829 TEMP1 = ZERO 830 TEMP2 = ZERO 831 DO 130 J = 1, N 832 TEMP1 = MAX( TEMP1, ABS( WR1( J ) )+ABS( WI1( J ) ), 833 $ ABS( WR2( J ) )+ABS( WI2( J ) ) ) 834 TEMP2 = MAX( TEMP2, ABS( WR1( J )-WR2( J ) )+ 835 $ ABS( WI1( J )-WI2( J ) ) ) 836 130 CONTINUE 837* 838 RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 839* 840* Compute the Left and Right Eigenvectors of T 841* 842* Compute the Right eigenvector Matrix: 843* 844 NTEST = 9 845 RESULT( 9 ) = ULPINV 846* 847* Select last max(N/4,1) real, max(N/4,1) complex eigenvectors 848* 849 NSELC = 0 850 NSELR = 0 851 J = N 852 140 CONTINUE 853 IF( WI1( J ).EQ.ZERO ) THEN 854 IF( NSELR.LT.MAX( N / 4, 1 ) ) THEN 855 NSELR = NSELR + 1 856 SELECT( J ) = .TRUE. 857 ELSE 858 SELECT( J ) = .FALSE. 859 END IF 860 J = J - 1 861 ELSE 862 IF( NSELC.LT.MAX( N / 4, 1 ) ) THEN 863 NSELC = NSELC + 1 864 SELECT( J ) = .TRUE. 865 SELECT( J-1 ) = .FALSE. 866 ELSE 867 SELECT( J ) = .FALSE. 868 SELECT( J-1 ) = .FALSE. 869 END IF 870 J = J - 2 871 END IF 872 IF( J.GT.0 ) 873 $ GO TO 140 874* 875 CALL STREVC( 'Right', 'All', SELECT, N, T1, LDA, DUMMA, LDU, 876 $ EVECTR, LDU, N, IN, WORK, IINFO ) 877 IF( IINFO.NE.0 ) THEN 878 WRITE( NOUNIT, FMT = 9999 )'STREVC(R,A)', IINFO, N, 879 $ JTYPE, IOLDSD 880 INFO = ABS( IINFO ) 881 GO TO 250 882 END IF 883* 884* Test 9: | TR - RW | / ( |T| |R| ulp ) 885* 886 CALL SGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, WR1, 887 $ WI1, WORK, DUMMA( 1 ) ) 888 RESULT( 9 ) = DUMMA( 1 ) 889 IF( DUMMA( 2 ).GT.THRESH ) THEN 890 WRITE( NOUNIT, FMT = 9998 )'Right', 'STREVC', 891 $ DUMMA( 2 ), N, JTYPE, IOLDSD 892 END IF 893* 894* Compute selected right eigenvectors and confirm that 895* they agree with previous right eigenvectors 896* 897 CALL STREVC( 'Right', 'Some', SELECT, N, T1, LDA, DUMMA, 898 $ LDU, EVECTL, LDU, N, IN, WORK, IINFO ) 899 IF( IINFO.NE.0 ) THEN 900 WRITE( NOUNIT, FMT = 9999 )'STREVC(R,S)', IINFO, N, 901 $ JTYPE, IOLDSD 902 INFO = ABS( IINFO ) 903 GO TO 250 904 END IF 905* 906 K = 1 907 MATCH = .TRUE. 908 DO 170 J = 1, N 909 IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN 910 DO 150 JJ = 1, N 911 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN 912 MATCH = .FALSE. 913 GO TO 180 914 END IF 915 150 CONTINUE 916 K = K + 1 917 ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN 918 DO 160 JJ = 1, N 919 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) .OR. 920 $ EVECTR( JJ, J+1 ).NE.EVECTL( JJ, K+1 ) ) THEN 921 MATCH = .FALSE. 922 GO TO 180 923 END IF 924 160 CONTINUE 925 K = K + 2 926 END IF 927 170 CONTINUE 928 180 CONTINUE 929 IF( .NOT.MATCH ) 930 $ WRITE( NOUNIT, FMT = 9997 )'Right', 'STREVC', N, JTYPE, 931 $ IOLDSD 932* 933* Compute the Left eigenvector Matrix: 934* 935 NTEST = 10 936 RESULT( 10 ) = ULPINV 937 CALL STREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU, 938 $ DUMMA, LDU, N, IN, WORK, IINFO ) 939 IF( IINFO.NE.0 ) THEN 940 WRITE( NOUNIT, FMT = 9999 )'STREVC(L,A)', IINFO, N, 941 $ JTYPE, IOLDSD 942 INFO = ABS( IINFO ) 943 GO TO 250 944 END IF 945* 946* Test 10: | LT - WL | / ( |T| |L| ulp ) 947* 948 CALL SGET22( 'Trans', 'N', 'Conj', N, T1, LDA, EVECTL, LDU, 949 $ WR1, WI1, WORK, DUMMA( 3 ) ) 950 RESULT( 10 ) = DUMMA( 3 ) 951 IF( DUMMA( 4 ).GT.THRESH ) THEN 952 WRITE( NOUNIT, FMT = 9998 )'Left', 'STREVC', DUMMA( 4 ), 953 $ N, JTYPE, IOLDSD 954 END IF 955* 956* Compute selected left eigenvectors and confirm that 957* they agree with previous left eigenvectors 958* 959 CALL STREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR, 960 $ LDU, DUMMA, LDU, N, IN, WORK, IINFO ) 961 IF( IINFO.NE.0 ) THEN 962 WRITE( NOUNIT, FMT = 9999 )'STREVC(L,S)', IINFO, N, 963 $ JTYPE, IOLDSD 964 INFO = ABS( IINFO ) 965 GO TO 250 966 END IF 967* 968 K = 1 969 MATCH = .TRUE. 970 DO 210 J = 1, N 971 IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN 972 DO 190 JJ = 1, N 973 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN 974 MATCH = .FALSE. 975 GO TO 220 976 END IF 977 190 CONTINUE 978 K = K + 1 979 ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN 980 DO 200 JJ = 1, N 981 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) .OR. 982 $ EVECTL( JJ, J+1 ).NE.EVECTR( JJ, K+1 ) ) THEN 983 MATCH = .FALSE. 984 GO TO 220 985 END IF 986 200 CONTINUE 987 K = K + 2 988 END IF 989 210 CONTINUE 990 220 CONTINUE 991 IF( .NOT.MATCH ) 992 $ WRITE( NOUNIT, FMT = 9997 )'Left', 'STREVC', N, JTYPE, 993 $ IOLDSD 994* 995* Call SHSEIN for Right eigenvectors of H, do test 11 996* 997 NTEST = 11 998 RESULT( 11 ) = ULPINV 999 DO 230 J = 1, N 1000 SELECT( J ) = .TRUE. 1001 230 CONTINUE 1002* 1003 CALL SHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, 1004 $ WR3, WI3, DUMMA, LDU, EVECTX, LDU, N1, IN, 1005 $ WORK, IWORK, IWORK, IINFO ) 1006 IF( IINFO.NE.0 ) THEN 1007 WRITE( NOUNIT, FMT = 9999 )'SHSEIN(R)', IINFO, N, JTYPE, 1008 $ IOLDSD 1009 INFO = ABS( IINFO ) 1010 IF( IINFO.LT.0 ) 1011 $ GO TO 250 1012 ELSE 1013* 1014* Test 11: | HX - XW | / ( |H| |X| ulp ) 1015* 1016* (from inverse iteration) 1017* 1018 CALL SGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, WR3, 1019 $ WI3, WORK, DUMMA( 1 ) ) 1020 IF( DUMMA( 1 ).LT.ULPINV ) 1021 $ RESULT( 11 ) = DUMMA( 1 )*ANINV 1022 IF( DUMMA( 2 ).GT.THRESH ) THEN 1023 WRITE( NOUNIT, FMT = 9998 )'Right', 'SHSEIN', 1024 $ DUMMA( 2 ), N, JTYPE, IOLDSD 1025 END IF 1026 END IF 1027* 1028* Call SHSEIN for Left eigenvectors of H, do test 12 1029* 1030 NTEST = 12 1031 RESULT( 12 ) = ULPINV 1032 DO 240 J = 1, N 1033 SELECT( J ) = .TRUE. 1034 240 CONTINUE 1035* 1036 CALL SHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, WR3, 1037 $ WI3, EVECTY, LDU, DUMMA, LDU, N1, IN, WORK, 1038 $ IWORK, IWORK, IINFO ) 1039 IF( IINFO.NE.0 ) THEN 1040 WRITE( NOUNIT, FMT = 9999 )'SHSEIN(L)', IINFO, N, JTYPE, 1041 $ IOLDSD 1042 INFO = ABS( IINFO ) 1043 IF( IINFO.LT.0 ) 1044 $ GO TO 250 1045 ELSE 1046* 1047* Test 12: | YH - WY | / ( |H| |Y| ulp ) 1048* 1049* (from inverse iteration) 1050* 1051 CALL SGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, WR3, 1052 $ WI3, WORK, DUMMA( 3 ) ) 1053 IF( DUMMA( 3 ).LT.ULPINV ) 1054 $ RESULT( 12 ) = DUMMA( 3 )*ANINV 1055 IF( DUMMA( 4 ).GT.THRESH ) THEN 1056 WRITE( NOUNIT, FMT = 9998 )'Left', 'SHSEIN', 1057 $ DUMMA( 4 ), N, JTYPE, IOLDSD 1058 END IF 1059 END IF 1060* 1061* Call SORMHR for Right eigenvectors of A, do test 13 1062* 1063 NTEST = 13 1064 RESULT( 13 ) = ULPINV 1065* 1066 CALL SORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 1067 $ LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO ) 1068 IF( IINFO.NE.0 ) THEN 1069 WRITE( NOUNIT, FMT = 9999 )'SORMHR(R)', IINFO, N, JTYPE, 1070 $ IOLDSD 1071 INFO = ABS( IINFO ) 1072 IF( IINFO.LT.0 ) 1073 $ GO TO 250 1074 ELSE 1075* 1076* Test 13: | AX - XW | / ( |A| |X| ulp ) 1077* 1078* (from inverse iteration) 1079* 1080 CALL SGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, WR3, 1081 $ WI3, WORK, DUMMA( 1 ) ) 1082 IF( DUMMA( 1 ).LT.ULPINV ) 1083 $ RESULT( 13 ) = DUMMA( 1 )*ANINV 1084 END IF 1085* 1086* Call SORMHR for Left eigenvectors of A, do test 14 1087* 1088 NTEST = 14 1089 RESULT( 14 ) = ULPINV 1090* 1091 CALL SORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 1092 $ LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO ) 1093 IF( IINFO.NE.0 ) THEN 1094 WRITE( NOUNIT, FMT = 9999 )'SORMHR(L)', IINFO, N, JTYPE, 1095 $ IOLDSD 1096 INFO = ABS( IINFO ) 1097 IF( IINFO.LT.0 ) 1098 $ GO TO 250 1099 ELSE 1100* 1101* Test 14: | YA - WY | / ( |A| |Y| ulp ) 1102* 1103* (from inverse iteration) 1104* 1105 CALL SGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, WR3, 1106 $ WI3, WORK, DUMMA( 3 ) ) 1107 IF( DUMMA( 3 ).LT.ULPINV ) 1108 $ RESULT( 14 ) = DUMMA( 3 )*ANINV 1109 END IF 1110* 1111* End of Loop -- Check for RESULT(j) > THRESH 1112* 1113 250 CONTINUE 1114* 1115 NTESTT = NTESTT + NTEST 1116 CALL SLAFTS( 'SHS', N, N, JTYPE, NTEST, RESULT, IOLDSD, 1117 $ THRESH, NOUNIT, NERRS ) 1118* 1119 260 CONTINUE 1120 270 CONTINUE 1121* 1122* Summary 1123* 1124 CALL SLASUM( 'SHS', NOUNIT, NERRS, NTESTT ) 1125* 1126 RETURN 1127* 1128 9999 FORMAT( ' SCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 1129 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 1130 9998 FORMAT( ' SCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ', 1131 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 1132 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 1133 $ ')' ) 1134 9997 FORMAT( ' SCHKHS: Selected ', A, ' Eigenvectors from ', A, 1135 $ ' do not match other eigenvectors ', 9X, 'N=', I6, 1136 $ ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 1137* 1138* End of SCHKHS 1139* 1140 END 1141