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