dget24.f

Go to the documentation of this file.

*> \brief \b DGET24
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
*                          H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
*                          LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
*                          RESULT, WORK, LWORK, IWORK, BWORK, INFO )
*
*       .. Scalar Arguments ..
*       LOGICAL            COMP
*       INTEGER            INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
*       DOUBLE PRECISION   RCDEIN, RCDVIN, THRESH
*       ..
*       .. Array Arguments ..
*       LOGICAL            BWORK( * )
*       INTEGER            ISEED( 4 ), ISLCT( * ), IWORK( * )
*       DOUBLE PRECISION   A( LDA, * ), H( LDA, * ), HT( LDA, * ),
*      $                   RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
*      $                   WI( * ), WIT( * ), WITMP( * ), WORK( * ),
*      $                   WR( * ), WRT( * ), WRTMP( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    DGET24 checks the nonsymmetric eigenvalue (Schur form) problem
*>    expert driver DGEESX.
*>
*>    If COMP = .FALSE., the first 13 of the following tests will be
*>    be performed on the input matrix A, and also tests 14 and 15
*>    if LWORK is sufficiently large.
*>    If COMP = .TRUE., all 17 test will be performed.
*>
*>    (1)     0 if T is in Schur form, 1/ulp otherwise
*>           (no sorting of eigenvalues)
*>
*>    (2)     | A - VS T VS' | / ( n |A| ulp )
*>
*>      Here VS is the matrix of Schur eigenvectors, and T is in Schur
*>      form  (no sorting of eigenvalues).
*>
*>    (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
*>
*>    (4)     0     if WR+sqrt(-1)*WI are eigenvalues of T
*>            1/ulp otherwise
*>            (no sorting of eigenvalues)
*>
*>    (5)     0     if T(with VS) = T(without VS),
*>            1/ulp otherwise
*>            (no sorting of eigenvalues)
*>
*>    (6)     0     if eigenvalues(with VS) = eigenvalues(without VS),
*>            1/ulp otherwise
*>            (no sorting of eigenvalues)
*>
*>    (7)     0 if T is in Schur form, 1/ulp otherwise
*>            (with sorting of eigenvalues)
*>
*>    (8)     | A - VS T VS' | / ( n |A| ulp )
*>
*>      Here VS is the matrix of Schur eigenvectors, and T is in Schur
*>      form  (with sorting of eigenvalues).
*>
*>    (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
*>
*>    (10)    0     if WR+sqrt(-1)*WI are eigenvalues of T
*>            1/ulp otherwise
*>            If workspace sufficient, also compare WR, WI with and
*>            without reciprocal condition numbers
*>            (with sorting of eigenvalues)
*>
*>    (11)    0     if T(with VS) = T(without VS),
*>            1/ulp otherwise
*>            If workspace sufficient, also compare T with and without
*>            reciprocal condition numbers
*>            (with sorting of eigenvalues)
*>
*>    (12)    0     if eigenvalues(with VS) = eigenvalues(without VS),
*>            1/ulp otherwise
*>            If workspace sufficient, also compare VS with and without
*>            reciprocal condition numbers
*>            (with sorting of eigenvalues)
*>
*>    (13)    if sorting worked and SDIM is the number of
*>            eigenvalues which were SELECTed
*>            If workspace sufficient, also compare SDIM with and
*>            without reciprocal condition numbers
*>
*>    (14)    if RCONDE the same no matter if VS and/or RCONDV computed
*>
*>    (15)    if RCONDV the same no matter if VS and/or RCONDE computed
*>
*>    (16)  |RCONDE - RCDEIN| / cond(RCONDE)
*>
*>       RCONDE is the reciprocal average eigenvalue condition number
*>       computed by DGEESX and RCDEIN (the precomputed true value)
*>       is supplied as input.  cond(RCONDE) is the condition number
*>       of RCONDE, and takes errors in computing RCONDE into account,
*>       so that the resulting quantity should be O(ULP). cond(RCONDE)
*>       is essentially given by norm(A)/RCONDV.
*>
*>    (17)  |RCONDV - RCDVIN| / cond(RCONDV)
*>
*>       RCONDV is the reciprocal right invariant subspace condition
*>       number computed by DGEESX and RCDVIN (the precomputed true
*>       value) is supplied as input. cond(RCONDV) is the condition
*>       number of RCONDV, and takes errors in computing RCONDV into
*>       account, so that the resulting quantity should be O(ULP).
*>       cond(RCONDV) is essentially given by norm(A)/RCONDE.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] COMP
*> \verbatim
*>          COMP is LOGICAL
*>          COMP describes which input tests to perform:
*>            = .FALSE. if the computed condition numbers are not to
*>                      be tested against RCDVIN and RCDEIN
*>            = .TRUE.  if they are to be compared
*> \endverbatim
*>
*> \param[in] JTYPE
*> \verbatim
*>          JTYPE is INTEGER
*>          Type of input matrix. Used to label output if error occurs.
*> \endverbatim
*>
*> \param[in] ISEED
*> \verbatim
*>          ISEED is INTEGER array, dimension (4)
*>          If COMP = .FALSE., the random number generator seed
*>          used to produce matrix.
*>          If COMP = .TRUE., ISEED(1) = the number of the example.
*>          Used to label output if error occurs.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*>          THRESH is DOUBLE PRECISION
*>          A test will count as "failed" if the "error", computed as
*>          described above, exceeds THRESH.  Note that the error
*>          is scaled to be O(1), so THRESH should be a reasonably
*>          small multiple of 1, e.g., 10 or 100.  In particular,
*>          it should not depend on the precision (single vs. double)
*>          or the size of the matrix.  It must be at least zero.
*> \endverbatim
*>
*> \param[in] NOUNIT
*> \verbatim
*>          NOUNIT is INTEGER
*>          The FORTRAN unit number for printing out error messages
*>          (e.g., if a routine returns INFO not equal to 0.)
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The dimension of A. N must be at least 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, N)
*>          Used to hold the matrix whose eigenvalues are to be
*>          computed.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of A, and H. LDA must be at
*>          least 1 and at least N.
*> \endverbatim
*>
*> \param[out] H
*> \verbatim
*>          H is DOUBLE PRECISION array, dimension (LDA, N)
*>          Another copy of the test matrix A, modified by DGEESX.
*> \endverbatim
*>
*> \param[out] HT
*> \verbatim
*>          HT is DOUBLE PRECISION array, dimension (LDA, N)
*>          Yet another copy of the test matrix A, modified by DGEESX.
*> \endverbatim
*>
*> \param[out] WR
*> \verbatim
*>          WR is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] WI
*> \verbatim
*>          WI is DOUBLE PRECISION array, dimension (N)
*>
*>          The real and imaginary parts of the eigenvalues of A.
*>          On exit, WR + WI*i are the eigenvalues of the matrix in A.
*> \endverbatim
*>
*> \param[out] WRT
*> \verbatim
*>          WRT is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] WIT
*> \verbatim
*>          WIT is DOUBLE PRECISION array, dimension (N)
*>
*>          Like WR, WI, these arrays contain the eigenvalues of A,
*>          but those computed when DGEESX only computes a partial
*>          eigendecomposition, i.e. not Schur vectors
*> \endverbatim
*>
*> \param[out] WRTMP
*> \verbatim
*>          WRTMP is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] WITMP
*> \verbatim
*>          WITMP is DOUBLE PRECISION array, dimension (N)
*>
*>          Like WR, WI, these arrays contain the eigenvalues of A,
*>          but sorted by increasing real part.
*> \endverbatim
*>
*> \param[out] VS
*> \verbatim
*>          VS is DOUBLE PRECISION array, dimension (LDVS, N)
*>          VS holds the computed Schur vectors.
*> \endverbatim
*>
*> \param[in] LDVS
*> \verbatim
*>          LDVS is INTEGER
*>          Leading dimension of VS. Must be at least max(1, N).
*> \endverbatim
*>
*> \param[out] VS1
*> \verbatim
*>          VS1 is DOUBLE PRECISION array, dimension (LDVS, N)
*>          VS1 holds another copy of the computed Schur vectors.
*> \endverbatim
*>
*> \param[in] RCDEIN
*> \verbatim
*>          RCDEIN is DOUBLE PRECISION
*>          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
*>          condition number for the average of selected eigenvalues.
*> \endverbatim
*>
*> \param[in] RCDVIN
*> \verbatim
*>          RCDVIN is DOUBLE PRECISION
*>          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
*>          condition number for the selected right invariant subspace.
*> \endverbatim
*>
*> \param[in] NSLCT
*> \verbatim
*>          NSLCT is INTEGER
*>          When COMP = .TRUE. the number of selected eigenvalues
*>          corresponding to the precomputed values RCDEIN and RCDVIN.
*> \endverbatim
*>
*> \param[in] ISLCT
*> \verbatim
*>          ISLCT is INTEGER array, dimension (NSLCT)
*>          When COMP = .TRUE. ISLCT selects the eigenvalues of the
*>          input matrix corresponding to the precomputed values RCDEIN
*>          and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
*>          eigenvalue with the J-th largest real part is selected.
*>          Not referenced if COMP = .FALSE.
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*>          RESULT is DOUBLE PRECISION array, dimension (17)
*>          The values computed by the 17 tests described above.
*>          The values are currently limited to 1/ulp, to avoid
*>          overflow.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (LWORK)
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The number of entries in WORK to be passed to DGEESX. This
*>          must be at least 3*N, and N+N**2 if tests 14--16 are to
*>          be performed.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (N*N)
*> \endverbatim
*>
*> \param[out] BWORK
*> \verbatim
*>          BWORK is LOGICAL array, dimension (N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          If 0,  successful exit.
*>          If <0, input parameter -INFO had an incorrect value.
*>          If >0, DGEESX returned an error code, the absolute
*>                 value of which is returned.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE dget24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
     $                   H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
     $                   LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
     $                   RESULT, WORK, LWORK, IWORK, BWORK, INFO )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            COMP
      INTEGER            INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
      DOUBLE PRECISION   RCDEIN, RCDVIN, THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * )
      INTEGER            ISEED( 4 ), ISLCT( * ), IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), H( LDA, * ), HT( LDA, * ),
     $                   result( 17 ), vs( ldvs, * ), vs1( ldvs, * ),
     $                   wi( * ), wit( * ), witmp( * ), work( * ),
     $                   wr( * ), wrt( * ), wrtmp( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0d0, one = 1.0d0 )
      DOUBLE PRECISION   EPSIN
      parameter( epsin = 5.9605d-8 )
