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, WR3, WI3, EVECTL, EVECTR, EVECTY, EVECTX, 14* UU, TAU, WORK, NWORK, IWORK, SELECT, RESULT, 15* 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( * ), WI3( * ), WORK( * ), WR1( * ), 30* $ 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*> WR3 - REAL array, dimension (max(NN)) 272*> WI3 - REAL array, dimension (max(NN)) 273*> Like WR1, WI1, these arrays contain the eigenvalues of A, 274*> but those computed when SHSEQR only computes the 275*> eigenvalues, i.e., not the Schur vectors and no more of the 276*> Schur form than is necessary for computing the 277*> eigenvalues. 278*> Modified. 279*> 280*> EVECTL - REAL array, dimension (LDU,max(NN)) 281*> The (upper triangular) left eigenvector matrix for the 282*> matrix in T1. For complex conjugate pairs, the real part 283*> is stored in one row and the imaginary part in the next. 284*> Modified. 285*> 286*> EVECTR - REAL array, dimension (LDU,max(NN)) 287*> The (upper triangular) right eigenvector matrix for the 288*> matrix in T1. For complex conjugate pairs, the real part 289*> is stored in one column and the imaginary part in the next. 290*> Modified. 291*> 292*> EVECTY - REAL array, dimension (LDU,max(NN)) 293*> The left eigenvector matrix for the 294*> matrix in H. For complex conjugate pairs, the real part 295*> is stored in one row and the imaginary part in the next. 296*> Modified. 297*> 298*> EVECTX - REAL array, dimension (LDU,max(NN)) 299*> The right eigenvector matrix for the 300*> matrix in H. For complex conjugate pairs, the real part 301*> is stored in one column and the imaginary part in the next. 302*> Modified. 303*> 304*> UU - REAL array, dimension (LDU,max(NN)) 305*> Details of the orthogonal matrix computed by SGEHRD. 306*> Modified. 307*> 308*> TAU - REAL array, dimension(max(NN)) 309*> Further details of the orthogonal matrix computed by SGEHRD. 310*> Modified. 311*> 312*> WORK - REAL array, dimension (NWORK) 313*> Workspace. 314*> Modified. 315*> 316*> NWORK - INTEGER 317*> The number of entries in WORK. NWORK >= 4*NN(j)*NN(j) + 2. 318*> 319*> IWORK - INTEGER array, dimension (max(NN)) 320*> Workspace. 321*> Modified. 322*> 323*> SELECT - LOGICAL array, dimension (max(NN)) 324*> Workspace. 325*> Modified. 326*> 327*> RESULT - REAL array, dimension (14) 328*> The values computed by the fourteen tests described above. 329*> The values are currently limited to 1/ulp, to avoid 330*> overflow. 331*> Modified. 332*> 333*> INFO - INTEGER 334*> If 0, then everything ran OK. 335*> -1: NSIZES < 0 336*> -2: Some NN(j) < 0 337*> -3: NTYPES < 0 338*> -6: THRESH < 0 339*> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). 340*> -14: LDU < 1 or LDU < NMAX. 341*> -28: NWORK too small. 342*> If SLATMR, SLATMS, or SLATME returns an error code, the 343*> absolute value of it is returned. 344*> If 1, then SHSEQR could not find all the shifts. 345*> If 2, then the EISPACK code (for small blocks) failed. 346*> If >2, then 30*N iterations were not enough to find an 347*> eigenvalue or to decompose the problem. 348*> Modified. 349*> 350*>----------------------------------------------------------------------- 351*> 352*> Some Local Variables and Parameters: 353*> ---- ----- --------- --- ---------- 354*> 355*> ZERO, ONE Real 0 and 1. 356*> MAXTYP The number of types defined. 357*> MTEST The number of tests defined: care must be taken 358*> that (1) the size of RESULT, (2) the number of 359*> tests actually performed, and (3) MTEST agree. 360*> NTEST The number of tests performed on this matrix 361*> so far. This should be less than MTEST, and 362*> equal to it by the last test. It will be less 363*> if any of the routines being tested indicates 364*> that it could not compute the matrices that 365*> would be tested. 366*> NMAX Largest value in NN. 367*> NMATS The number of matrices generated so far. 368*> NERRS The number of tests which have exceeded THRESH 369*> so far (computed by SLAFTS). 370*> COND, CONDS, 371*> IMODE Values to be passed to the matrix generators. 372*> ANORM Norm of A; passed to matrix generators. 373*> 374*> OVFL, UNFL Overflow and underflow thresholds. 375*> ULP, ULPINV Finest relative precision and its inverse. 376*> RTOVFL, RTUNFL, 377*> RTULP, RTULPI Square roots of the previous 4 values. 378*> 379*> The following four arrays decode JTYPE: 380*> KTYPE(j) The general type (1-10) for type "j". 381*> KMODE(j) The MODE value to be passed to the matrix 382*> generator for type "j". 383*> KMAGN(j) The order of magnitude ( O(1), 384*> O(overflow^(1/2) ), O(underflow^(1/2) ) 385*> KCONDS(j) Selects whether CONDS is to be 1 or 386*> 1/sqrt(ulp). (0 means irrelevant.) 387*> \endverbatim 388* 389* Authors: 390* ======== 391* 392*> \author Univ. of Tennessee 393*> \author Univ. of California Berkeley 394*> \author Univ. of Colorado Denver 395*> \author NAG Ltd. 396* 397*> \date November 2011 398* 399*> \ingroup single_eig 400* 401* ===================================================================== 402 SUBROUTINE SCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 403 $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, WR1, 404 $ WI1, WR3, WI3, EVECTL, EVECTR, EVECTY, EVECTX, 405 $ UU, TAU, WORK, NWORK, IWORK, SELECT, RESULT, 406 $ INFO ) 407* 408* -- LAPACK test routine (version 3.4.0) -- 409* -- LAPACK is a software package provided by Univ. of Tennessee, -- 410* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 411* November 2011 412* 413* .. Scalar Arguments .. 414 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK 415 REAL THRESH 416* .. 417* .. Array Arguments .. 418 LOGICAL DOTYPE( * ), SELECT( * ) 419 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 420 REAL A( LDA, * ), EVECTL( LDU, * ), 421 $ EVECTR( LDU, * ), EVECTX( LDU, * ), 422 $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), 423 $ T1( LDA, * ), T2( LDA, * ), TAU( * ), 424 $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), 425 $ WI1( * ), WI3( * ), WORK( * ), WR1( * ), 426 $ WR3( * ), Z( LDU, * ) 427* .. 428* 429* ===================================================================== 430* 431* .. Parameters .. 432 REAL ZERO, ONE 433 PARAMETER ( ZERO = 0.0, ONE = 1.0 ) 434 INTEGER MAXTYP 435 PARAMETER ( MAXTYP = 21 ) 436* .. 437* .. Local Scalars .. 438 LOGICAL BADNN, MATCH 439 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL, 440 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS, 441 $ NMATS, NMAX, NSELC, NSELR, NTEST, NTESTT 442 REAL ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP, 443 $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL 444* .. 445* .. Local Arrays .. 446 CHARACTER ADUMMA( 1 ) 447 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), 448 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 449 $ KTYPE( MAXTYP ) 450 REAL DUMMA( 6 ) 451* .. 452* .. External Functions .. 453 REAL SLAMCH 454 EXTERNAL SLAMCH 455* .. 456* .. External Subroutines .. 457 EXTERNAL SCOPY, SGEHRD, SGEMM, SGET10, SGET22, SHSEIN, 458 $ SHSEQR, SHST01, SLABAD, SLACPY, SLAFTS, SLASET, 459 $ SLASUM, SLATME, SLATMR, SLATMS, SORGHR, SORMHR, 460 $ STREVC, XERBLA 461* .. 462* .. Intrinsic Functions .. 463 INTRINSIC ABS, MAX, MIN, REAL, SQRT 464* .. 465* .. Data statements .. 466 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 467 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 468 $ 3, 1, 2, 3 / 469 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 470 $ 1, 5, 5, 5, 4, 3, 1 / 471 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 472* .. 473* .. Executable Statements .. 474* 475* Check for errors 476* 477 NTESTT = 0 478 INFO = 0 479* 480 BADNN = .FALSE. 481 NMAX = 0 482 DO 10 J = 1, NSIZES 483 NMAX = MAX( NMAX, NN( J ) ) 484 IF( NN( J ).LT.0 ) 485 $ BADNN = .TRUE. 486 10 CONTINUE 487* 488* Check for errors 489* 490 IF( NSIZES.LT.0 ) THEN 491 INFO = -1 492 ELSE IF( BADNN ) THEN 493 INFO = -2 494 ELSE IF( NTYPES.LT.0 ) THEN 495 INFO = -3 496 ELSE IF( THRESH.LT.ZERO ) THEN 497 INFO = -6 498 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN 499 INFO = -9 500 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN 501 INFO = -14 502 ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN 503 INFO = -28 504 END IF 505* 506 IF( INFO.NE.0 ) THEN 507 CALL XERBLA( 'SCHKHS', -INFO ) 508 RETURN 509 END IF 510* 511* Quick return if possible 512* 513 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 514 $ RETURN 515* 516* More important constants 517* 518 UNFL = SLAMCH( 'Safe minimum' ) 519 OVFL = SLAMCH( 'Overflow' ) 520 CALL SLABAD( UNFL, OVFL ) 521 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 522 ULPINV = ONE / ULP 523 RTUNFL = SQRT( UNFL ) 524 RTOVFL = SQRT( OVFL ) 525 RTULP = SQRT( ULP ) 526 RTULPI = ONE / RTULP 527* 528* Loop over sizes, types 529* 530 NERRS = 0 531 NMATS = 0 532* 533 DO 270 JSIZE = 1, NSIZES 534 N = NN( JSIZE ) 535 IF( N.EQ.0 ) 536 $ GO TO 270 537 N1 = MAX( 1, N ) 538 ANINV = ONE / REAL( N1 ) 539* 540 IF( NSIZES.NE.1 ) THEN 541 MTYPES = MIN( MAXTYP, NTYPES ) 542 ELSE 543 MTYPES = MIN( MAXTYP+1, NTYPES ) 544 END IF 545* 546 DO 260 JTYPE = 1, MTYPES 547 IF( .NOT.DOTYPE( JTYPE ) ) 548 $ GO TO 260 549 NMATS = NMATS + 1 550 NTEST = 0 551* 552* Save ISEED in case of an error. 553* 554 DO 20 J = 1, 4 555 IOLDSD( J ) = ISEED( J ) 556 20 CONTINUE 557* 558* Initialize RESULT 559* 560 DO 30 J = 1, 14 561 RESULT( J ) = ZERO 562 30 CONTINUE 563* 564* Compute "A" 565* 566* Control parameters: 567* 568* KMAGN KCONDS KMODE KTYPE 569* =1 O(1) 1 clustered 1 zero 570* =2 large large clustered 2 identity 571* =3 small exponential Jordan 572* =4 arithmetic diagonal, (w/ eigenvalues) 573* =5 random log symmetric, w/ eigenvalues 574* =6 random general, w/ eigenvalues 575* =7 random diagonal 576* =8 random symmetric 577* =9 random general 578* =10 random triangular 579* 580 IF( MTYPES.GT.MAXTYP ) 581 $ GO TO 100 582* 583 ITYPE = KTYPE( JTYPE ) 584 IMODE = KMODE( JTYPE ) 585* 586* Compute norm 587* 588 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 589* 590 40 CONTINUE 591 ANORM = ONE 592 GO TO 70 593* 594 50 CONTINUE 595 ANORM = ( RTOVFL*ULP )*ANINV 596 GO TO 70 597* 598 60 CONTINUE 599 ANORM = RTUNFL*N*ULPINV 600 GO TO 70 601* 602 70 CONTINUE 603* 604 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 605 IINFO = 0 606 COND = ULPINV 607* 608* Special Matrices 609* 610 IF( ITYPE.EQ.1 ) THEN 611* 612* Zero 613* 614 IINFO = 0 615* 616 ELSE IF( ITYPE.EQ.2 ) THEN 617* 618* Identity 619* 620 DO 80 JCOL = 1, N 621 A( JCOL, JCOL ) = ANORM 622 80 CONTINUE 623* 624 ELSE IF( ITYPE.EQ.3 ) THEN 625* 626* Jordan Block 627* 628 DO 90 JCOL = 1, N 629 A( JCOL, JCOL ) = ANORM 630 IF( JCOL.GT.1 ) 631 $ A( JCOL, JCOL-1 ) = ONE 632 90 CONTINUE 633* 634 ELSE IF( ITYPE.EQ.4 ) THEN 635* 636* Diagonal Matrix, [Eigen]values Specified 637* 638 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 639 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 640 $ IINFO ) 641* 642 ELSE IF( ITYPE.EQ.5 ) THEN 643* 644* Symmetric, eigenvalues specified 645* 646 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 647 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 648 $ IINFO ) 649* 650 ELSE IF( ITYPE.EQ.6 ) THEN 651* 652* General, eigenvalues specified 653* 654 IF( KCONDS( JTYPE ).EQ.1 ) THEN 655 CONDS = ONE 656 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 657 CONDS = RTULPI 658 ELSE 659 CONDS = ZERO 660 END IF 661* 662 ADUMMA( 1 ) = ' ' 663 CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 664 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 665 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 666 $ IINFO ) 667* 668 ELSE IF( ITYPE.EQ.7 ) THEN 669* 670* Diagonal, random eigenvalues 671* 672 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 673 $ 'T', 'N', WORK( N+1 ), 1, ONE, 674 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 675 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 676* 677 ELSE IF( ITYPE.EQ.8 ) THEN 678* 679* Symmetric, random eigenvalues 680* 681 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 682 $ 'T', 'N', WORK( N+1 ), 1, ONE, 683 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 684 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 685* 686 ELSE IF( ITYPE.EQ.9 ) THEN 687* 688* General, random eigenvalues 689* 690 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 691 $ 'T', 'N', WORK( N+1 ), 1, ONE, 692 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 693 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 694* 695 ELSE IF( ITYPE.EQ.10 ) THEN 696* 697* Triangular, random eigenvalues 698* 699 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 700 $ 'T', 'N', WORK( N+1 ), 1, ONE, 701 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 702 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 703* 704 ELSE 705* 706 IINFO = 1 707 END IF 708* 709 IF( IINFO.NE.0 ) THEN 710 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 711 $ IOLDSD 712 INFO = ABS( IINFO ) 713 RETURN 714 END IF 715* 716 100 CONTINUE 717* 718* Call SGEHRD to compute H and U, do tests. 719* 720 CALL SLACPY( ' ', N, N, A, LDA, H, LDA ) 721* 722 NTEST = 1 723* 724 ILO = 1 725 IHI = N 726* 727 CALL SGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ), 728 $ NWORK-N, IINFO ) 729* 730 IF( IINFO.NE.0 ) THEN 731 RESULT( 1 ) = ULPINV 732 WRITE( NOUNIT, FMT = 9999 )'SGEHRD', IINFO, N, JTYPE, 733 $ IOLDSD 734 INFO = ABS( IINFO ) 735 GO TO 250 736 END IF 737* 738 DO 120 J = 1, N - 1 739 UU( J+1, J ) = ZERO 740 DO 110 I = J + 2, N 741 U( I, J ) = H( I, J ) 742 UU( I, J ) = H( I, J ) 743 H( I, J ) = ZERO 744 110 CONTINUE 745 120 CONTINUE 746 CALL SCOPY( N-1, WORK, 1, TAU, 1 ) 747 CALL SORGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ), 748 $ NWORK-N, IINFO ) 749 NTEST = 2 750* 751 CALL SHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK, 752 $ NWORK, RESULT( 1 ) ) 753* 754* Call SHSEQR to compute T1, T2 and Z, do tests. 755* 756* Eigenvalues only (WR3,WI3) 757* 758 CALL SLACPY( ' ', N, N, H, LDA, T2, LDA ) 759 NTEST = 3 760 RESULT( 3 ) = ULPINV 761* 762 CALL SHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, WR3, WI3, UZ, 763 $ LDU, WORK, NWORK, IINFO ) 764 IF( IINFO.NE.0 ) THEN 765 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(E)', IINFO, N, JTYPE, 766 $ IOLDSD 767 IF( IINFO.LE.N+2 ) THEN 768 INFO = ABS( IINFO ) 769 GO TO 250 770 END IF 771 END IF 772* 773* Eigenvalues (WR1,WI1) and Full Schur Form (T2) 774* 775 CALL SLACPY( ' ', N, N, H, LDA, T2, LDA ) 776* 777 CALL SHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, WR1, WI1, UZ, 778 $ LDU, WORK, NWORK, IINFO ) 779 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 780 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(S)', IINFO, N, JTYPE, 781 $ IOLDSD 782 INFO = ABS( IINFO ) 783 GO TO 250 784 END IF 785* 786* Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors 787* (UZ) 788* 789 CALL SLACPY( ' ', N, N, H, LDA, T1, LDA ) 790 CALL SLACPY( ' ', N, N, U, LDU, UZ, LDA ) 791* 792 CALL SHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, WR1, WI1, UZ, 793 $ LDU, WORK, NWORK, IINFO ) 794 IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN 795 WRITE( NOUNIT, FMT = 9999 )'SHSEQR(V)', IINFO, N, JTYPE, 796 $ IOLDSD 797 INFO = ABS( IINFO ) 798 GO TO 250 799 END IF 800* 801* Compute Z = U' UZ 802* 803 CALL SGEMM( 'T', 'N', N, N, N, ONE, U, LDU, UZ, LDU, ZERO, 804 $ Z, LDU ) 805 NTEST = 8 806* 807* Do Tests 3: | H - Z T Z' | / ( |H| n ulp ) 808* and 4: | I - Z Z' | / ( n ulp ) 809* 810 CALL SHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK, 811 $ NWORK, RESULT( 3 ) ) 812* 813* Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp ) 814* and 6: | I - UZ (UZ)' | / ( n ulp ) 815* 816 CALL SHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK, 817 $ NWORK, RESULT( 5 ) ) 818* 819* Do Test 7: | T2 - T1 | / ( |T| n ulp ) 820* 821 CALL SGET10( N, N, T2, LDA, T1, LDA, WORK, RESULT( 7 ) ) 822* 823* Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp ) 824* 825 TEMP1 = ZERO 826 TEMP2 = ZERO 827 DO 130 J = 1, N 828 TEMP1 = MAX( TEMP1, ABS( WR1( J ) )+ABS( WI1( J ) ), 829 $ ABS( WR3( J ) )+ABS( WI3( J ) ) ) 830 TEMP2 = MAX( TEMP2, ABS( WR1( J )-WR3( J ) )+ 831 $ ABS( WR1( J )-WR3( J ) ) ) 832 130 CONTINUE 833* 834 RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 835* 836* Compute the Left and Right Eigenvectors of T 837* 838* Compute the Right eigenvector Matrix: 839* 840 NTEST = 9 841 RESULT( 9 ) = ULPINV 842* 843* Select last max(N/4,1) real, max(N/4,1) complex eigenvectors 844* 845 NSELC = 0 846 NSELR = 0 847 J = N 848 140 CONTINUE 849 IF( WI1( J ).EQ.ZERO ) THEN 850 IF( NSELR.LT.MAX( N / 4, 1 ) ) THEN 851 NSELR = NSELR + 1 852 SELECT( J ) = .TRUE. 853 ELSE 854 SELECT( J ) = .FALSE. 855 END IF 856 J = J - 1 857 ELSE 858 IF( NSELC.LT.MAX( N / 4, 1 ) ) THEN 859 NSELC = NSELC + 1 860 SELECT( J ) = .TRUE. 861 SELECT( J-1 ) = .FALSE. 862 ELSE 863 SELECT( J ) = .FALSE. 864 SELECT( J-1 ) = .FALSE. 865 END IF 866 J = J - 2 867 END IF 868 IF( J.GT.0 ) 869 $ GO TO 140 870* 871 CALL STREVC( 'Right', 'All', SELECT, N, T1, LDA, DUMMA, LDU, 872 $ EVECTR, LDU, N, IN, WORK, IINFO ) 873 IF( IINFO.NE.0 ) THEN 874 WRITE( NOUNIT, FMT = 9999 )'STREVC(R,A)', IINFO, N, 875 $ JTYPE, IOLDSD 876 INFO = ABS( IINFO ) 877 GO TO 250 878 END IF 879* 880* Test 9: | TR - RW | / ( |T| |R| ulp ) 881* 882 CALL SGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, WR1, 883 $ WI1, WORK, DUMMA( 1 ) ) 884 RESULT( 9 ) = DUMMA( 1 ) 885 IF( DUMMA( 2 ).GT.THRESH ) THEN 886 WRITE( NOUNIT, FMT = 9998 )'Right', 'STREVC', 887 $ DUMMA( 2 ), N, JTYPE, IOLDSD 888 END IF 889* 890* Compute selected right eigenvectors and confirm that 891* they agree with previous right eigenvectors 892* 893 CALL STREVC( 'Right', 'Some', SELECT, N, T1, LDA, DUMMA, 894 $ LDU, EVECTL, LDU, N, IN, WORK, IINFO ) 895 IF( IINFO.NE.0 ) THEN 896 WRITE( NOUNIT, FMT = 9999 )'STREVC(R,S)', IINFO, N, 897 $ JTYPE, IOLDSD 898 INFO = ABS( IINFO ) 899 GO TO 250 900 END IF 901* 902 K = 1 903 MATCH = .TRUE. 904 DO 170 J = 1, N 905 IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN 906 DO 150 JJ = 1, N 907 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN 908 MATCH = .FALSE. 909 GO TO 180 910 END IF 911 150 CONTINUE 912 K = K + 1 913 ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN 914 DO 160 JJ = 1, N 915 IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) .OR. 916 $ EVECTR( JJ, J+1 ).NE.EVECTL( JJ, K+1 ) ) THEN 917 MATCH = .FALSE. 918 GO TO 180 919 END IF 920 160 CONTINUE 921 K = K + 2 922 END IF 923 170 CONTINUE 924 180 CONTINUE 925 IF( .NOT.MATCH ) 926 $ WRITE( NOUNIT, FMT = 9997 )'Right', 'STREVC', N, JTYPE, 927 $ IOLDSD 928* 929* Compute the Left eigenvector Matrix: 930* 931 NTEST = 10 932 RESULT( 10 ) = ULPINV 933 CALL STREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU, 934 $ DUMMA, LDU, N, IN, WORK, IINFO ) 935 IF( IINFO.NE.0 ) THEN 936 WRITE( NOUNIT, FMT = 9999 )'STREVC(L,A)', IINFO, N, 937 $ JTYPE, IOLDSD 938 INFO = ABS( IINFO ) 939 GO TO 250 940 END IF 941* 942* Test 10: | LT - WL | / ( |T| |L| ulp ) 943* 944 CALL SGET22( 'Trans', 'N', 'Conj', N, T1, LDA, EVECTL, LDU, 945 $ WR1, WI1, WORK, DUMMA( 3 ) ) 946 RESULT( 10 ) = DUMMA( 3 ) 947 IF( DUMMA( 4 ).GT.THRESH ) THEN 948 WRITE( NOUNIT, FMT = 9998 )'Left', 'STREVC', DUMMA( 4 ), 949 $ N, JTYPE, IOLDSD 950 END IF 951* 952* Compute selected left eigenvectors and confirm that 953* they agree with previous left eigenvectors 954* 955 CALL STREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR, 956 $ LDU, DUMMA, LDU, N, IN, WORK, IINFO ) 957 IF( IINFO.NE.0 ) THEN 958 WRITE( NOUNIT, FMT = 9999 )'STREVC(L,S)', IINFO, N, 959 $ JTYPE, IOLDSD 960 INFO = ABS( IINFO ) 961 GO TO 250 962 END IF 963* 964 K = 1 965 MATCH = .TRUE. 966 DO 210 J = 1, N 967 IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN 968 DO 190 JJ = 1, N 969 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN 970 MATCH = .FALSE. 971 GO TO 220 972 END IF 973 190 CONTINUE 974 K = K + 1 975 ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN 976 DO 200 JJ = 1, N 977 IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) .OR. 978 $ EVECTL( JJ, J+1 ).NE.EVECTR( JJ, K+1 ) ) THEN 979 MATCH = .FALSE. 980 GO TO 220 981 END IF 982 200 CONTINUE 983 K = K + 2 984 END IF 985 210 CONTINUE 986 220 CONTINUE 987 IF( .NOT.MATCH ) 988 $ WRITE( NOUNIT, FMT = 9997 )'Left', 'STREVC', N, JTYPE, 989 $ IOLDSD 990* 991* Call SHSEIN for Right eigenvectors of H, do test 11 992* 993 NTEST = 11 994 RESULT( 11 ) = ULPINV 995 DO 230 J = 1, N 996 SELECT( J ) = .TRUE. 997 230 CONTINUE 998* 999 CALL SHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, 1000 $ WR3, WI3, DUMMA, LDU, EVECTX, LDU, N1, IN, 1001 $ WORK, IWORK, IWORK, IINFO ) 1002 IF( IINFO.NE.0 ) THEN 1003 WRITE( NOUNIT, FMT = 9999 )'SHSEIN(R)', IINFO, N, JTYPE, 1004 $ IOLDSD 1005 INFO = ABS( IINFO ) 1006 IF( IINFO.LT.0 ) 1007 $ GO TO 250 1008 ELSE 1009* 1010* Test 11: | HX - XW | / ( |H| |X| ulp ) 1011* 1012* (from inverse iteration) 1013* 1014 CALL SGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, WR3, 1015 $ WI3, WORK, DUMMA( 1 ) ) 1016 IF( DUMMA( 1 ).LT.ULPINV ) 1017 $ RESULT( 11 ) = DUMMA( 1 )*ANINV 1018 IF( DUMMA( 2 ).GT.THRESH ) THEN 1019 WRITE( NOUNIT, FMT = 9998 )'Right', 'SHSEIN', 1020 $ DUMMA( 2 ), N, JTYPE, IOLDSD 1021 END IF 1022 END IF 1023* 1024* Call SHSEIN for Left eigenvectors of H, do test 12 1025* 1026 NTEST = 12 1027 RESULT( 12 ) = ULPINV 1028 DO 240 J = 1, N 1029 SELECT( J ) = .TRUE. 1030 240 CONTINUE 1031* 1032 CALL SHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, WR3, 1033 $ WI3, EVECTY, LDU, DUMMA, LDU, N1, IN, WORK, 1034 $ IWORK, IWORK, IINFO ) 1035 IF( IINFO.NE.0 ) THEN 1036 WRITE( NOUNIT, FMT = 9999 )'SHSEIN(L)', IINFO, N, JTYPE, 1037 $ IOLDSD 1038 INFO = ABS( IINFO ) 1039 IF( IINFO.LT.0 ) 1040 $ GO TO 250 1041 ELSE 1042* 1043* Test 12: | YH - WY | / ( |H| |Y| ulp ) 1044* 1045* (from inverse iteration) 1046* 1047 CALL SGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, WR3, 1048 $ WI3, WORK, DUMMA( 3 ) ) 1049 IF( DUMMA( 3 ).LT.ULPINV ) 1050 $ RESULT( 12 ) = DUMMA( 3 )*ANINV 1051 IF( DUMMA( 4 ).GT.THRESH ) THEN 1052 WRITE( NOUNIT, FMT = 9998 )'Left', 'SHSEIN', 1053 $ DUMMA( 4 ), N, JTYPE, IOLDSD 1054 END IF 1055 END IF 1056* 1057* Call SORMHR for Right eigenvectors of A, do test 13 1058* 1059 NTEST = 13 1060 RESULT( 13 ) = ULPINV 1061* 1062 CALL SORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 1063 $ LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO ) 1064 IF( IINFO.NE.0 ) THEN 1065 WRITE( NOUNIT, FMT = 9999 )'SORMHR(R)', IINFO, N, JTYPE, 1066 $ IOLDSD 1067 INFO = ABS( IINFO ) 1068 IF( IINFO.LT.0 ) 1069 $ GO TO 250 1070 ELSE 1071* 1072* Test 13: | AX - XW | / ( |A| |X| ulp ) 1073* 1074* (from inverse iteration) 1075* 1076 CALL SGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, WR3, 1077 $ WI3, WORK, DUMMA( 1 ) ) 1078 IF( DUMMA( 1 ).LT.ULPINV ) 1079 $ RESULT( 13 ) = DUMMA( 1 )*ANINV 1080 END IF 1081* 1082* Call SORMHR for Left eigenvectors of A, do test 14 1083* 1084 NTEST = 14 1085 RESULT( 14 ) = ULPINV 1086* 1087 CALL SORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU, 1088 $ LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO ) 1089 IF( IINFO.NE.0 ) THEN 1090 WRITE( NOUNIT, FMT = 9999 )'SORMHR(L)', IINFO, N, JTYPE, 1091 $ IOLDSD 1092 INFO = ABS( IINFO ) 1093 IF( IINFO.LT.0 ) 1094 $ GO TO 250 1095 ELSE 1096* 1097* Test 14: | YA - WY | / ( |A| |Y| ulp ) 1098* 1099* (from inverse iteration) 1100* 1101 CALL SGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, WR3, 1102 $ WI3, WORK, DUMMA( 3 ) ) 1103 IF( DUMMA( 3 ).LT.ULPINV ) 1104 $ RESULT( 14 ) = DUMMA( 3 )*ANINV 1105 END IF 1106* 1107* End of Loop -- Check for RESULT(j) > THRESH 1108* 1109 250 CONTINUE 1110* 1111 NTESTT = NTESTT + NTEST 1112 CALL SLAFTS( 'SHS', N, N, JTYPE, NTEST, RESULT, IOLDSD, 1113 $ THRESH, NOUNIT, NERRS ) 1114* 1115 260 CONTINUE 1116 270 CONTINUE 1117* 1118* Summary 1119* 1120 CALL SLASUM( 'SHS', NOUNIT, NERRS, NTESTT ) 1121* 1122 RETURN 1123* 1124 9999 FORMAT( ' SCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 1125 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 1126 9998 FORMAT( ' SCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ', 1127 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X, 1128 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, 1129 $ ')' ) 1130 9997 FORMAT( ' SCHKHS: Selected ', A, ' Eigenvectors from ', A, 1131 $ ' do not match other eigenvectors ', 9X, 'N=', I6, 1132 $ ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 1133* 1134* End of SCHKHS 1135* 1136 END 1137