1*> \brief \b SGET24 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 SGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, 12* H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS, 13* LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, 14* RESULT, WORK, LWORK, IWORK, BWORK, INFO ) 15* 16* .. Scalar Arguments .. 17* LOGICAL COMP 18* INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT 19* REAL RCDEIN, RCDVIN, THRESH 20* .. 21* .. Array Arguments .. 22* LOGICAL BWORK( * ) 23* INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * ) 24* REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ), 25* $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ), 26* $ WI( * ), WIT( * ), WITMP( * ), WORK( * ), 27* $ WR( * ), WRT( * ), WRTMP( * ) 28* .. 29* 30* 31*> \par Purpose: 32* ============= 33*> 34*> \verbatim 35*> 36*> SGET24 checks the nonsymmetric eigenvalue (Schur form) problem 37*> expert driver SGEESX. 38*> 39*> If COMP = .FALSE., the first 13 of the following tests will be 40*> be performed on the input matrix A, and also tests 14 and 15 41*> if LWORK is sufficiently large. 42*> If COMP = .TRUE., all 17 test will be performed. 43*> 44*> (1) 0 if T is in Schur form, 1/ulp otherwise 45*> (no sorting of eigenvalues) 46*> 47*> (2) | A - VS T VS' | / ( n |A| ulp ) 48*> 49*> Here VS is the matrix of Schur eigenvectors, and T is in Schur 50*> form (no sorting of eigenvalues). 51*> 52*> (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues). 53*> 54*> (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T 55*> 1/ulp otherwise 56*> (no sorting of eigenvalues) 57*> 58*> (5) 0 if T(with VS) = T(without VS), 59*> 1/ulp otherwise 60*> (no sorting of eigenvalues) 61*> 62*> (6) 0 if eigenvalues(with VS) = eigenvalues(without VS), 63*> 1/ulp otherwise 64*> (no sorting of eigenvalues) 65*> 66*> (7) 0 if T is in Schur form, 1/ulp otherwise 67*> (with sorting of eigenvalues) 68*> 69*> (8) | A - VS T VS' | / ( n |A| ulp ) 70*> 71*> Here VS is the matrix of Schur eigenvectors, and T is in Schur 72*> form (with sorting of eigenvalues). 73*> 74*> (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues). 75*> 76*> (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T 77*> 1/ulp otherwise 78*> If workspace sufficient, also compare WR, WI with and 79*> without reciprocal condition numbers 80*> (with sorting of eigenvalues) 81*> 82*> (11) 0 if T(with VS) = T(without VS), 83*> 1/ulp otherwise 84*> If workspace sufficient, also compare T with and without 85*> reciprocal condition numbers 86*> (with sorting of eigenvalues) 87*> 88*> (12) 0 if eigenvalues(with VS) = eigenvalues(without VS), 89*> 1/ulp otherwise 90*> If workspace sufficient, also compare VS with and without 91*> reciprocal condition numbers 92*> (with sorting of eigenvalues) 93*> 94*> (13) if sorting worked and SDIM is the number of 95*> eigenvalues which were SELECTed 96*> If workspace sufficient, also compare SDIM with and 97*> without reciprocal condition numbers 98*> 99*> (14) if RCONDE the same no matter if VS and/or RCONDV computed 100*> 101*> (15) if RCONDV the same no matter if VS and/or RCONDE computed 102*> 103*> (16) |RCONDE - RCDEIN| / cond(RCONDE) 104*> 105*> RCONDE is the reciprocal average eigenvalue condition number 106*> computed by SGEESX and RCDEIN (the precomputed true value) 107*> is supplied as input. cond(RCONDE) is the condition number 108*> of RCONDE, and takes errors in computing RCONDE into account, 109*> so that the resulting quantity should be O(ULP). cond(RCONDE) 110*> is essentially given by norm(A)/RCONDV. 111*> 112*> (17) |RCONDV - RCDVIN| / cond(RCONDV) 113*> 114*> RCONDV is the reciprocal right invariant subspace condition 115*> number computed by SGEESX and RCDVIN (the precomputed true 116*> value) is supplied as input. cond(RCONDV) is the condition 117*> number of RCONDV, and takes errors in computing RCONDV into 118*> account, so that the resulting quantity should be O(ULP). 119*> cond(RCONDV) is essentially given by norm(A)/RCONDE. 120*> \endverbatim 121* 122* Arguments: 123* ========== 124* 125*> \param[in] COMP 126*> \verbatim 127*> COMP is LOGICAL 128*> COMP describes which input tests to perform: 129*> = .FALSE. if the computed condition numbers are not to 130*> be tested against RCDVIN and RCDEIN 131*> = .TRUE. if they are to be compared 132*> \endverbatim 133*> 134*> \param[in] JTYPE 135*> \verbatim 136*> JTYPE is INTEGER 137*> Type of input matrix. Used to label output if error occurs. 138*> \endverbatim 139*> 140*> \param[in] ISEED 141*> \verbatim 142*> ISEED is INTEGER array, dimension (4) 143*> If COMP = .FALSE., the random number generator seed 144*> used to produce matrix. 145*> If COMP = .TRUE., ISEED(1) = the number of the example. 146*> Used to label output if error occurs. 147*> \endverbatim 148*> 149*> \param[in] THRESH 150*> \verbatim 151*> THRESH is REAL 152*> A test will count as "failed" if the "error", computed as 153*> described above, exceeds THRESH. Note that the error 154*> is scaled to be O(1), so THRESH should be a reasonably 155*> small multiple of 1, e.g., 10 or 100. In particular, 156*> it should not depend on the precision (single vs. double) 157*> or the size of the matrix. It must be at least zero. 158*> \endverbatim 159*> 160*> \param[in] NOUNIT 161*> \verbatim 162*> NOUNIT is INTEGER 163*> The FORTRAN unit number for printing out error messages 164*> (e.g., if a routine returns INFO not equal to 0.) 165*> \endverbatim 166*> 167*> \param[in] N 168*> \verbatim 169*> N is INTEGER 170*> The dimension of A. N must be at least 0. 171*> \endverbatim 172*> 173*> \param[in,out] A 174*> \verbatim 175*> A is REAL array, dimension (LDA, N) 176*> Used to hold the matrix whose eigenvalues are to be 177*> computed. 178*> \endverbatim 179*> 180*> \param[in] LDA 181*> \verbatim 182*> LDA is INTEGER 183*> The leading dimension of A, and H. LDA must be at 184*> least 1 and at least N. 185*> \endverbatim 186*> 187*> \param[out] H 188*> \verbatim 189*> H is REAL array, dimension (LDA, N) 190*> Another copy of the test matrix A, modified by SGEESX. 191*> \endverbatim 192*> 193*> \param[out] HT 194*> \verbatim 195*> HT is REAL array, dimension (LDA, N) 196*> Yet another copy of the test matrix A, modified by SGEESX. 197*> \endverbatim 198*> 199*> \param[out] WR 200*> \verbatim 201*> WR is REAL array, dimension (N) 202*> \endverbatim 203*> 204*> \param[out] WI 205*> \verbatim 206*> WI is REAL array, dimension (N) 207*> 208*> The real and imaginary parts of the eigenvalues of A. 209*> On exit, WR + WI*i are the eigenvalues of the matrix in A. 210*> \endverbatim 211*> 212*> \param[out] WRT 213*> \verbatim 214*> WRT is REAL array, dimension (N) 215*> \endverbatim 216*> 217*> \param[out] WIT 218*> \verbatim 219*> WIT is REAL array, dimension (N) 220*> 221*> Like WR, WI, these arrays contain the eigenvalues of A, 222*> but those computed when SGEESX only computes a partial 223*> eigendecomposition, i.e. not Schur vectors 224*> \endverbatim 225*> 226*> \param[out] WRTMP 227*> \verbatim 228*> WRTMP is REAL array, dimension (N) 229*> \endverbatim 230*> 231*> \param[out] WITMP 232*> \verbatim 233*> WITMP is REAL array, dimension (N) 234*> 235*> Like WR, WI, these arrays contain the eigenvalues of A, 236*> but sorted by increasing real part. 237*> \endverbatim 238*> 239*> \param[out] VS 240*> \verbatim 241*> VS is REAL array, dimension (LDVS, N) 242*> VS holds the computed Schur vectors. 243*> \endverbatim 244*> 245*> \param[in] LDVS 246*> \verbatim 247*> LDVS is INTEGER 248*> Leading dimension of VS. Must be at least max(1, N). 249*> \endverbatim 250*> 251*> \param[out] VS1 252*> \verbatim 253*> VS1 is REAL array, dimension (LDVS, N) 254*> VS1 holds another copy of the computed Schur vectors. 255*> \endverbatim 256*> 257*> \param[in] RCDEIN 258*> \verbatim 259*> RCDEIN is REAL 260*> When COMP = .TRUE. RCDEIN holds the precomputed reciprocal 261*> condition number for the average of selected eigenvalues. 262*> \endverbatim 263*> 264*> \param[in] RCDVIN 265*> \verbatim 266*> RCDVIN is REAL 267*> When COMP = .TRUE. RCDVIN holds the precomputed reciprocal 268*> condition number for the selected right invariant subspace. 269*> \endverbatim 270*> 271*> \param[in] NSLCT 272*> \verbatim 273*> NSLCT is INTEGER 274*> When COMP = .TRUE. the number of selected eigenvalues 275*> corresponding to the precomputed values RCDEIN and RCDVIN. 276*> \endverbatim 277*> 278*> \param[in] ISLCT 279*> \verbatim 280*> ISLCT is INTEGER array, dimension (NSLCT) 281*> When COMP = .TRUE. ISLCT selects the eigenvalues of the 282*> input matrix corresponding to the precomputed values RCDEIN 283*> and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the 284*> eigenvalue with the J-th largest real part is selected. 285*> Not referenced if COMP = .FALSE. 286*> \endverbatim 287*> 288*> \param[out] RESULT 289*> \verbatim 290*> RESULT is REAL array, dimension (17) 291*> The values computed by the 17 tests described above. 292*> The values are currently limited to 1/ulp, to avoid 293*> overflow. 294*> \endverbatim 295*> 296*> \param[out] WORK 297*> \verbatim 298*> WORK is REAL array, dimension (LWORK) 299*> \endverbatim 300*> 301*> \param[in] LWORK 302*> \verbatim 303*> LWORK is INTEGER 304*> The number of entries in WORK to be passed to SGEESX. This 305*> must be at least 3*N, and N+N**2 if tests 14--16 are to 306*> be performed. 307*> \endverbatim 308*> 309*> \param[out] IWORK 310*> \verbatim 311*> IWORK is INTEGER array, dimension (N*N) 312*> \endverbatim 313*> 314*> \param[out] BWORK 315*> \verbatim 316*> BWORK is LOGICAL array, dimension (N) 317*> \endverbatim 318*> 319*> \param[out] INFO 320*> \verbatim 321*> INFO is INTEGER 322*> If 0, successful exit. 323*> If <0, input parameter -INFO had an incorrect value. 324*> If >0, SGEESX returned an error code, the absolute 325*> value of which is returned. 326*> \endverbatim 327* 328* Authors: 329* ======== 330* 331*> \author Univ. of Tennessee 332*> \author Univ. of California Berkeley 333*> \author Univ. of Colorado Denver 334*> \author NAG Ltd. 335* 336*> \date November 2011 337* 338*> \ingroup single_eig 339* 340* ===================================================================== 341 SUBROUTINE SGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, 342 $ H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS, 343 $ LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, 344 $ RESULT, WORK, LWORK, IWORK, BWORK, INFO ) 345* 346* -- LAPACK test routine (version 3.4.0) -- 347* -- LAPACK is a software package provided by Univ. of Tennessee, -- 348* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 349* November 2011 350* 351* .. Scalar Arguments .. 352 LOGICAL COMP 353 INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT 354 REAL RCDEIN, RCDVIN, THRESH 355* .. 356* .. Array Arguments .. 357 LOGICAL BWORK( * ) 358 INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * ) 359 REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ), 360 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ), 361 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ), 362 $ WR( * ), WRT( * ), WRTMP( * ) 363* .. 364* 365* ===================================================================== 366* 367* .. Parameters .. 368 REAL ZERO, ONE 369 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 370 REAL EPSIN 371 PARAMETER ( EPSIN = 5.9605E-8 ) 372* .. 373* .. Local Scalars .. 374 CHARACTER SORT 375 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK, 376 $ RSUB, SDIM, SDIM1 377 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV, 378 $ SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN, 379 $ VRMIN, WNORM 380* .. 381* .. Local Arrays .. 382 INTEGER IPNT( 20 ) 383* .. 384* .. Arrays in Common .. 385 LOGICAL SELVAL( 20 ) 386 REAL SELWI( 20 ), SELWR( 20 ) 387* .. 388* .. Scalars in Common .. 389 INTEGER SELDIM, SELOPT 390* .. 391* .. Common blocks .. 392 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI 393* .. 394* .. External Functions .. 395 LOGICAL SSLECT 396 REAL SLAMCH, SLANGE 397 EXTERNAL SSLECT, SLAMCH, SLANGE 398* .. 399* .. External Subroutines .. 400 EXTERNAL SCOPY, SGEESX, SGEMM, SLACPY, SORT01, XERBLA 401* .. 402* .. Intrinsic Functions .. 403 INTRINSIC ABS, MAX, MIN, REAL, SIGN, SQRT 404* .. 405* .. Executable Statements .. 406* 407* Check for errors 408* 409 INFO = 0 410 IF( THRESH.LT.ZERO ) THEN 411 INFO = -3 412 ELSE IF( NOUNIT.LE.0 ) THEN 413 INFO = -5 414 ELSE IF( N.LT.0 ) THEN 415 INFO = -6 416 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN 417 INFO = -8 418 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN 419 INFO = -18 420 ELSE IF( LWORK.LT.3*N ) THEN 421 INFO = -26 422 END IF 423* 424 IF( INFO.NE.0 ) THEN 425 CALL XERBLA( 'SGET24', -INFO ) 426 RETURN 427 END IF 428* 429* Quick return if nothing to do 430* 431 DO 10 I = 1, 17 432 RESULT( I ) = -ONE 433 10 CONTINUE 434* 435 IF( N.EQ.0 ) 436 $ RETURN 437* 438* Important constants 439* 440 SMLNUM = SLAMCH( 'Safe minimum' ) 441 ULP = SLAMCH( 'Precision' ) 442 ULPINV = ONE / ULP 443* 444* Perform tests (1)-(13) 445* 446 SELOPT = 0 447 LIWORK = N*N 448 DO 120 ISORT = 0, 1 449 IF( ISORT.EQ.0 ) THEN 450 SORT = 'N' 451 RSUB = 0 452 ELSE 453 SORT = 'S' 454 RSUB = 6 455 END IF 456* 457* Compute Schur form and Schur vectors, and test them 458* 459 CALL SLACPY( 'F', N, N, A, LDA, H, LDA ) 460 CALL SGEESX( 'V', SORT, SSLECT, 'N', N, H, LDA, SDIM, WR, WI, 461 $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, 462 $ LIWORK, BWORK, IINFO ) 463 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 464 RESULT( 1+RSUB ) = ULPINV 465 IF( JTYPE.NE.22 ) THEN 466 WRITE( NOUNIT, FMT = 9998 )'SGEESX1', IINFO, N, JTYPE, 467 $ ISEED 468 ELSE 469 WRITE( NOUNIT, FMT = 9999 )'SGEESX1', IINFO, N, 470 $ ISEED( 1 ) 471 END IF 472 INFO = ABS( IINFO ) 473 RETURN 474 END IF 475 IF( ISORT.EQ.0 ) THEN 476 CALL SCOPY( N, WR, 1, WRTMP, 1 ) 477 CALL SCOPY( N, WI, 1, WITMP, 1 ) 478 END IF 479* 480* Do Test (1) or Test (7) 481* 482 RESULT( 1+RSUB ) = ZERO 483 DO 30 J = 1, N - 2 484 DO 20 I = J + 2, N 485 IF( H( I, J ).NE.ZERO ) 486 $ RESULT( 1+RSUB ) = ULPINV 487 20 CONTINUE 488 30 CONTINUE 489 DO 40 I = 1, N - 2 490 IF( H( I+1, I ).NE.ZERO .AND. H( I+2, I+1 ).NE.ZERO ) 491 $ RESULT( 1+RSUB ) = ULPINV 492 40 CONTINUE 493 DO 50 I = 1, N - 1 494 IF( H( I+1, I ).NE.ZERO ) THEN 495 IF( H( I, I ).NE.H( I+1, I+1 ) .OR. H( I, I+1 ).EQ. 496 $ ZERO .OR. SIGN( ONE, H( I+1, I ) ).EQ. 497 $ SIGN( ONE, H( I, I+1 ) ) )RESULT( 1+RSUB ) = ULPINV 498 END IF 499 50 CONTINUE 500* 501* Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP) 502* 503* Copy A to VS1, used as workspace 504* 505 CALL SLACPY( ' ', N, N, A, LDA, VS1, LDVS ) 506* 507* Compute Q*H and store in HT. 508* 509 CALL SGEMM( 'No transpose', 'No transpose', N, N, N, ONE, VS, 510 $ LDVS, H, LDA, ZERO, HT, LDA ) 511* 512* Compute A - Q*H*Q' 513* 514 CALL SGEMM( 'No transpose', 'Transpose', N, N, N, -ONE, HT, 515 $ LDA, VS, LDVS, ONE, VS1, LDVS ) 516* 517 ANORM = MAX( SLANGE( '1', N, N, A, LDA, WORK ), SMLNUM ) 518 WNORM = SLANGE( '1', N, N, VS1, LDVS, WORK ) 519* 520 IF( ANORM.GT.WNORM ) THEN 521 RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP ) 522 ELSE 523 IF( ANORM.LT.ONE ) THEN 524 RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / 525 $ ( N*ULP ) 526 ELSE 527 RESULT( 2+RSUB ) = MIN( WNORM / ANORM, REAL( N ) ) / 528 $ ( N*ULP ) 529 END IF 530 END IF 531* 532* Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP ) 533* 534 CALL SORT01( 'Columns', N, N, VS, LDVS, WORK, LWORK, 535 $ RESULT( 3+RSUB ) ) 536* 537* Do Test (4) or Test (10) 538* 539 RESULT( 4+RSUB ) = ZERO 540 DO 60 I = 1, N 541 IF( H( I, I ).NE.WR( I ) ) 542 $ RESULT( 4+RSUB ) = ULPINV 543 60 CONTINUE 544 IF( N.GT.1 ) THEN 545 IF( H( 2, 1 ).EQ.ZERO .AND. WI( 1 ).NE.ZERO ) 546 $ RESULT( 4+RSUB ) = ULPINV 547 IF( H( N, N-1 ).EQ.ZERO .AND. WI( N ).NE.ZERO ) 548 $ RESULT( 4+RSUB ) = ULPINV 549 END IF 550 DO 70 I = 1, N - 1 551 IF( H( I+1, I ).NE.ZERO ) THEN 552 TMP = SQRT( ABS( H( I+1, I ) ) )* 553 $ SQRT( ABS( H( I, I+1 ) ) ) 554 RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ), 555 $ ABS( WI( I )-TMP ) / 556 $ MAX( ULP*TMP, SMLNUM ) ) 557 RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ), 558 $ ABS( WI( I+1 )+TMP ) / 559 $ MAX( ULP*TMP, SMLNUM ) ) 560 ELSE IF( I.GT.1 ) THEN 561 IF( H( I+1, I ).EQ.ZERO .AND. H( I, I-1 ).EQ.ZERO .AND. 562 $ WI( I ).NE.ZERO )RESULT( 4+RSUB ) = ULPINV 563 END IF 564 70 CONTINUE 565* 566* Do Test (5) or Test (11) 567* 568 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 569 CALL SGEESX( 'N', SORT, SSLECT, 'N', N, HT, LDA, SDIM, WRT, 570 $ WIT, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, 571 $ LIWORK, BWORK, IINFO ) 572 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 573 RESULT( 5+RSUB ) = ULPINV 574 IF( JTYPE.NE.22 ) THEN 575 WRITE( NOUNIT, FMT = 9998 )'SGEESX2', IINFO, N, JTYPE, 576 $ ISEED 577 ELSE 578 WRITE( NOUNIT, FMT = 9999 )'SGEESX2', IINFO, N, 579 $ ISEED( 1 ) 580 END IF 581 INFO = ABS( IINFO ) 582 GO TO 250 583 END IF 584* 585 RESULT( 5+RSUB ) = ZERO 586 DO 90 J = 1, N 587 DO 80 I = 1, N 588 IF( H( I, J ).NE.HT( I, J ) ) 589 $ RESULT( 5+RSUB ) = ULPINV 590 80 CONTINUE 591 90 CONTINUE 592* 593* Do Test (6) or Test (12) 594* 595 RESULT( 6+RSUB ) = ZERO 596 DO 100 I = 1, N 597 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) ) 598 $ RESULT( 6+RSUB ) = ULPINV 599 100 CONTINUE 600* 601* Do Test (13) 602* 603 IF( ISORT.EQ.1 ) THEN 604 RESULT( 13 ) = ZERO 605 KNTEIG = 0 606 DO 110 I = 1, N 607 IF( SSLECT( WR( I ), WI( I ) ) .OR. 608 $ SSLECT( WR( I ), -WI( I ) ) )KNTEIG = KNTEIG + 1 609 IF( I.LT.N ) THEN 610 IF( ( SSLECT( WR( I+1 ), WI( I+1 ) ) .OR. 611 $ SSLECT( WR( I+1 ), -WI( I+1 ) ) ) .AND. 612 $ ( .NOT.( SSLECT( WR( I ), 613 $ WI( I ) ) .OR. SSLECT( WR( I ), 614 $ -WI( I ) ) ) ) .AND. IINFO.NE.N+2 )RESULT( 13 ) 615 $ = ULPINV 616 END IF 617 110 CONTINUE 618 IF( SDIM.NE.KNTEIG ) 619 $ RESULT( 13 ) = ULPINV 620 END IF 621* 622 120 CONTINUE 623* 624* If there is enough workspace, perform tests (14) and (15) 625* as well as (10) through (13) 626* 627 IF( LWORK.GE.N+( N*N ) / 2 ) THEN 628* 629* Compute both RCONDE and RCONDV with VS 630* 631 SORT = 'S' 632 RESULT( 14 ) = ZERO 633 RESULT( 15 ) = ZERO 634 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 635 CALL SGEESX( 'V', SORT, SSLECT, 'B', N, HT, LDA, SDIM1, WRT, 636 $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK, 637 $ IWORK, LIWORK, BWORK, IINFO ) 638 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 639 RESULT( 14 ) = ULPINV 640 RESULT( 15 ) = ULPINV 641 IF( JTYPE.NE.22 ) THEN 642 WRITE( NOUNIT, FMT = 9998 )'SGEESX3', IINFO, N, JTYPE, 643 $ ISEED 644 ELSE 645 WRITE( NOUNIT, FMT = 9999 )'SGEESX3', IINFO, N, 646 $ ISEED( 1 ) 647 END IF 648 INFO = ABS( IINFO ) 649 GO TO 250 650 END IF 651* 652* Perform tests (10), (11), (12), and (13) 653* 654 DO 140 I = 1, N 655 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) ) 656 $ RESULT( 10 ) = ULPINV 657 DO 130 J = 1, N 658 IF( H( I, J ).NE.HT( I, J ) ) 659 $ RESULT( 11 ) = ULPINV 660 IF( VS( I, J ).NE.VS1( I, J ) ) 661 $ RESULT( 12 ) = ULPINV 662 130 CONTINUE 663 140 CONTINUE 664 IF( SDIM.NE.SDIM1 ) 665 $ RESULT( 13 ) = ULPINV 666* 667* Compute both RCONDE and RCONDV without VS, and compare 668* 669 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 670 CALL SGEESX( 'N', SORT, SSLECT, 'B', N, HT, LDA, SDIM1, WRT, 671 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, 672 $ IWORK, LIWORK, BWORK, IINFO ) 673 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 674 RESULT( 14 ) = ULPINV 675 RESULT( 15 ) = ULPINV 676 IF( JTYPE.NE.22 ) THEN 677 WRITE( NOUNIT, FMT = 9998 )'SGEESX4', IINFO, N, JTYPE, 678 $ ISEED 679 ELSE 680 WRITE( NOUNIT, FMT = 9999 )'SGEESX4', IINFO, N, 681 $ ISEED( 1 ) 682 END IF 683 INFO = ABS( IINFO ) 684 GO TO 250 685 END IF 686* 687* Perform tests (14) and (15) 688* 689 IF( RCNDE1.NE.RCONDE ) 690 $ RESULT( 14 ) = ULPINV 691 IF( RCNDV1.NE.RCONDV ) 692 $ RESULT( 15 ) = ULPINV 693* 694* Perform tests (10), (11), (12), and (13) 695* 696 DO 160 I = 1, N 697 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) ) 698 $ RESULT( 10 ) = ULPINV 699 DO 150 J = 1, N 700 IF( H( I, J ).NE.HT( I, J ) ) 701 $ RESULT( 11 ) = ULPINV 702 IF( VS( I, J ).NE.VS1( I, J ) ) 703 $ RESULT( 12 ) = ULPINV 704 150 CONTINUE 705 160 CONTINUE 706 IF( SDIM.NE.SDIM1 ) 707 $ RESULT( 13 ) = ULPINV 708* 709* Compute RCONDE with VS, and compare 710* 711 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 712 CALL SGEESX( 'V', SORT, SSLECT, 'E', N, HT, LDA, SDIM1, WRT, 713 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, 714 $ IWORK, LIWORK, BWORK, IINFO ) 715 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 716 RESULT( 14 ) = ULPINV 717 IF( JTYPE.NE.22 ) THEN 718 WRITE( NOUNIT, FMT = 9998 )'SGEESX5', IINFO, N, JTYPE, 719 $ ISEED 720 ELSE 721 WRITE( NOUNIT, FMT = 9999 )'SGEESX5', IINFO, N, 722 $ ISEED( 1 ) 723 END IF 724 INFO = ABS( IINFO ) 725 GO TO 250 726 END IF 727* 728* Perform test (14) 729* 730 IF( RCNDE1.NE.RCONDE ) 731 $ RESULT( 14 ) = ULPINV 732* 733* Perform tests (10), (11), (12), and (13) 734* 735 DO 180 I = 1, N 736 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) ) 737 $ RESULT( 10 ) = ULPINV 738 DO 170 J = 1, N 739 IF( H( I, J ).NE.HT( I, J ) ) 740 $ RESULT( 11 ) = ULPINV 741 IF( VS( I, J ).NE.VS1( I, J ) ) 742 $ RESULT( 12 ) = ULPINV 743 170 CONTINUE 744 180 CONTINUE 745 IF( SDIM.NE.SDIM1 ) 746 $ RESULT( 13 ) = ULPINV 747* 748* Compute RCONDE without VS, and compare 749* 750 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 751 CALL SGEESX( 'N', SORT, SSLECT, 'E', N, HT, LDA, SDIM1, WRT, 752 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, 753 $ IWORK, LIWORK, BWORK, IINFO ) 754 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 755 RESULT( 14 ) = ULPINV 756 IF( JTYPE.NE.22 ) THEN 757 WRITE( NOUNIT, FMT = 9998 )'SGEESX6', IINFO, N, JTYPE, 758 $ ISEED 759 ELSE 760 WRITE( NOUNIT, FMT = 9999 )'SGEESX6', IINFO, N, 761 $ ISEED( 1 ) 762 END IF 763 INFO = ABS( IINFO ) 764 GO TO 250 765 END IF 766* 767* Perform test (14) 768* 769 IF( RCNDE1.NE.RCONDE ) 770 $ RESULT( 14 ) = ULPINV 771* 772* Perform tests (10), (11), (12), and (13) 773* 774 DO 200 I = 1, N 775 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) ) 776 $ RESULT( 10 ) = ULPINV 777 DO 190 J = 1, N 778 IF( H( I, J ).NE.HT( I, J ) ) 779 $ RESULT( 11 ) = ULPINV 780 IF( VS( I, J ).NE.VS1( I, J ) ) 781 $ RESULT( 12 ) = ULPINV 782 190 CONTINUE 783 200 CONTINUE 784 IF( SDIM.NE.SDIM1 ) 785 $ RESULT( 13 ) = ULPINV 786* 787* Compute RCONDV with VS, and compare 788* 789 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 790 CALL SGEESX( 'V', SORT, SSLECT, 'V', N, HT, LDA, SDIM1, WRT, 791 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, 792 $ IWORK, LIWORK, BWORK, IINFO ) 793 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 794 RESULT( 15 ) = ULPINV 795 IF( JTYPE.NE.22 ) THEN 796 WRITE( NOUNIT, FMT = 9998 )'SGEESX7', IINFO, N, JTYPE, 797 $ ISEED 798 ELSE 799 WRITE( NOUNIT, FMT = 9999 )'SGEESX7', IINFO, N, 800 $ ISEED( 1 ) 801 END IF 802 INFO = ABS( IINFO ) 803 GO TO 250 804 END IF 805* 806* Perform test (15) 807* 808 IF( RCNDV1.NE.RCONDV ) 809 $ RESULT( 15 ) = ULPINV 810* 811* Perform tests (10), (11), (12), and (13) 812* 813 DO 220 I = 1, N 814 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) ) 815 $ RESULT( 10 ) = ULPINV 816 DO 210 J = 1, N 817 IF( H( I, J ).NE.HT( I, J ) ) 818 $ RESULT( 11 ) = ULPINV 819 IF( VS( I, J ).NE.VS1( I, J ) ) 820 $ RESULT( 12 ) = ULPINV 821 210 CONTINUE 822 220 CONTINUE 823 IF( SDIM.NE.SDIM1 ) 824 $ RESULT( 13 ) = ULPINV 825* 826* Compute RCONDV without VS, and compare 827* 828 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 829 CALL SGEESX( 'N', SORT, SSLECT, 'V', N, HT, LDA, SDIM1, WRT, 830 $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, 831 $ IWORK, LIWORK, BWORK, IINFO ) 832 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 833 RESULT( 15 ) = ULPINV 834 IF( JTYPE.NE.22 ) THEN 835 WRITE( NOUNIT, FMT = 9998 )'SGEESX8', IINFO, N, JTYPE, 836 $ ISEED 837 ELSE 838 WRITE( NOUNIT, FMT = 9999 )'SGEESX8', IINFO, N, 839 $ ISEED( 1 ) 840 END IF 841 INFO = ABS( IINFO ) 842 GO TO 250 843 END IF 844* 845* Perform test (15) 846* 847 IF( RCNDV1.NE.RCONDV ) 848 $ RESULT( 15 ) = ULPINV 849* 850* Perform tests (10), (11), (12), and (13) 851* 852 DO 240 I = 1, N 853 IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) ) 854 $ RESULT( 10 ) = ULPINV 855 DO 230 J = 1, N 856 IF( H( I, J ).NE.HT( I, J ) ) 857 $ RESULT( 11 ) = ULPINV 858 IF( VS( I, J ).NE.VS1( I, J ) ) 859 $ RESULT( 12 ) = ULPINV 860 230 CONTINUE 861 240 CONTINUE 862 IF( SDIM.NE.SDIM1 ) 863 $ RESULT( 13 ) = ULPINV 864* 865 END IF 866* 867 250 CONTINUE 868* 869* If there are precomputed reciprocal condition numbers, compare 870* computed values with them. 871* 872 IF( COMP ) THEN 873* 874* First set up SELOPT, SELDIM, SELVAL, SELWR, and SELWI so that 875* the logical function SSLECT selects the eigenvalues specified 876* by NSLCT and ISLCT. 877* 878 SELDIM = N 879 SELOPT = 1 880 EPS = MAX( ULP, EPSIN ) 881 DO 260 I = 1, N 882 IPNT( I ) = I 883 SELVAL( I ) = .FALSE. 884 SELWR( I ) = WRTMP( I ) 885 SELWI( I ) = WITMP( I ) 886 260 CONTINUE 887 DO 280 I = 1, N - 1 888 KMIN = I 889 VRMIN = WRTMP( I ) 890 VIMIN = WITMP( I ) 891 DO 270 J = I + 1, N 892 IF( WRTMP( J ).LT.VRMIN ) THEN 893 KMIN = J 894 VRMIN = WRTMP( J ) 895 VIMIN = WITMP( J ) 896 END IF 897 270 CONTINUE 898 WRTMP( KMIN ) = WRTMP( I ) 899 WITMP( KMIN ) = WITMP( I ) 900 WRTMP( I ) = VRMIN 901 WITMP( I ) = VIMIN 902 ITMP = IPNT( I ) 903 IPNT( I ) = IPNT( KMIN ) 904 IPNT( KMIN ) = ITMP 905 280 CONTINUE 906 DO 290 I = 1, NSLCT 907 SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE. 908 290 CONTINUE 909* 910* Compute condition numbers 911* 912 CALL SLACPY( 'F', N, N, A, LDA, HT, LDA ) 913 CALL SGEESX( 'N', 'S', SSLECT, 'B', N, HT, LDA, SDIM1, WRT, 914 $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK, 915 $ IWORK, LIWORK, BWORK, IINFO ) 916 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN 917 RESULT( 16 ) = ULPINV 918 RESULT( 17 ) = ULPINV 919 WRITE( NOUNIT, FMT = 9999 )'SGEESX9', IINFO, N, ISEED( 1 ) 920 INFO = ABS( IINFO ) 921 GO TO 300 922 END IF 923* 924* Compare condition number for average of selected eigenvalues 925* taking its condition number into account 926* 927 ANORM = SLANGE( '1', N, N, A, LDA, WORK ) 928 V = MAX( REAL( N )*EPS*ANORM, SMLNUM ) 929 IF( ANORM.EQ.ZERO ) 930 $ V = ONE 931 IF( V.GT.RCONDV ) THEN 932 TOL = ONE 933 ELSE 934 TOL = V / RCONDV 935 END IF 936 IF( V.GT.RCDVIN ) THEN 937 TOLIN = ONE 938 ELSE 939 TOLIN = V / RCDVIN 940 END IF 941 TOL = MAX( TOL, SMLNUM / EPS ) 942 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 943 IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN 944 RESULT( 16 ) = ULPINV 945 ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN 946 RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL ) 947 ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN 948 RESULT( 16 ) = ULPINV 949 ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN 950 RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN ) 951 ELSE 952 RESULT( 16 ) = ONE 953 END IF 954* 955* Compare condition numbers for right invariant subspace 956* taking its condition number into account 957* 958 IF( V.GT.RCONDV*RCONDE ) THEN 959 TOL = RCONDV 960 ELSE 961 TOL = V / RCONDE 962 END IF 963 IF( V.GT.RCDVIN*RCDEIN ) THEN 964 TOLIN = RCDVIN 965 ELSE 966 TOLIN = V / RCDEIN 967 END IF 968 TOL = MAX( TOL, SMLNUM / EPS ) 969 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 970 IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN 971 RESULT( 17 ) = ULPINV 972 ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN 973 RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL ) 974 ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN 975 RESULT( 17 ) = ULPINV 976 ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN 977 RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN ) 978 ELSE 979 RESULT( 17 ) = ONE 980 END IF 981* 982 300 CONTINUE 983* 984 END IF 985* 986 9999 FORMAT( ' SGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 987 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 ) 988 9998 FORMAT( ' SGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 989 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 990* 991 RETURN 992* 993* End of SGET24 994* 995 END 996