*     ..
*     .. Local Scalars ..
      CHARACTER          SORT
      INTEGER            I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
     $                   RSUB, SDIM, SDIM1
      DOUBLE PRECISION   ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
     $                   smlnum, tmp, tol, tolin, ulp, ulpinv, v, vimin,
     $                   vrmin, wnorm
*     ..
*     .. Local Arrays ..
      INTEGER            IPNT( 20 )
*     ..
*     .. Arrays in Common ..
      LOGICAL            SELVAL( 20 )
      DOUBLE PRECISION   SELWI( 20 ), SELWR( 20 )
*     ..
*     .. Scalars in Common ..
      INTEGER            SELDIM, SELOPT
*     ..
*     .. Common blocks ..
      COMMON             / sslct / selopt, seldim, selval, selwr, selwi
*     ..
*     .. External Functions ..
      LOGICAL            DSLECT
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           DSLECT, DLAMCH, DLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dgeesx, dgemm, dlacpy, dort01, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, max, min, sign, sqrt
*     ..
*     .. Executable Statements ..
*
*     Check for errors
*
      info = 0
      IF( thresh.LT.zero ) THEN
         info = -3
      ELSE IF( nounit.LE.0 ) THEN
         info = -5
      ELSE IF( n.LT.0 ) THEN
         info = -6
      ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
         info = -8
      ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n ) THEN
         info = -18
      ELSE IF( lwork.LT.3*n ) THEN
         info = -26
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DGET24', -info )
         RETURN
      END IF
