LAPACK: TESTING/EIG/cget24.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
00002      $                   H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN,
00003      $                   RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK,
00004      $                   LWORK, RWORK, BWORK, INFO )
00005 *
00006 *  -- LAPACK test routine (version 3.1) --
00007 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       LOGICAL            COMP
00012       INTEGER            INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
00013      $                   NSLCT
00014       REAL               RCDEIN, RCDVIN, THRESH
00015 *     ..
00016 *     .. Array Arguments ..
00017       LOGICAL            BWORK( * )
00018       INTEGER            ISEED( 4 ), ISLCT( * )
00019       REAL               RESULT( 17 ), RWORK( * )
00020       COMPLEX            A( LDA, * ), H( LDA, * ), HT( LDA, * ),
00021      $                   VS( LDVS, * ), VS1( LDVS, * ), W( * ),
00022      $                   WORK( * ), WT( * ), WTMP( * )
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *     CGET24 checks the nonsymmetric eigenvalue (Schur form) problem
00029 *     expert driver CGEESX.
00030 *
00031 *     If COMP = .FALSE., the first 13 of the following tests will be
00032 *     be performed on the input matrix A, and also tests 14 and 15
00033 *     if LWORK is sufficiently large.
00034 *     If COMP = .TRUE., all 17 test will be performed.
00035 *
00036 *     (1)     0 if T is in Schur form, 1/ulp otherwise
00037 *            (no sorting of eigenvalues)
00038 *
00039 *     (2)     | A - VS T VS' | / ( n |A| ulp )
00040 *
00041 *       Here VS is the matrix of Schur eigenvectors, and T is in Schur
00042 *       form  (no sorting of eigenvalues).
00043 *
00044 *     (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
00045 *
00046 *     (4)     0     if W are eigenvalues of T
00047 *             1/ulp otherwise
00048 *             (no sorting of eigenvalues)
00049 *
00050 *     (5)     0     if T(with VS) = T(without VS),
00051 *             1/ulp otherwise
00052 *             (no sorting of eigenvalues)
00053 *
00054 *     (6)     0     if eigenvalues(with VS) = eigenvalues(without VS),
00055 *             1/ulp otherwise
00056 *             (no sorting of eigenvalues)
00057 *
00058 *     (7)     0 if T is in Schur form, 1/ulp otherwise
00059 *             (with sorting of eigenvalues)
00060 *
00061 *     (8)     | A - VS T VS' | / ( n |A| ulp )
00062 *
00063 *       Here VS is the matrix of Schur eigenvectors, and T is in Schur
00064 *       form  (with sorting of eigenvalues).
00065 *
00066 *     (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
00067 *
00068 *     (10)    0     if W are eigenvalues of T
00069 *             1/ulp otherwise
00070 *             If workspace sufficient, also compare W with and
00071 *             without reciprocal condition numbers
00072 *             (with sorting of eigenvalues)
00073 *
00074 *     (11)    0     if T(with VS) = T(without VS),
00075 *             1/ulp otherwise
00076 *             If workspace sufficient, also compare T with and without
00077 *             reciprocal condition numbers
00078 *             (with sorting of eigenvalues)
00079 *
00080 *     (12)    0     if eigenvalues(with VS) = eigenvalues(without VS),
00081 *             1/ulp otherwise
00082 *             If workspace sufficient, also compare VS with and without
00083 *             reciprocal condition numbers
00084 *             (with sorting of eigenvalues)
00085 *
00086 *     (13)    if sorting worked and SDIM is the number of
00087 *             eigenvalues which were SELECTed
00088 *             If workspace sufficient, also compare SDIM with and
00089 *             without reciprocal condition numbers
00090 *
00091 *     (14)    if RCONDE the same no matter if VS and/or RCONDV computed
00092 *
00093 *     (15)    if RCONDV the same no matter if VS and/or RCONDE computed
00094 *
00095 *     (16)  |RCONDE - RCDEIN| / cond(RCONDE)
00096 *
00097 *        RCONDE is the reciprocal average eigenvalue condition number
00098 *        computed by CGEESX and RCDEIN (the precomputed true value)
00099 *        is supplied as input.  cond(RCONDE) is the condition number
00100 *        of RCONDE, and takes errors in computing RCONDE into account,
00101 *        so that the resulting quantity should be O(ULP). cond(RCONDE)
00102 *        is essentially given by norm(A)/RCONDV.
00103 *
00104 *     (17)  |RCONDV - RCDVIN| / cond(RCONDV)
00105 *
00106 *        RCONDV is the reciprocal right invariant subspace condition
00107 *        number computed by CGEESX and RCDVIN (the precomputed true
00108 *        value) is supplied as input. cond(RCONDV) is the condition
00109 *        number of RCONDV, and takes errors in computing RCONDV into
00110 *        account, so that the resulting quantity should be O(ULP).
00111 *        cond(RCONDV) is essentially given by norm(A)/RCONDE.
00112 *
00113 *  Arguments
00114 *  =========
00115 *
00116 *  COMP    (input) LOGICAL
00117 *          COMP describes which input tests to perform:
00118 *            = .FALSE. if the computed condition numbers are not to
00119 *                      be tested against RCDVIN and RCDEIN
00120 *            = .TRUE.  if they are to be compared
00121 *
00122 *  JTYPE   (input) INTEGER
00123 *          Type of input matrix. Used to label output if error occurs.
00124 *
00125 *  ISEED   (input) INTEGER array, dimension (4)
00126 *          If COMP = .FALSE., the random number generator seed
00127 *          used to produce matrix.
00128 *          If COMP = .TRUE., ISEED(1) = the number of the example.
00129 *          Used to label output if error occurs.
00130 *
00131 *  THRESH  (input) REAL
00132 *          A test will count as "failed" if the "error", computed as
00133 *          described above, exceeds THRESH.  Note that the error
00134 *          is scaled to be O(1), so THRESH should be a reasonably
00135 *          small multiple of 1, e.g., 10 or 100.  In particular,
00136 *          it should not depend on the precision (single vs. double)
00137 *          or the size of the matrix.  It must be at least zero.
00138 *
00139 *  NOUNIT  (input) INTEGER
00140 *          The FORTRAN unit number for printing out error messages
00141 *          (e.g., if a routine returns INFO not equal to 0.)
00142 *
00143 *  N       (input) INTEGER
00144 *          The dimension of A. N must be at least 0.
00145 *
00146 *  A       (input/output) COMPLEX array, dimension (LDA, N)
00147 *          Used to hold the matrix whose eigenvalues are to be
00148 *          computed.
00149 *
00150 *  LDA     (input) INTEGER
00151 *          The leading dimension of A, and H. LDA must be at
00152 *          least 1 and at least N.
00153 *
00154 *  H       (workspace) COMPLEX array, dimension (LDA, N)
00155 *          Another copy of the test matrix A, modified by CGEESX.
00156 *
00157 *  HT      (workspace) COMPLEX array, dimension (LDA, N)
00158 *          Yet another copy of the test matrix A, modified by CGEESX.
00159 *
00160 *  W       (workspace) COMPLEX array, dimension (N)
00161 *          The computed eigenvalues of A.
00162 *
00163 *  WT      (workspace) COMPLEX array, dimension (N)
00164 *          Like W, this array contains the eigenvalues of A,
00165 *          but those computed when CGEESX only computes a partial
00166 *          eigendecomposition, i.e. not Schur vectors
00167 *
00168 *  WTMP    (workspace) COMPLEX array, dimension (N)
00169 *          Like W, this array contains the eigenvalues of A,
00170 *          but sorted by increasing real or imaginary part.
00171 *
00172 *  VS      (workspace) COMPLEX array, dimension (LDVS, N)
00173 *          VS holds the computed Schur vectors.
00174 *
00175 *  LDVS    (input) INTEGER
00176 *          Leading dimension of VS. Must be at least max(1, N).
00177 *
00178 *  VS1     (workspace) COMPLEX array, dimension (LDVS, N)
00179 *          VS1 holds another copy of the computed Schur vectors.
00180 *
00181 *  RCDEIN  (input) REAL
00182 *          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
00183 *          condition number for the average of selected eigenvalues.
00184 *
00185 *  RCDVIN  (input) REAL
00186 *          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
00187 *          condition number for the selected right invariant subspace.
00188 *
00189 *  NSLCT   (input) INTEGER
00190 *          When COMP = .TRUE. the number of selected eigenvalues
00191 *          corresponding to the precomputed values RCDEIN and RCDVIN.
00192 *
00193 *  ISLCT   (input) INTEGER array, dimension (NSLCT)
00194 *          When COMP = .TRUE. ISLCT selects the eigenvalues of the
00195 *          input matrix corresponding to the precomputed values RCDEIN
00196 *          and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
00197 *          eigenvalue with the J-th largest real or imaginary part is
00198 *          selected. The real part is used if ISRT = 0, and the
00199 *          imaginary part if ISRT = 1.
00200 *          Not referenced if COMP = .FALSE.
00201 *
00202 *  ISRT    (input) INTEGER
00203 *          When COMP = .TRUE., ISRT describes how ISLCT is used to
00204 *          choose a subset of the spectrum.
00205 *          Not referenced if COMP = .FALSE.
00206 *
00207 *  RESULT  (output) REAL array, dimension (17)
00208 *          The values computed by the 17 tests described above.
00209 *          The values are currently limited to 1/ulp, to avoid
00210 *          overflow.
00211 *
00212 *  WORK    (workspace) COMPLEX array, dimension (2*N*N)
00213 *
00214 *  LWORK   (input) INTEGER
00215 *          The number of entries in WORK to be passed to CGEESX. This
00216 *          must be at least 2*N, and N*(N+1)/2 if tests 14--16 are to
00217 *          be performed.
00218 *
00219 *  RWORK   (workspace) REAL array, dimension (N)
00220 *
00221 *  BWORK   (workspace) LOGICAL array, dimension (N)
00222 *
00223 *  INFO    (output) INTEGER
00224 *          If 0,  successful exit.
00225 *          If <0, input parameter -INFO had an incorrect value.
00226 *          If >0, CGEESX returned an error code, the absolute
00227 *                 value of which is returned.
00228 *
00229 *  =====================================================================
00230 *
00231 *     .. Parameters ..
00232       COMPLEX            CZERO, CONE
00233       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00234      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00235       REAL               ZERO, ONE
00236       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00237       REAL               EPSIN
00238       PARAMETER          ( EPSIN = 5.9605E-8 )
00239 *     ..
00240 *     .. Local Scalars ..
00241       CHARACTER          SORT
00242       INTEGER            I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB,
00243      $                   SDIM, SDIM1
00244       REAL               ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
00245      $                   SMLNUM, TOL, TOLIN, ULP, ULPINV, V, VRICMP,
00246      $                   VRIMIN, WNORM
00247       COMPLEX            CTMP
00248 *     ..
00249 *     .. Local Arrays ..
00250       INTEGER            IPNT( 20 )
00251 *     ..
00252 *     .. External Functions ..
00253       LOGICAL            CSLECT
00254       REAL               CLANGE, SLAMCH
00255       EXTERNAL           CSLECT, CLANGE, SLAMCH
00256 *     ..
00257 *     .. External Subroutines ..
00258       EXTERNAL           CCOPY, CGEESX, CGEMM, CLACPY, CUNT01, XERBLA
00259 *     ..
00260 *     .. Intrinsic Functions ..
00261       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
00262 *     ..
00263 *     .. Arrays in Common ..
00264       LOGICAL            SELVAL( 20 )
00265       REAL               SELWI( 20 ), SELWR( 20 )
00266 *     ..
00267 *     .. Scalars in Common ..
00268       INTEGER            SELDIM, SELOPT
00269 *     ..
00270 *     .. Common blocks ..
00271       COMMON             / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
00272 *     ..
00273 *     .. Executable Statements ..
00274 *
00275 *     Check for errors
00276 *
00277       INFO = 0
00278       IF( THRESH.LT.ZERO ) THEN
00279          INFO = -3
00280       ELSE IF( NOUNIT.LE.0 ) THEN
00281          INFO = -5
00282       ELSE IF( N.LT.0 ) THEN
00283          INFO = -6
00284       ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
00285          INFO = -8
00286       ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN
00287          INFO = -15
00288       ELSE IF( LWORK.LT.2*N ) THEN
00289          INFO = -24
00290       END IF
00291 *
00292       IF( INFO.NE.0 ) THEN
00293          CALL XERBLA( 'CGET24', -INFO )
00294          RETURN
00295       END IF
00296 *
00297 *     Quick return if nothing to do
00298 *
00299       DO 10 I = 1, 17
00300          RESULT( I ) = -ONE
00301    10 CONTINUE
00302 *
00303       IF( N.EQ.0 )
00304      $   RETURN
00305 *
00306 *     Important constants
00307 *
00308       SMLNUM = SLAMCH( 'Safe minimum' )
00309       ULP = SLAMCH( 'Precision' )
00310       ULPINV = ONE / ULP
00311 *
00312 *     Perform tests (1)-(13)
00313 *
00314       SELOPT = 0
00315       DO 90 ISORT = 0, 1
00316          IF( ISORT.EQ.0 ) THEN
00317             SORT = 'N'
00318             RSUB = 0
00319          ELSE
00320             SORT = 'S'
00321             RSUB = 6
00322          END IF
00323 *
00324 *        Compute Schur form and Schur vectors, and test them
00325 *
00326          CALL CLACPY( 'F', N, N, A, LDA, H, LDA )
00327          CALL CGEESX( 'V', SORT, CSLECT, 'N', N, H, LDA, SDIM, W, VS,
00328      $                LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
00329      $                IINFO )
00330          IF( IINFO.NE.0 ) THEN
00331             RESULT( 1+RSUB ) = ULPINV
00332             IF( JTYPE.NE.22 ) THEN
00333                WRITE( NOUNIT, FMT = 9998 )'CGEESX1', IINFO, N, JTYPE,
00334      $            ISEED
00335             ELSE
00336                WRITE( NOUNIT, FMT = 9999 )'CGEESX1', IINFO, N,
00337      $            ISEED( 1 )
00338             END IF
00339             INFO = ABS( IINFO )
00340             RETURN
00341          END IF
00342          IF( ISORT.EQ.0 ) THEN
00343             CALL CCOPY( N, W, 1, WTMP, 1 )
00344          END IF
00345 *
00346 *        Do Test (1) or Test (7)
00347 *
00348          RESULT( 1+RSUB ) = ZERO
00349          DO 30 J = 1, N - 1
00350             DO 20 I = J + 1, N
00351                IF( H( I, J ).NE.CZERO )
00352      $            RESULT( 1+RSUB ) = ULPINV
00353    20       CONTINUE
00354    30    CONTINUE
00355 *
00356 *        Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
00357 *
00358 *        Copy A to VS1, used as workspace
00359 *
00360          CALL CLACPY( ' ', N, N, A, LDA, VS1, LDVS )
00361 *
00362 *        Compute Q*H and store in HT.
00363 *
00364          CALL CGEMM( 'No transpose', 'No transpose', N, N, N, CONE, VS,
00365      $               LDVS, H, LDA, CZERO, HT, LDA )
00366 *
00367 *        Compute A - Q*H*Q'
00368 *
00369          CALL CGEMM( 'No transpose', 'Conjugate transpose', N, N, N,
00370      $               -CONE, HT, LDA, VS, LDVS, CONE, VS1, LDVS )
00371 *
00372          ANORM = MAX( CLANGE( '1', N, N, A, LDA, RWORK ), SMLNUM )
00373          WNORM = CLANGE( '1', N, N, VS1, LDVS, RWORK )
00374 *
00375          IF( ANORM.GT.WNORM ) THEN
00376             RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP )
00377          ELSE
00378             IF( ANORM.LT.ONE ) THEN
00379                RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) /
00380      $                            ( N*ULP )
00381             ELSE
00382                RESULT( 2+RSUB ) = MIN( WNORM / ANORM, REAL( N ) ) /
00383      $                            ( N*ULP )
00384             END IF
00385          END IF
00386 *
00387 *        Test (3) or (9):  Compute norm( I - Q'*Q ) / ( N * ULP )
00388 *
00389          CALL CUNT01( 'Columns', N, N, VS, LDVS, WORK, LWORK, RWORK,
00390      $                RESULT( 3+RSUB ) )
00391 *
00392 *        Do Test (4) or Test (10)
00393 *
00394          RESULT( 4+RSUB ) = ZERO
00395          DO 40 I = 1, N
00396             IF( H( I, I ).NE.W( I ) )
00397      $         RESULT( 4+RSUB ) = ULPINV
00398    40    CONTINUE
00399 *
00400 *        Do Test (5) or Test (11)
00401 *
00402          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00403          CALL CGEESX( 'N', SORT, CSLECT, 'N', N, HT, LDA, SDIM, WT, VS,
00404      $                LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
00405      $                IINFO )
00406          IF( IINFO.NE.0 ) THEN
00407             RESULT( 5+RSUB ) = ULPINV
00408             IF( JTYPE.NE.22 ) THEN
00409                WRITE( NOUNIT, FMT = 9998 )'CGEESX2', IINFO, N, JTYPE,
00410      $            ISEED
00411             ELSE
00412                WRITE( NOUNIT, FMT = 9999 )'CGEESX2', IINFO, N,
00413      $            ISEED( 1 )
00414             END IF
00415             INFO = ABS( IINFO )
00416             GO TO 220
00417          END IF
00418 *
00419          RESULT( 5+RSUB ) = ZERO
00420          DO 60 J = 1, N
00421             DO 50 I = 1, N
00422                IF( H( I, J ).NE.HT( I, J ) )
00423      $            RESULT( 5+RSUB ) = ULPINV
00424    50       CONTINUE
00425    60    CONTINUE
00426 *
00427 *        Do Test (6) or Test (12)
00428 *
00429          RESULT( 6+RSUB ) = ZERO
00430          DO 70 I = 1, N
00431             IF( W( I ).NE.WT( I ) )
00432      $         RESULT( 6+RSUB ) = ULPINV
00433    70    CONTINUE
00434 *
00435 *        Do Test (13)
00436 *
00437          IF( ISORT.EQ.1 ) THEN
00438             RESULT( 13 ) = ZERO
00439             KNTEIG = 0
00440             DO 80 I = 1, N
00441                IF( CSLECT( W( I ) ) )
00442      $            KNTEIG = KNTEIG + 1
00443                IF( I.LT.N ) THEN
00444                   IF( CSLECT( W( I+1 ) ) .AND.
00445      $                ( .NOT.CSLECT( W( I ) ) ) )RESULT( 13 ) = ULPINV
00446                END IF
00447    80       CONTINUE
00448             IF( SDIM.NE.KNTEIG )
00449      $         RESULT( 13 ) = ULPINV
00450          END IF
00451 *
00452    90 CONTINUE
00453 *
00454 *     If there is enough workspace, perform tests (14) and (15)
00455 *     as well as (10) through (13)
00456 *
00457       IF( LWORK.GE.( N*( N+1 ) ) / 2 ) THEN
00458 *
00459 *        Compute both RCONDE and RCONDV with VS
00460 *
00461          SORT = 'S'
00462          RESULT( 14 ) = ZERO
00463          RESULT( 15 ) = ZERO
00464          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00465          CALL CGEESX( 'V', SORT, CSLECT, 'B', N, HT, LDA, SDIM1, WT,
00466      $                VS1, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
00467      $                BWORK, IINFO )
00468          IF( IINFO.NE.0 ) THEN
00469             RESULT( 14 ) = ULPINV
00470             RESULT( 15 ) = ULPINV
00471             IF( JTYPE.NE.22 ) THEN
00472                WRITE( NOUNIT, FMT = 9998 )'CGEESX3', IINFO, N, JTYPE,
00473      $            ISEED
00474             ELSE
00475                WRITE( NOUNIT, FMT = 9999 )'CGEESX3', IINFO, N,
00476      $            ISEED( 1 )
00477             END IF
00478             INFO = ABS( IINFO )
00479             GO TO 220
00480          END IF
00481 *
00482 *        Perform tests (10), (11), (12), and (13)
00483 *
00484          DO 110 I = 1, N
00485             IF( W( I ).NE.WT( I ) )
00486      $         RESULT( 10 ) = ULPINV
00487             DO 100 J = 1, N
00488                IF( H( I, J ).NE.HT( I, J ) )
00489      $            RESULT( 11 ) = ULPINV
00490                IF( VS( I, J ).NE.VS1( I, J ) )
00491      $            RESULT( 12 ) = ULPINV
00492   100       CONTINUE
00493   110    CONTINUE
00494          IF( SDIM.NE.SDIM1 )
00495      $      RESULT( 13 ) = ULPINV
00496 *
00497 *        Compute both RCONDE and RCONDV without VS, and compare
00498 *
00499          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00500          CALL CGEESX( 'N', SORT, CSLECT, 'B', N, HT, LDA, SDIM1, WT,
00501      $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
00502      $                BWORK, IINFO )
00503          IF( IINFO.NE.0 ) THEN
00504             RESULT( 14 ) = ULPINV
00505             RESULT( 15 ) = ULPINV
00506             IF( JTYPE.NE.22 ) THEN
00507                WRITE( NOUNIT, FMT = 9998 )'CGEESX4', IINFO, N, JTYPE,
00508      $            ISEED
00509             ELSE
00510                WRITE( NOUNIT, FMT = 9999 )'CGEESX4', IINFO, N,
00511      $            ISEED( 1 )
00512             END IF
00513             INFO = ABS( IINFO )
00514             GO TO 220
00515          END IF
00516 *
00517 *        Perform tests (14) and (15)
00518 *
00519          IF( RCNDE1.NE.RCONDE )
00520      $      RESULT( 14 ) = ULPINV
00521          IF( RCNDV1.NE.RCONDV )
00522      $      RESULT( 15 ) = ULPINV
00523 *
00524 *        Perform tests (10), (11), (12), and (13)
00525 *
00526          DO 130 I = 1, N
00527             IF( W( I ).NE.WT( I ) )
00528      $         RESULT( 10 ) = ULPINV
00529             DO 120 J = 1, N
00530                IF( H( I, J ).NE.HT( I, J ) )
00531      $            RESULT( 11 ) = ULPINV
00532                IF( VS( I, J ).NE.VS1( I, J ) )
00533      $            RESULT( 12 ) = ULPINV
00534   120       CONTINUE
00535   130    CONTINUE
00536          IF( SDIM.NE.SDIM1 )
00537      $      RESULT( 13 ) = ULPINV
00538 *
00539 *        Compute RCONDE with VS, and compare
00540 *
00541          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00542          CALL CGEESX( 'V', SORT, CSLECT, 'E', N, HT, LDA, SDIM1, WT,
00543      $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
00544      $                BWORK, IINFO )
00545          IF( IINFO.NE.0 ) THEN
00546             RESULT( 14 ) = ULPINV
00547             IF( JTYPE.NE.22 ) THEN
00548                WRITE( NOUNIT, FMT = 9998 )'CGEESX5', IINFO, N, JTYPE,
00549      $            ISEED
00550             ELSE
00551                WRITE( NOUNIT, FMT = 9999 )'CGEESX5', IINFO, N,
00552      $            ISEED( 1 )
00553             END IF
00554             INFO = ABS( IINFO )
00555             GO TO 220
00556          END IF
00557 *
00558 *        Perform test (14)
00559 *
00560          IF( RCNDE1.NE.RCONDE )
00561      $      RESULT( 14 ) = ULPINV
00562 *
00563 *        Perform tests (10), (11), (12), and (13)
00564 *
00565          DO 150 I = 1, N
00566             IF( W( I ).NE.WT( I ) )
00567      $         RESULT( 10 ) = ULPINV
00568             DO 140 J = 1, N
00569                IF( H( I, J ).NE.HT( I, J ) )
00570      $            RESULT( 11 ) = ULPINV
00571                IF( VS( I, J ).NE.VS1( I, J ) )
00572      $            RESULT( 12 ) = ULPINV
00573   140       CONTINUE
00574   150    CONTINUE
00575          IF( SDIM.NE.SDIM1 )
00576      $      RESULT( 13 ) = ULPINV
00577 *
00578 *        Compute RCONDE without VS, and compare
00579 *
00580          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00581          CALL CGEESX( 'N', SORT, CSLECT, 'E', N, HT, LDA, SDIM1, WT,
00582      $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
00583      $                BWORK, IINFO )
00584          IF( IINFO.NE.0 ) THEN
00585             RESULT( 14 ) = ULPINV
00586             IF( JTYPE.NE.22 ) THEN
00587                WRITE( NOUNIT, FMT = 9998 )'CGEESX6', IINFO, N, JTYPE,
00588      $            ISEED
00589             ELSE
00590                WRITE( NOUNIT, FMT = 9999 )'CGEESX6', IINFO, N,
00591      $            ISEED( 1 )
00592             END IF
00593             INFO = ABS( IINFO )
00594             GO TO 220
00595          END IF
00596 *
00597 *        Perform test (14)
00598 *
00599          IF( RCNDE1.NE.RCONDE )
00600      $      RESULT( 14 ) = ULPINV
00601 *
00602 *        Perform tests (10), (11), (12), and (13)
00603 *
00604          DO 170 I = 1, N
00605             IF( W( I ).NE.WT( I ) )
00606      $         RESULT( 10 ) = ULPINV
00607             DO 160 J = 1, N
00608                IF( H( I, J ).NE.HT( I, J ) )
00609      $            RESULT( 11 ) = ULPINV
00610                IF( VS( I, J ).NE.VS1( I, J ) )
00611      $            RESULT( 12 ) = ULPINV
00612   160       CONTINUE
00613   170    CONTINUE
00614          IF( SDIM.NE.SDIM1 )
00615      $      RESULT( 13 ) = ULPINV
00616 *
00617 *        Compute RCONDV with VS, and compare
00618 *
00619          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00620          CALL CGEESX( 'V', SORT, CSLECT, 'V', N, HT, LDA, SDIM1, WT,
00621      $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
00622      $                BWORK, IINFO )
00623          IF( IINFO.NE.0 ) THEN
00624             RESULT( 15 ) = ULPINV
00625             IF( JTYPE.NE.22 ) THEN
00626                WRITE( NOUNIT, FMT = 9998 )'CGEESX7', IINFO, N, JTYPE,
00627      $            ISEED
00628             ELSE
00629                WRITE( NOUNIT, FMT = 9999 )'CGEESX7', IINFO, N,
00630      $            ISEED( 1 )
00631             END IF
00632             INFO = ABS( IINFO )
00633             GO TO 220
00634          END IF
00635 *
00636 *        Perform test (15)
00637 *
00638          IF( RCNDV1.NE.RCONDV )
00639      $      RESULT( 15 ) = ULPINV
00640 *
00641 *        Perform tests (10), (11), (12), and (13)
00642 *
00643          DO 190 I = 1, N
00644             IF( W( I ).NE.WT( I ) )
00645      $         RESULT( 10 ) = ULPINV
00646             DO 180 J = 1, N
00647                IF( H( I, J ).NE.HT( I, J ) )
00648      $            RESULT( 11 ) = ULPINV
00649                IF( VS( I, J ).NE.VS1( I, J ) )
00650      $            RESULT( 12 ) = ULPINV
00651   180       CONTINUE
00652   190    CONTINUE
00653          IF( SDIM.NE.SDIM1 )
00654      $      RESULT( 13 ) = ULPINV
00655 *
00656 *        Compute RCONDV without VS, and compare
00657 *
00658          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00659          CALL CGEESX( 'N', SORT, CSLECT, 'V', N, HT, LDA, SDIM1, WT,
00660      $                VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
00661      $                BWORK, IINFO )
00662          IF( IINFO.NE.0 ) THEN
00663             RESULT( 15 ) = ULPINV
00664             IF( JTYPE.NE.22 ) THEN
00665                WRITE( NOUNIT, FMT = 9998 )'CGEESX8', IINFO, N, JTYPE,
00666      $            ISEED
00667             ELSE
00668                WRITE( NOUNIT, FMT = 9999 )'CGEESX8', IINFO, N,
00669      $            ISEED( 1 )
00670             END IF
00671             INFO = ABS( IINFO )
00672             GO TO 220
00673          END IF
00674 *
00675 *        Perform test (15)
00676 *
00677          IF( RCNDV1.NE.RCONDV )
00678      $      RESULT( 15 ) = ULPINV
00679 *
00680 *        Perform tests (10), (11), (12), and (13)
00681 *
00682          DO 210 I = 1, N
00683             IF( W( I ).NE.WT( I ) )
00684      $         RESULT( 10 ) = ULPINV
00685             DO 200 J = 1, N
00686                IF( H( I, J ).NE.HT( I, J ) )
00687      $            RESULT( 11 ) = ULPINV
00688                IF( VS( I, J ).NE.VS1( I, J ) )
00689      $            RESULT( 12 ) = ULPINV
00690   200       CONTINUE
00691   210    CONTINUE
00692          IF( SDIM.NE.SDIM1 )
00693      $      RESULT( 13 ) = ULPINV
00694 *
00695       END IF
00696 *
00697   220 CONTINUE
00698 *
00699 *     If there are precomputed reciprocal condition numbers, compare
00700 *     computed values with them.
00701 *
00702       IF( COMP ) THEN
00703 *
00704 *        First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that
00705 *        the logical function CSLECT selects the eigenvalues specified
00706 *        by NSLCT, ISLCT and ISRT.
00707 *
00708          SELDIM = N
00709          SELOPT = 1
00710          EPS = MAX( ULP, EPSIN )
00711          DO 230 I = 1, N
00712             IPNT( I ) = I
00713             SELVAL( I ) = .FALSE.
00714             SELWR( I ) = REAL( WTMP( I ) )
00715             SELWI( I ) = AIMAG( WTMP( I ) )
00716   230    CONTINUE
00717          DO 250 I = 1, N - 1
00718             KMIN = I
00719             IF( ISRT.EQ.0 ) THEN
00720                VRIMIN = REAL( WTMP( I ) )
00721             ELSE
00722                VRIMIN = AIMAG( WTMP( I ) )
00723             END IF
00724             DO 240 J = I + 1, N
00725                IF( ISRT.EQ.0 ) THEN
00726                   VRICMP = REAL( WTMP( J ) )
00727                ELSE
00728                   VRICMP = AIMAG( WTMP( J ) )
00729                END IF
00730                IF( VRICMP.LT.VRIMIN ) THEN
00731                   KMIN = J
00732                   VRIMIN = VRICMP
00733                END IF
00734   240       CONTINUE
00735             CTMP = WTMP( KMIN )
00736             WTMP( KMIN ) = WTMP( I )
00737             WTMP( I ) = CTMP
00738             ITMP = IPNT( I )
00739             IPNT( I ) = IPNT( KMIN )
00740             IPNT( KMIN ) = ITMP
00741   250    CONTINUE
00742          DO 260 I = 1, NSLCT
00743             SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE.
00744   260    CONTINUE
00745 *
00746 *        Compute condition numbers
00747 *
00748          CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
00749          CALL CGEESX( 'N', 'S', CSLECT, 'B', N, HT, LDA, SDIM1, WT, VS1,
00750      $                LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
00751      $                IINFO )
00752          IF( IINFO.NE.0 ) THEN
00753             RESULT( 16 ) = ULPINV
00754             RESULT( 17 ) = ULPINV
00755             WRITE( NOUNIT, FMT = 9999 )'CGEESX9', IINFO, N, ISEED( 1 )
00756             INFO = ABS( IINFO )
00757             GO TO 270
00758          END IF
00759 *
00760 *        Compare condition number for average of selected eigenvalues
00761 *        taking its condition number into account
00762 *
00763          ANORM = CLANGE( '1', N, N, A, LDA, RWORK )
00764          V = MAX( REAL( N )*EPS*ANORM, SMLNUM )
00765          IF( ANORM.EQ.ZERO )
00766      $      V = ONE
00767          IF( V.GT.RCONDV ) THEN
00768             TOL = ONE
00769          ELSE
00770             TOL = V / RCONDV
00771          END IF
00772          IF( V.GT.RCDVIN ) THEN
00773             TOLIN = ONE
00774          ELSE
00775             TOLIN = V / RCDVIN
00776          END IF
00777          TOL = MAX( TOL, SMLNUM / EPS )
00778          TOLIN = MAX( TOLIN, SMLNUM / EPS )
00779          IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN
00780             RESULT( 16 ) = ULPINV
00781          ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN
00782             RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
00783          ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN
00784             RESULT( 16 ) = ULPINV
00785          ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN
00786             RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN )
00787          ELSE
00788             RESULT( 16 ) = ONE
00789          END IF
00790 *
00791 *        Compare condition numbers for right invariant subspace
00792 *        taking its condition number into account
00793 *
00794          IF( V.GT.RCONDV*RCONDE ) THEN
00795             TOL = RCONDV
00796          ELSE
00797             TOL = V / RCONDE
00798          END IF
00799          IF( V.GT.RCDVIN*RCDEIN ) THEN
00800             TOLIN = RCDVIN
00801          ELSE
00802             TOLIN = V / RCDEIN
00803          END IF
00804          TOL = MAX( TOL, SMLNUM / EPS )
00805          TOLIN = MAX( TOLIN, SMLNUM / EPS )
00806          IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN
00807             RESULT( 17 ) = ULPINV
00808          ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN
00809             RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
00810          ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN
00811             RESULT( 17 ) = ULPINV
00812          ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN
00813             RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN )
00814          ELSE
00815             RESULT( 17 ) = ONE
00816          END IF
00817 *
00818   270    CONTINUE
00819 *
00820       END IF
00821 *
00822  9999 FORMAT( ' CGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00823      $      I6, ', INPUT EXAMPLE NUMBER = ', I4 )
00824  9998 FORMAT( ' CGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00825      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
00826 *
00827       RETURN
00828 *
00829 *     End of CGET24
00830 *
00831       END