LAPACK 3.3.0

dget24.f

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