*
*     Quick return if nothing to do
*
      DO 10 i = 1, 17
         result( i ) = -one
   10 CONTINUE
*
      IF( n.EQ.0 )
     $   RETURN
*
*     Important constants
*
      smlnum = dlamch( 'Safe minimum' )
      ulp = dlamch( 'Precision' )
      ulpinv = one / ulp
*
*     Perform tests (1)-(13)
*
      selopt = 0
      liwork = n*n
      DO 120 isort = 0, 1
         IF( isort.EQ.0 ) THEN
            sort = 'N'
            rsub = 0
         ELSE
            sort = 'S'
            rsub = 6
         END IF
*
*        Compute Schur form and Schur vectors, and test them
*
         CALL dlacpy( 'F', n, n, a, lda, h, lda )
         CALL dgeesx( 'V', sort, dslect, 'N', n, h, lda, sdim, wr, wi,
     $                vs, ldvs, rconde, rcondv, work, lwork, iwork,
     $                liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 1+rsub ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX1', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX1', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            RETURN
         END IF
         IF( isort.EQ.0 ) THEN
            CALL dcopy( n, wr, 1, wrtmp, 1 )
            CALL dcopy( n, wi, 1, witmp, 1 )
         END IF
*
*        Do Test (1) or Test (7)
*
         result( 1+rsub ) = zero
         DO 30 j = 1, n - 2
            DO 20 i = j + 2, n
               IF( h( i, j ).NE.zero )
     $            result( 1+rsub ) = ulpinv
   20       CONTINUE
   30    CONTINUE
         DO 40 i = 1, n - 2
            IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.zero )
     $         result( 1+rsub ) = ulpinv
   40    CONTINUE
         DO 50 i = 1, n - 1
            IF( h( i+1, i ).NE.zero ) THEN
               IF( h( i, i ).NE.h( i+1, i+1 ) .OR. h( i, i+1 ).EQ.
     $             zero .OR. sign( one, h( i+1, i ) ).EQ.
     $             sign( one, h( i, i+1 ) ) )result( 1+rsub ) = ulpinv
            END IF
   50    CONTINUE
