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