testdata/lapack/TESTING/EIG/dget24.f
*> \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.
*
*> \date December 2016
*
*> \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 (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
* .. 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