*
*        Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
*
*        Copy A to VS1, used as workspace
*
         CALL dlacpy( ' ', n, n, a, lda, vs1, ldvs )
*
*        Compute Q*H and store in HT.
*
         CALL dgemm( 'No transpose', 'No transpose', n, n, n, one, vs,
     $               ldvs, h, lda, zero, ht, lda )
*
*        Compute A - Q*H*Q'
*
         CALL dgemm( 'No transpose', 'Transpose', n, n, n, -one, ht,
     $               lda, vs, ldvs, one, vs1, ldvs )
*
         anorm = max( dlange( '1', n, n, a, lda, work ), smlnum )
         wnorm = dlange( '1', n, n, vs1, ldvs, work )
*
         IF( anorm.GT.wnorm ) THEN
            result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
         ELSE
            IF( anorm.LT.one ) THEN
               result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
     $                            ( n*ulp )
            ELSE
               result( 2+rsub ) = min( wnorm / anorm, dble( n ) ) /
     $                            ( n*ulp )
            END IF
         END IF
*
*        Test (3) or (9):  Compute norm( I - Q'*Q ) / ( N * ULP )
*
         CALL dort01( 'Columns', n, n, vs, ldvs, work, lwork,
     $                result( 3+rsub ) )
*
*        Do Test (4) or Test (10)
*
         result( 4+rsub ) = zero
         DO 60 i = 1, n
            IF( h( i, i ).NE.wr( i ) )
     $         result( 4+rsub ) = ulpinv
   60    CONTINUE
         IF( n.GT.1 ) THEN
            IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
     $         result( 4+rsub ) = ulpinv
            IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
     $         result( 4+rsub ) = ulpinv
         END IF
         DO 70 i = 1, n - 1
            IF( h( i+1, i ).NE.zero ) THEN
               tmp = sqrt( abs( h( i+1, i ) ) )*
     $               sqrt( abs( h( i, i+1 ) ) )
               result( 4+rsub ) = max( result( 4+rsub ),
     $                            abs( wi( i )-tmp ) /
     $                            max( ulp*tmp, smlnum ) )
               result( 4+rsub ) = max( result( 4+rsub ),
     $                            abs( wi( i+1 )+tmp ) /
     $                            max( ulp*tmp, smlnum ) )
            ELSE IF( i.GT.1 ) THEN
               IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.zero .AND.
     $             wi( i ).NE.zero )result( 4+rsub ) = ulpinv
            END IF
   70    CONTINUE
