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*> \date November 2015 405* 406*> \ingroup single_eig 407* 408* ===================================================================== 409 SUBROUTINE SCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 410 $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, WR1, 411 $ WI1, WR2, WI2, WR3, WI3, EVECTL, EVECTR, 412 $ EVECTY, EVECTX, UU, TAU, WORK, NWORK, IWORK, 413 $ SELECT, RESULT, INFO ) 414* 415* -- LAPACK test routine (version 3.6.0) -- 416* -- LAPACK is a software package provided by Univ. of Tennessee, -- 417* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 418* November 2015 419* 420* .. Scalar Arguments .. 421 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK 422 REAL THRESH 423* .. 424* .. Array Arguments .. 425 LOGICAL DOTYPE( * ), SELECT( * ) 426 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 427 REAL A( LDA, * ), EVECTL( LDU, * ), 428 $ EVECTR( LDU, * ), EVECTX( LDU, * ), 429 $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), 430 $ T1( LDA, * ), T2( LDA, * ), TAU( * ), 431 $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), 432 $ WI1( * ), WI2( * ), WI3( * ), WORK( * ), 433 $ WR1( * ), WR2( * ), WR3( * ), Z( LDU, * ) 434* .. 435* 436* ===================================================================== 437* 438* .. Parameters .. 439 REAL ZERO, ONE 440 PARAMETER ( ZERO = 0.0, ONE = 1.0 ) 441 INTEGER MAXTYP 442 PARAMETER ( MAXTYP = 21 ) 443* .. 444* .. Local Scalars .. 445 LOGICAL BADNN, MATCH 446 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL, 447 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS, 448 $ NMATS, NMAX, NSELC, NSELR, NTEST, NTESTT 449 REAL ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP, 450 $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL 451* .. 452* .. Local Arrays .. 453 CHARACTER ADUMMA( 1 ) 454 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), 455 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 456 $ KTYPE( MAXTYP ) 457 REAL DUMMA( 6 ) 458* .. 459* .. External Functions .. 460 REAL SLAMCH 461 EXTERNAL SLAMCH 462* .. 463* .. External Subroutines .. 464 EXTERNAL SCOPY, SGEHRD, SGEMM, SGET10, SGET22, SHSEIN, 465 $ SHSEQR, SHST01, SLABAD, SLACPY, SLAFTS, SLASET, 466 $ SLASUM, SLATME, SLATMR, SLATMS, SORGHR, SORMHR, 467 $ STREVC, XERBLA 468* .. 469* .. Intrinsic Functions .. 470 INTRINSIC ABS, MAX, MIN, REAL, SQRT 471* .. 472* .. Data statements .. 473 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 474 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 475 $ 3, 1, 2, 3 / 476 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 477 $ 1, 5, 5, 5, 4, 3, 1 / 478 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 479* .. 480* .. Executable Statements .. 481* 482* Check for errors 483* 484 NTESTT = 0 485 INFO = 0 486* 487 BADNN = .FALSE. 488 NMAX = 0 489 DO 10 J = 1, NSIZES 490 NMAX = MAX( NMAX, NN( J ) ) 491 IF( NN( J ).LT.0 ) 492 $ BADNN = .TRUE. 493 10 CONTINUE 494* 495* Check for errors 496* 497 IF( NSIZES.LT.0 ) THEN 498 INFO = -1 499 ELSE IF( BADNN ) THEN 500 INFO = -2 501 ELSE IF( NTYPES.LT.0 ) THEN 502 INFO = -3 503 ELSE IF( THRESH.LT.ZERO ) THEN 504 INFO = -6 505 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 506 INFO = -9 507 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN 508 INFO = -14 509 ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN 510 INFO = -28 511 END IF 512* 513 IF( INFO.NE.0 ) THEN 514 CALL XERBLA( 'SCHKHS', -INFO ) 515 RETURN 516 END IF 517* 518* Quick return if possible 519* 520 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 521 $ RETURN 522* 523* More important constants 524* 525 UNFL = SLAMCH( 'Safe minimum' ) 526 OVFL = SLAMCH( 'Overflow' ) 527 CALL SLABAD( UNFL, OVFL ) 528 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 529 ULPINV = ONE / ULP 530 RTUNFL = SQRT( UNFL ) 531 RTOVFL = SQRT( OVFL ) 532 RTULP = SQRT( ULP ) 533 RTULPI = ONE / RTULP 534* 535* Loop over sizes, types 536* 537 NERRS = 0 538 NMATS = 0 539* 540 DO 270 JSIZE = 1, NSIZES 541 N = NN( JSIZE ) 542 IF( N.EQ.0 ) 543 $ GO TO 270 544 N1 = MAX( 1, N ) 545 ANINV = ONE / REAL( N1 ) 546* 547 IF( NSIZES.NE.1 ) THEN 548 MTYPES = MIN( MAXTYP, NTYPES ) 549 ELSE 550 MTYPES = MIN( MAXTYP+1, NTYPES ) 551 END IF 552* 553 DO 260 JTYPE = 1, MTYPES 554 IF( .NOT.DOTYPE( JTYPE ) ) 555 $ GO TO 260 556 NMATS = NMATS + 1 557 NTEST = 0 558* 559* Save ISEED in case of an error. 560* 561 DO 20 J = 1, 4 562 IOLDSD( J ) = ISEED( J ) 563 20 CONTINUE 564* 565* Initialize RESULT 566* 567 DO 30 J = 1, 14 568 RESULT( J ) = ZERO 569 30 CONTINUE 570* 571* Compute "A" 572* 573* Control parameters: 574* 575* KMAGN KCONDS KMODE KTYPE 576* =1 O(1) 1 clustered 1 zero 577* =2 large large clustered 2 identity 578* =3 small exponential Jordan 579* =4 arithmetic diagonal, (w/ eigenvalues) 580* =5 random log symmetric, w/ eigenvalues 581* =6 random general, w/ eigenvalues 582* =7 random diagonal 583* =8 random symmetric 584* =9 random general 585* =10 random triangular 586* 587 IF( MTYPES.GT.MAXTYP ) 588 $ GO TO 100 589* 590 ITYPE = KTYPE( JTYPE ) 591 IMODE = KMODE( JTYPE ) 592* 593* Compute norm 594* 595 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 596* 597 40 CONTINUE 598 ANORM = ONE 599 GO TO 70 600* 601 50 CONTINUE 602 ANORM = ( RTOVFL*ULP )*ANINV 603 GO TO 70 604* 605 60 CONTINUE 606 ANORM = RTUNFL*N*ULPINV 607 GO TO 70 608* 609 70 CONTINUE 610* 611 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 612 IINFO = 0 613 COND = ULPINV 614* 615* Special Matrices 616* 617 IF( ITYPE.EQ.1 ) THEN 618* 619* Zero 620* 621 IINFO = 0 622* 623 ELSE IF( ITYPE.EQ.2 ) THEN 624* 625* Identity 626* 627 DO 80 JCOL = 1, N 628 A( JCOL, JCOL ) = ANORM 629 80 CONTINUE 630* 631 ELSE IF( ITYPE.EQ.3 ) THEN 632* 633* Jordan Block 634* 635 DO 90 JCOL = 1, N 636 A( JCOL, JCOL ) = ANORM 637 IF( JCOL.GT.1 ) 638 $ A( JCOL, JCOL-1 ) = ONE 639 90 CONTINUE 640* 641 ELSE IF( ITYPE.EQ.4 ) THEN 642* 643* Diagonal Matrix, [Eigen]values Specified 644* 645 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 646 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 647 $ IINFO ) 648* 649 ELSE IF( ITYPE.EQ.5 ) THEN 650* 651* Symmetric, eigenvalues specified 652* 653 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 654 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 655 $ IINFO ) 656* 657 ELSE IF( ITYPE.EQ.6 ) THEN 658* 659* General, eigenvalues specified 660* 661 IF( KCONDS( JTYPE ).EQ.1 ) THEN 662 CONDS = ONE 663 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 664 CONDS = RTULPI 665 ELSE 666 CONDS = ZERO 667 END IF 668* 669 ADUMMA( 1 ) = ' ' 670 CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 671 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 672 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 673 $ IINFO ) 674* 675 ELSE IF( ITYPE.EQ.7 ) THEN 676* 677* Diagonal, random eigenvalues 678* 679 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 680 $ 'T', 'N', WORK( N+1 ), 1, ONE, 681 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 682 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 683* 684 ELSE IF( ITYPE.EQ.8 ) THEN 685* 686* Symmetric, random eigenvalues 687* 688 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 689 $ 'T', 'N', WORK( N+1 ), 1, ONE, 690 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 691 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 692* 693 ELSE IF( ITYPE.EQ.9 ) THEN 694* 695* General, random eigenvalues 696* 697 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 698 $ 'T', 'N', WORK( N+1 ), 1, ONE, 699 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 700 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 701* 702 ELSE IF( ITYPE.EQ.10 ) THEN 703* 704* Triangular, random eigenvalues 705* 706 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 707 $ 'T', 'N', WORK( N+1 ), 1, ONE, 708 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 709 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 710* 711 ELSE 712* 713 IINFO = 1 714 END IF 715* 716 IF( IINFO.NE.0 ) THEN 717 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 718 $ IOLDSD 719 INFO = ABS( IINFO ) 720 RETURN 721 END IF 722* 723 100 CONTINUE 724* 725* Call SGEHRD to compute H and U, do tests. 726* 727 CALL SLACPY( ' ', N, N, A, LDA, H, LDA ) 728* 729 NTEST = 1 730* 731 ILO = 1 732 IHI = N 733* 734 CALL SGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ), 735 $ NWORK-N, IINFO ) 736* 737 IF( IINFO.NE.0 ) THEN 738 RESULT( 1 ) = ULPINV 739 WRITE( NOUNIT, FMT = 9999 )'SGEHRD', IINFO, N, JTYPE, 740 $ IOLDSD 741 INFO = ABS( IINFO ) 742 GO TO 250 743 END IF 744* 745 DO 120 J = 1, N - 1 746 UU( J+1, J ) = ZERO 747 DO 110 I = J + 2, N 748 U( I, J ) = H( I, J ) 749 UU( I, J ) = H( I, J ) 750 H( I, J ) = ZERO 751 110 CONTINUE 752 120 CONTINUE 753 CALL SCOPY( N-1, WORK, 1, TAU, 1 ) 754 CALL SORGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ), 755 $ NWORK-N, IINFO ) 756 NTEST = 2 757* 758 CALL SHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK, 759 $ NWORK, RESULT( 1 ) ) 760* 761* Call SHSEQR to compute T1, T2 and Z, do tests. 762* 763* Eigenvalues only (WR3,WI3) 764* 765 CALL SLACPY( ' ', N, N, H, LDA, T2, LDA ) 766 NTEST = 3 767 RESULT( 3 ) = ULPINV 768* 769 CALL SHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, WR3, WI3, UZ, 770 $ LDU, WORK, NWORK, IINFO ) 771 IF( IINFO.NE.0 ) THEN 772 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(E)', IINFO, N, JTYPE, 773 $ IOLDSD 774 IF( IINFO.LE.N+2 ) THEN 775 INFO = ABS( IINFO ) 776 GO TO 250 777 END IF 778 END IF 779* 780* Eigenvalues (WR2,WI2) and Full Schur Form (T2) 781* 782 CALL SLACPY( ' ', N, N, H, LDA, T2, LDA ) 783* 784 CALL SHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, WR2, WI2, UZ, 785 $ LDU, WORK, NWORK, IINFO ) 786 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 787 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(S)', IINFO, N, JTYPE, 788 $ IOLDSD 789 INFO = ABS( IINFO ) 790 GO TO 250 791 END IF 792* 793* Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors 794* (UZ) 795* 796 CALL SLACPY( ' ', N, N, H, LDA, T1, LDA ) 797 CALL SLACPY( ' ', N, N, U, LDU, UZ, LDU ) 798* 799 CALL SHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, WR1, WI1, UZ, 800 $ LDU, WORK, NWORK, IINFO ) 801 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 802 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(V)', IINFO, N, JTYPE, 803 $ IOLDSD 804 INFO = ABS( IINFO ) 805 GO TO 250 806 END IF 807* 808* Compute Z = U' UZ 809* 810 CALL SGEMM( 'T', 'N', N, N, N, ONE, U, LDU, UZ, LDU, ZERO, 811 $ Z, LDU ) 812 NTEST = 8 813* 814* Do Tests 3: | H - Z T Z' | / ( |H| n ulp ) 815* and 4: | I - Z Z' | / ( n ulp ) 816* 817 CALL SHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK, 818 $ NWORK, RESULT( 3 ) ) 819* 820* Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp ) 821* and 6: | I - UZ (UZ)' | / ( n ulp ) 822* 823 CALL SHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK, 824 $ NWORK, RESULT( 5 ) ) 825* 826* Do Test 7: | T2 - T1 | / ( |T| n ulp ) 827* 828 CALL SGET10( N, N, T2, LDA, T1, LDA, WORK, RESULT( 7 ) ) 829* 830* Do Test 8: | W2 - W1 | / ( max(|W1|,|W2|) ulp ) 831* 832 TEMP1 = ZERO 833 TEMP2 = ZERO 834 DO 130 J = 1, N 835 TEMP1 = MAX( TEMP1, ABS( WR1( J ) )+ABS( WI1( J ) ), 836 $ ABS( WR2( J ) )+ABS( WI2( J ) ) ) 837 TEMP2 = MAX( TEMP2, ABS( WR1( J )-WR2( J ) )+ 838 $ ABS( WI1( J )-WI2( J ) ) ) 839 130 CONTINUE 840* 841 RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 842* 843* Compute the Left and Right Eigenvectors of T 844* 845* Compute the Right eigenvector Matrix: 846* 847 NTEST = 9 848 RESULT( 9 ) = ULPINV 849* 850* Select last max(N/4,1) real, max(N/4,1) complex eigenvectors 851* 852 NSELC = 0 853 NSELR = 0 854 J = N 855 140 CONTINUE 856 IF( WI1( J ).EQ.ZERO ) THEN 857 IF( NSELR.LT.MAX( N / 4, 1 ) ) THEN 858 NSELR = NSELR + 1 859 SELECT( J ) = .TRUE. 860 ELSE 861 SELECT( J ) = .FALSE. 862 END IF 863 J = J - 1 864 ELSE 865 IF( NSELC.LT.MAX( N / 4, 1 ) ) THEN 866 NSELC = NSELC + 1 867 SELECT( J ) = .TRUE. 868 SELECT( J-1 ) = .FALSE. 869 ELSE 870 SELECT( J ) = .FALSE. 871 SELECT( J-1 ) = .FALSE. 872 END IF 873 J = J - 2 874 END IF 875 IF( J.GT.0 ) 876 $ GO TO 140 877* 878 CALL STREVC( 'Right', 'All', SELECT, N, T1, LDA, DUMMA, LDU, 879 $ EVECTR, LDU, N, IN, WORK, IINFO ) 880 IF( IINFO.NE.0 ) THEN 881 WRITE( NOUNIT, FMT = 9999 )'STREVC(R,A)', IINFO, N, 882 $ JTYPE, IOLDSD 883 INFO = ABS( IINFO ) 884 GO TO 250 885 END IF 886* 887* Test 9: | TR - RW | / ( |T| |R| ulp ) 888* 889 CALL SGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, WR1, 890 $ WI1, WORK, DUMMA( 1 ) ) 891 RESULT( 9 ) = DUMMA( 1 ) 892 IF( DUMMA( 2 ).GT.THRESH ) THEN 893 WRITE( NOUNIT, FMT = 9998 )'Right', 'STREVC', 894 $ DUMMA( 2 ), N, JTYPE, IOLDSD 895 END IF 896* 897* Compute selected right eigenvectors and confirm that 898* they agree with previous right eigenvectors 899* 900 CALL STREVC( 'Right', 'Some', SELECT, N, T1, LDA, DUMMA, 901 $ LDU, EVECTL, LDU, N, IN, WORK, IINFO ) 902 IF( IINFO.NE.0 ) THEN 903 WRITE( NOUNIT, FMT = 9999 )'STREVC(R,S)', IINFO, N, 904 $ JTYPE, IOLDSD 905 INFO = ABS( IINFO ) 906 GO TO 250 907 END IF 908* 909 K = 1 910 MATCH = .TRUE. 911 DO 170 J = 1, N 912 IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN 913 DO 150 JJ = 1, N 914 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN 915 MATCH = .FALSE. 916 GO TO 180 917 END IF 918 150 CONTINUE 919 K = K + 1 920 ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN 921 DO 160 JJ = 1, N 922 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) .OR. 923 $ EVECTR( JJ, J+1 ).NE.EVECTL( JJ, K+1 ) ) THEN 924 MATCH = .FALSE. 925 GO TO 180 926 END IF 927 160 CONTINUE 928 K = K + 2 929 END IF 930 170 CONTINUE 931 180 CONTINUE 932 IF( .NOT.MATCH ) 933 $ WRITE( NOUNIT, FMT = 9997 )'Right', 'STREVC', N, JTYPE, 934 $ IOLDSD 935* 936* Compute the Left eigenvector Matrix: 937* 938 NTEST = 10 939 RESULT( 10 ) = ULPINV 940 CALL STREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU, 941 $ DUMMA, LDU, N, IN, WORK, IINFO ) 942 IF( IINFO.NE.0 ) THEN 943 WRITE( NOUNIT, FMT = 9999 )'STREVC(L,A)', IINFO, N, 944 $ JTYPE, IOLDSD 945 INFO = ABS( IINFO ) 946 GO TO 250 947 END IF 948* 949* Test 10: | LT - WL | / ( |T| |L| ulp ) 950* 951 CALL SGET22( 'Trans', 'N', 'Conj', N, T1, LDA, EVECTL, LDU, 952 $ WR1, WI1, WORK, DUMMA( 3 ) ) 953 RESULT( 10 ) = DUMMA( 3 ) 954 IF( DUMMA( 4 ).GT.THRESH ) THEN 955 WRITE( NOUNIT, FMT = 9998 )'Left', 'STREVC', DUMMA( 4 ), 956 $ N, JTYPE, IOLDSD 957 END IF 958* 959* Compute selected left eigenvectors and confirm that 960* they agree with previous left eigenvectors 961* 962 CALL STREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR, 963 $ LDU, DUMMA, LDU, N, IN, WORK, IINFO ) 964 IF( IINFO.NE.0 ) THEN 965 WRITE( NOUNIT, FMT = 9999 )'STREVC(L,S)', IINFO, N, 966 $ JTYPE, IOLDSD 967 INFO = ABS( IINFO ) 968 GO TO 250 969 END IF 970* 971 K = 1 972 MATCH = .TRUE. 973 DO 210 J = 1, N 974 IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN 975 DO 190 JJ = 1, N 976 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN 977 MATCH = .FALSE. 978 GO TO 220 979 END IF 980 190 CONTINUE 981 K = K + 1 982 ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN 983 DO 200 JJ = 1, N 984 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) .OR. 985 $ EVECTL( JJ, J+1 ).NE.EVECTR( JJ, K+1 ) ) THEN 986 MATCH = .FALSE. 987 GO TO 220 988 END IF 989 200 CONTINUE 990 K = K + 2 991 END IF 992 210 CONTINUE 993 220 CONTINUE 994 IF( .NOT.MATCH ) 995 $ WRITE( NOUNIT, FMT = 9997 )'Left', 'STREVC', N, JTYPE, 996 $ IOLDSD 997* 998* Call SHSEIN for Right eigenvectors of H, do test 11 999* 1000 NTEST = 11 1001 RESULT( 11 ) = ULPINV 1002 DO 230 J = 1, N 1003 SELECT( J ) = .TRUE. 1004 230 CONTINUE 1005* 1006 CALL SHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, 1007 $ WR3, WI3, DUMMA, LDU, EVECTX, LDU, N1, IN, 1008 $ WORK, IWORK, IWORK, IINFO ) 1009 IF( IINFO.NE.0 ) THEN 1010 WRITE( NOUNIT, FMT = 9999 )'SHSEIN(R)', IINFO, N, JTYPE, 1011 $ IOLDSD 1012 INFO = ABS( IINFO ) 1013 IF( IINFO.LT.0 ) 1014 $ GO TO 250 1015 ELSE 1016* 1017* Test 11: | HX - XW | / ( |H| |X| ulp ) 1018* 1019* (from inverse iteration) 1020* 1021 CALL SGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, WR3, 1022 $ WI3, WORK, DUMMA( 1 ) ) 1023 IF( DUMMA( 1 ).LT.ULPINV ) 1024 $ RESULT( 11 ) = DUMMA( 1 )*ANINV 1025 IF( DUMMA( 2 ).GT.THRESH ) THEN 1026 WRITE( NOUNIT, FMT = 9998 )'Right', 'SHSEIN', 1027 $ DUMMA( 2 ), N, JTYPE, IOLDSD 1028 END IF 1029 END IF 1030* 1031* Call SHSEIN for Left eigenvectors of H, do test 12 1032* 1033 NTEST = 12 1034 RESULT( 12 ) = ULPINV 1035 DO 240 J = 1, N 1036 SELECT( J ) = .TRUE. 1037 240 CONTINUE 1038* 1039 CALL SHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, WR3, 1040 $ WI3, EVECTY, LDU, DUMMA, LDU, N1, IN, WORK, 1041 $ IWORK, IWORK, IINFO ) 1042 IF( IINFO.NE.0 ) THEN 1043 WRITE( NOUNIT, FMT = 9999 )'SHSEIN(L)', IINFO, N, JTYPE, 1044 $ IOLDSD 1045 INFO = ABS( IINFO ) 1046 IF( IINFO.LT.0 ) 1047 $ GO TO 250 1048 ELSE 1049* 1050* Test 12: | YH - WY | / ( |H| |Y| ulp ) 1051* 1052* (from inverse iteration) 1053* 1054 CALL SGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, WR3, 1055 $ WI3, WORK, DUMMA( 3 ) ) 1056 IF( DUMMA( 3 ).LT.ULPINV ) 1057 $ RESULT( 12 ) = DUMMA( 3 )*ANINV 1058 IF( DUMMA( 4 ).GT.THRESH ) THEN 1059 WRITE( NOUNIT, FMT = 9998 )'Left', 'SHSEIN', 1060 $ DUMMA( 4 ), N, JTYPE, IOLDSD 1061 END IF 1062 END IF 1063* 1064* Call SORMHR for Right eigenvectors of A, do test 13 1065* 1066 NTEST = 13 1067 RESULT( 13 ) = ULPINV 1068* 1069 CALL SORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 1070 $ LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO ) 1071 IF( IINFO.NE.0 ) THEN 1072 WRITE( NOUNIT, FMT = 9999 )'SORMHR(R)', IINFO, N, JTYPE, 1073 $ IOLDSD 1074 INFO = ABS( IINFO ) 1075 IF( IINFO.LT.0 ) 1076 $ GO TO 250 1077 ELSE 1078* 1079* Test 13: | AX - XW | / ( |A| |X| ulp ) 1080* 1081* (from inverse iteration) 1082* 1083 CALL SGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, WR3, 1084 $ WI3, WORK, DUMMA( 1 ) ) 1085 IF( DUMMA( 1 ).LT.ULPINV ) 1086 $ RESULT( 13 ) = DUMMA( 1 )*ANINV 1087 END IF 1088* 1089* Call SORMHR for Left eigenvectors of A, do test 14 1090* 1091 NTEST = 14 1092 RESULT( 14 ) = ULPINV 1093* 1094 CALL SORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 1095 $ LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO ) 1096 IF( IINFO.NE.0 ) THEN 1097 WRITE( NOUNIT, FMT = 9999 )'SORMHR(L)', IINFO, N, JTYPE, 1098 $ IOLDSD 1099 INFO = ABS( IINFO ) 1100 IF( IINFO.LT.0 ) 1101 $ GO TO 250 1102 ELSE 1103* 1104* Test 14: | YA - WY | / ( |A| |Y| ulp ) 1105* 1106* (from inverse iteration) 1107* 1108 CALL SGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, WR3, 1109 $ WI3, WORK, DUMMA( 3 ) ) 1110 IF( DUMMA( 3 ).LT.ULPINV ) 1111 $ RESULT( 14 ) = DUMMA( 3 )*ANINV 1112 END IF 1113* 1114* End of Loop -- Check for RESULT(j) > THRESH 1115* 1116 250 CONTINUE 1117* 1118 NTESTT = NTESTT + NTEST 1119 CALL SLAFTS( 'SHS', N, N, JTYPE, NTEST, RESULT, IOLDSD, 1120 $ THRESH, NOUNIT, NERRS ) 1121* 1122 260 CONTINUE 1123 270 CONTINUE 1124* 1125* Summary 1126* 1127 CALL SLASUM( 'SHS', NOUNIT, NERRS, NTESTT ) 1128* 1129 RETURN 1130* 1131 9999 FORMAT( ' SCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 1132 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 1133 9998 FORMAT( ' SCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ', 1134 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 1135 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 1136 $ ')' ) 1137 9997 FORMAT( ' SCHKHS: Selected ', A, ' Eigenvectors from ', A, 1138 $ ' do not match other eigenvectors ', 9X, 'N=', I6, 1139 $ ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 1140* 1141* End of SCHKHS 1142* 1143 END 1144