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