*
*        Do Test (5) or Test (11)
*
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'N', sort, dslect, 'N', n, ht, lda, sdim, wrt,
     $                wit, vs, ldvs, rconde, rcondv, work, lwork, iwork,
     $                liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 5+rsub ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX2', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX2', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            GO TO 250
         END IF
*
         result( 5+rsub ) = zero
         DO 90 j = 1, n
            DO 80 i = 1, n
               IF( h( i, j ).NE.ht( i, j ) )
     $            result( 5+rsub ) = ulpinv
   80       CONTINUE
   90    CONTINUE
*
*        Do Test (6) or Test (12)
*
         result( 6+rsub ) = zero
         DO 100 i = 1, n
            IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
     $         result( 6+rsub ) = ulpinv
  100    CONTINUE
*
*        Do Test (13)
*
         IF( isort.EQ.1 ) THEN
            result( 13 ) = zero
            knteig = 0
            DO 110 i = 1, n
               IF( dslect( wr( i ), wi( i ) ) .OR.
     $             dslect( wr( i ), -wi( i ) ) )knteig = knteig + 1
               IF( i.LT.n ) THEN
                  IF( ( dslect( wr( i+1 ), wi( i+1 ) ) .OR.
     $                dslect( wr( i+1 ), -wi( i+1 ) ) ) .AND.
     $                ( .NOT.( dslect( wr( i ),
     $                wi( i ) ) .OR. dslect( wr( i ),
     $                -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )result( 13 )
     $                = ulpinv
               END IF
  110       CONTINUE
            IF( sdim.NE.knteig )
     $         result( 13 ) = ulpinv
         END IF
*
  120 CONTINUE
*
*     If there is enough workspace, perform tests (14) and (15)
*     as well as (10) through (13)
*
      IF( lwork.GE.n+( n*n ) / 2 ) THEN
*
*        Compute both RCONDE and RCONDV with VS
*
         sort = 'S'
         result( 14 ) = zero
         result( 15 ) = zero
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'V', sort, dslect, 'B', n, ht, lda, sdim1, wrt,
     $                wit, vs1, ldvs, rconde, rcondv, work, lwork,
     $                iwork, liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 14 ) = ulpinv
            result( 15 ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX3', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX3', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            GO TO 250
         END IF
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 140 i = 1, n
            IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
     $         result( 10 ) = ulpinv
            DO 130 j = 1, n
               IF( h( i, j ).NE.ht( i, j ) )
     $            result( 11 ) = ulpinv
               IF( vs( i, j ).NE.vs1( i, j ) )
     $            result( 12 ) = ulpinv
  130       CONTINUE
  140    CONTINUE
         IF( sdim.NE.sdim1 )
     $      result( 13 ) = ulpinv
*
*        Compute both RCONDE and RCONDV without VS, and compare
*
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'N', sort, dslect, 'B', n, ht, lda, sdim1, wrt,
     $                wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
     $                iwork, liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 14 ) = ulpinv
            result( 15 ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX4', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX4', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            GO TO 250
         END IF
*
*        Perform tests (14) and (15)
*
         IF( rcnde1.NE.rconde )
     $      result( 14 ) = ulpinv
         IF( rcndv1.NE.rcondv )
     $      result( 15 ) = ulpinv
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 160 i = 1, n
            IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
     $         result( 10 ) = ulpinv
            DO 150 j = 1, n
               IF( h( i, j ).NE.ht( i, j ) )
     $            result( 11 ) = ulpinv
               IF( vs( i, j ).NE.vs1( i, j ) )
     $            result( 12 ) = ulpinv
  150       CONTINUE
  160    CONTINUE
         IF( sdim.NE.sdim1 )
     $      result( 13 ) = ulpinv
*
*        Compute RCONDE with VS, and compare
*
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'V', sort, dslect, 'E', n, ht, lda, sdim1, wrt,
     $                wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
     $                iwork, liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 14 ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX5', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX5', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            GO TO 250
         END IF
*
*        Perform test (14)
*
         IF( rcnde1.NE.rconde )
     $      result( 14 ) = ulpinv
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 180 i = 1, n
            IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
     $         result( 10 ) = ulpinv
            DO 170 j = 1, n
               IF( h( i, j ).NE.ht( i, j ) )
     $            result( 11 ) = ulpinv
               IF( vs( i, j ).NE.vs1( i, j ) )
     $            result( 12 ) = ulpinv
  170       CONTINUE
  180    CONTINUE
         IF( sdim.NE.sdim1 )
     $      result( 13 ) = ulpinv
*
*        Compute RCONDE without VS, and compare
*
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'N', sort, dslect, 'E', n, ht, lda, sdim1, wrt,
     $                wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
     $                iwork, liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 14 ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX6', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX6', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            GO TO 250
         END IF
*
*        Perform test (14)
*
         IF( rcnde1.NE.rconde )
     $      result( 14 ) = ulpinv
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 200 i = 1, n
            IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
     $         result( 10 ) = ulpinv
            DO 190 j = 1, n
               IF( h( i, j ).NE.ht( i, j ) )
     $            result( 11 ) = ulpinv
               IF( vs( i, j ).NE.vs1( i, j ) )
     $            result( 12 ) = ulpinv
  190       CONTINUE
  200    CONTINUE
         IF( sdim.NE.sdim1 )
     $      result( 13 ) = ulpinv
*
*        Compute RCONDV with VS, and compare
*
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'V', sort, dslect, 'V', n, ht, lda, sdim1, wrt,
     $                wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
     $                iwork, liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 15 ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX7', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX7', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            GO TO 250
         END IF
*
*        Perform test (15)
*
         IF( rcndv1.NE.rcondv )
     $      result( 15 ) = ulpinv
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 220 i = 1, n
            IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
     $         result( 10 ) = ulpinv
            DO 210 j = 1, n
               IF( h( i, j ).NE.ht( i, j ) )
     $            result( 11 ) = ulpinv
               IF( vs( i, j ).NE.vs1( i, j ) )
     $            result( 12 ) = ulpinv
  210       CONTINUE
  220    CONTINUE
         IF( sdim.NE.sdim1 )
     $      result( 13 ) = ulpinv
*
*        Compute RCONDV without VS, and compare
*
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'N', sort, dslect, 'V', n, ht, lda, sdim1, wrt,
     $                wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
     $                iwork, liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 15 ) = ulpinv
            IF( jtype.NE.22 ) THEN
               WRITE( nounit, fmt = 9998 )'DGEESX8', iinfo, n, jtype,
     $            iseed
            ELSE
               WRITE( nounit, fmt = 9999 )'DGEESX8', iinfo, n,
     $            iseed( 1 )
            END IF
            info = abs( iinfo )
            GO TO 250
         END IF
*
*        Perform test (15)
*
         IF( rcndv1.NE.rcondv )
     $      result( 15 ) = ulpinv
*
*        Perform tests (10), (11), (12), and (13)
*
         DO 240 i = 1, n
            IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
     $         result( 10 ) = ulpinv
            DO 230 j = 1, n
               IF( h( i, j ).NE.ht( i, j ) )
     $            result( 11 ) = ulpinv
               IF( vs( i, j ).NE.vs1( i, j ) )
     $            result( 12 ) = ulpinv
  230       CONTINUE
  240    CONTINUE
         IF( sdim.NE.sdim1 )
     $      result( 13 ) = ulpinv
*
      END IF
*
  250 CONTINUE
*
*     If there are precomputed reciprocal condition numbers, compare
*     computed values with them.
*
      IF( comp ) THEN
*
*        First set up SELOPT, SELDIM, SELVAL, SELWR, and SELWI so that
*        the logical function DSLECT selects the eigenvalues specified
*        by NSLCT and ISLCT.
*
         seldim = n
         selopt = 1
         eps = max( ulp, epsin )
         DO 260 i = 1, n
            ipnt( i ) = i
            selval( i ) = .false.
            selwr( i ) = wrtmp( i )
            selwi( i ) = witmp( i )
  260    CONTINUE
         DO 280 i = 1, n - 1
            kmin = i
            vrmin = wrtmp( i )
            vimin = witmp( i )
            DO 270 j = i + 1, n
               IF( wrtmp( j ).LT.vrmin ) THEN
                  kmin = j
                  vrmin = wrtmp( j )
                  vimin = witmp( j )
               END IF
  270       CONTINUE
            wrtmp( kmin ) = wrtmp( i )
            witmp( kmin ) = witmp( i )
            wrtmp( i ) = vrmin
            witmp( i ) = vimin
            itmp = ipnt( i )
            ipnt( i ) = ipnt( kmin )
            ipnt( kmin ) = itmp
  280    CONTINUE
         DO 290 i = 1, nslct
            selval( ipnt( islct( i ) ) ) = .true.
  290    CONTINUE
*
*        Compute condition numbers
*
         CALL dlacpy( 'F', n, n, a, lda, ht, lda )
         CALL dgeesx( 'N', 'S', dslect, 'B', n, ht, lda, sdim1, wrt,
     $                wit, vs1, ldvs, rconde, rcondv, work, lwork,
     $                iwork, liwork, bwork, iinfo )
         IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
            result( 16 ) = ulpinv
            result( 17 ) = ulpinv
            WRITE( nounit, fmt = 9999 )'DGEESX9', iinfo, n, iseed( 1 )
            info = abs( iinfo )
            GO TO 300
         END IF
*
*        Compare condition number for average of selected eigenvalues
*        taking its condition number into account
*
         anorm = dlange( '1', n, n, a, lda, work )
         v = max( dble( n )*eps*anorm, smlnum )
         IF( anorm.EQ.zero )
     $      v = one
         IF( v.GT.rcondv ) THEN
            tol = one
         ELSE
            tol = v / rcondv
         END IF
         IF( v.GT.rcdvin ) THEN
            tolin = one
         ELSE
            tolin = v / rcdvin
         END IF
         tol = max( tol, smlnum / eps )
         tolin = max( tolin, smlnum / eps )
         IF( eps*( rcdein-tolin ).GT.rconde+tol ) THEN
            result( 16 ) = ulpinv
         ELSE IF( rcdein-tolin.GT.rconde+tol ) THEN
            result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
         ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) ) THEN
            result( 16 ) = ulpinv
         ELSE IF( rcdein+tolin.LT.rconde-tol ) THEN
            result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
         ELSE
            result( 16 ) = one
         END IF
*
*        Compare condition numbers for right invariant subspace
*        taking its condition number into account
*
         IF( v.GT.rcondv*rconde ) THEN
            tol = rcondv
         ELSE
            tol = v / rconde
         END IF
         IF( v.GT.rcdvin*rcdein ) THEN
            tolin = rcdvin
         ELSE
            tolin = v / rcdein
         END IF
         tol = max( tol, smlnum / eps )
         tolin = max( tolin, smlnum / eps )
         IF( eps*( rcdvin-tolin ).GT.rcondv+tol ) THEN
            result( 17 ) = ulpinv
         ELSE IF( rcdvin-tolin.GT.rcondv+tol ) THEN
            result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
         ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) ) THEN
            result( 17 ) = ulpinv
         ELSE IF( rcdvin+tolin.LT.rcondv-tol ) THEN
            result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
         ELSE
            result( 17 ) = one
         END IF
*
  300    CONTINUE
*
      END IF
*
 9999 FORMAT( ' DGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
     $      i6, ', INPUT EXAMPLE NUMBER = ', i4 )
 9998 FORMAT( ' DGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
     $      i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
*
      RETURN
*
*     End of DGET24
*
      END