You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
509 lines
18 KiB
509 lines
18 KiB
2 years ago
|
*> \brief \b ZEBCHVXX
|
||
|
|
*
|
||
|
|
* =========== DOCUMENTATION ===========
|
||
|
|
*
|
||
|
|
* Online html documentation available at
|
||
|
|
* http://www.netlib.org/lapack/explore-html/
|
||
|
|
*
|
||
|
|
* Definition:
|
||
|
|
* ===========
|
||
|
|
*
|
||
|
|
* SUBROUTINE ZEBCHVXX( THRESH, PATH )
|
||
|
|
*
|
||
|
|
* .. Scalar Arguments ..
|
||
|
|
* DOUBLE PRECISION THRESH
|
||
|
|
* CHARACTER*3 PATH
|
||
|
|
* ..
|
||
|
|
*
|
||
|
|
* Purpose
|
||
|
|
* ======
|
||
|
|
*
|
||
|
|
*> \details \b Purpose:
|
||
|
|
*> \verbatim
|
||
|
|
*>
|
||
|
|
*> ZEBCHVXX will run Z**SVXX on a series of Hilbert matrices and then
|
||
|
|
*> compare the error bounds returned by Z**SVXX to see if the returned
|
||
|
|
*> answer indeed falls within those bounds.
|
||
|
|
*>
|
||
|
|
*> Eight test ratios will be computed. The tests will pass if they are .LT.
|
||
|
|
*> THRESH. There are two cases that are determined by 1 / (SQRT( N ) * EPS).
|
||
|
|
*> If that value is .LE. to the component wise reciprocal condition number,
|
||
|
|
*> it uses the guaranteed case, other wise it uses the unguaranteed case.
|
||
|
|
*>
|
||
|
|
*> Test ratios:
|
||
|
|
*> Let Xc be X_computed and Xt be X_truth.
|
||
|
|
*> The norm used is the infinity norm.
|
||
|
|
*>
|
||
|
|
*> Let A be the guaranteed case and B be the unguaranteed case.
|
||
|
|
*>
|
||
|
|
*> 1. Normwise guaranteed forward error bound.
|
||
|
|
*> A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and
|
||
|
|
*> ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS.
|
||
|
|
*> If these conditions are met, the test ratio is set to be
|
||
|
|
*> ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
|
||
|
|
*> B: For this case, CGESVXX should just return 1. If it is less than
|
||
|
|
*> one, treat it the same as in 1A. Otherwise it fails. (Set test
|
||
|
|
*> ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?)
|
||
|
|
*>
|
||
|
|
*> 2. Componentwise guaranteed forward error bound.
|
||
|
|
*> A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i )
|
||
|
|
*> for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS.
|
||
|
|
*> If these conditions are met, the test ratio is set to be
|
||
|
|
*> ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
|
||
|
|
*> B: Same as normwise test ratio.
|
||
|
|
*>
|
||
|
|
*> 3. Backwards error.
|
||
|
|
*> A: The test ratio is set to BERR/EPS.
|
||
|
|
*> B: Same test ratio.
|
||
|
|
*>
|
||
|
|
*> 4. Reciprocal condition number.
|
||
|
|
*> A: A condition number is computed with Xt and compared with the one
|
||
|
|
*> returned from CGESVXX. Let RCONDc be the RCOND returned by CGESVXX
|
||
|
|
*> and RCONDt be the RCOND from the truth value. Test ratio is set to
|
||
|
|
*> MAX(RCONDc/RCONDt, RCONDt/RCONDc).
|
||
|
|
*> B: Test ratio is set to 1 / (EPS * RCONDc).
|
||
|
|
*>
|
||
|
|
*> 5. Reciprocal normwise condition number.
|
||
|
|
*> A: The test ratio is set to
|
||
|
|
*> MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )).
|
||
|
|
*> B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )).
|
||
|
|
*>
|
||
|
|
*> 6. Reciprocal componentwise condition number.
|
||
|
|
*> A: Test ratio is set to
|
||
|
|
*> MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )).
|
||
|
|
*> B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )).
|
||
|
|
*>
|
||
|
|
*> .. Parameters ..
|
||
|
|
*> NMAX is determined by the largest number in the inverse of the hilbert
|
||
|
|
*> matrix. Precision is exhausted when the largest entry in it is greater
|
||
|
|
*> than 2 to the power of the number of bits in the fraction of the data
|
||
|
|
*> type used plus one, which is 24 for single precision.
|
||
|
|
*> NMAX should be 6 for single and 11 for double.
|
||
|
|
*> \endverbatim
|
||
|
|
*
|
||
|
|
* Authors:
|
||
|
|
* ========
|
||
|
|
*
|
||
|
|
*> \author Univ. of Tennessee
|
||
|
|
*> \author Univ. of California Berkeley
|
||
|
|
*> \author Univ. of Colorado Denver
|
||
|
|
*> \author NAG Ltd.
|
||
|
|
*
|
||
|
|
*> \ingroup complex16_lin
|
||
|
|
*
|
||
|
|
* =====================================================================
|
||
|
|
SUBROUTINE ZEBCHVXX( THRESH, PATH )
|
||
|
|
IMPLICIT NONE
|
||
|
|
* .. Scalar Arguments ..
|
||
|
|
DOUBLE PRECISION THRESH
|
||
|
|
CHARACTER*3 PATH
|
||
|
|
|
||
|
|
INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
|
||
|
|
PARAMETER (NMAX = 10, NPARAMS = 2, NERRBND = 3,
|
||
|
|
$ NTESTS = 6)
|
||
|
|
|
||
|
|
* .. Local Scalars ..
|
||
|
|
INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA,
|
||
|
|
$ N_AUX_TESTS, LDAB, LDAFB
|
||
|
|
CHARACTER FACT, TRANS, UPLO, EQUED
|
||
|
|
CHARACTER*2 C2
|
||
|
|
CHARACTER(3) NGUAR, CGUAR
|
||
|
|
LOGICAL printed_guide
|
||
|
|
DOUBLE PRECISION NCOND, CCOND, M, NORMDIF, NORMT, RCOND,
|
||
|
|
$ RNORM, RINORM, SUMR, SUMRI, EPS,
|
||
|
|
$ BERR(NMAX), RPVGRW, ORCOND,
|
||
|
|
$ CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND,
|
||
|
|
$ CWISE_RCOND, NWISE_RCOND,
|
||
|
|
$ CONDTHRESH, ERRTHRESH
|
||
|
|
COMPLEX*16 ZDUM
|
||
|
|
|
||
|
|
* .. Local Arrays ..
|
||
|
|
DOUBLE PRECISION TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
|
||
|
|
$ S(NMAX),R(NMAX),C(NMAX),RWORK(3*NMAX),
|
||
|
|
$ DIFF(NMAX, NMAX),
|
||
|
|
$ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
|
||
|
|
INTEGER IPIV(NMAX)
|
||
|
|
COMPLEX*16 A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
|
||
|
|
$ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
|
||
|
|
$ ACOPY(NMAX, NMAX),
|
||
|
|
$ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
|
||
|
|
$ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
|
||
|
|
$ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX )
|
||
|
|
|
||
|
|
* .. External Functions ..
|
||
|
|
DOUBLE PRECISION DLAMCH
|
||
|
|
|
||
|
|
* .. External Subroutines ..
|
||
|
|
EXTERNAL ZLAHILB, ZGESVXX, ZPOSVXX, ZSYSVXX,
|
||
|
|
$ ZGBSVXX, ZLACPY, LSAMEN
|
||
|
|
LOGICAL LSAMEN
|
||
|
|
|
||
|
|
* .. Intrinsic Functions ..
|
||
|
|
INTRINSIC SQRT, MAX, ABS, DBLE, DIMAG
|
||
|
|
|
||
|
|
* .. Statement Functions ..
|
||
|
|
DOUBLE PRECISION CABS1
|
||
|
|
|
||
|
|
* .. Statement Function Definitions ..
|
||
|
|
CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
|
||
|
|
|
||
|
|
* .. Parameters ..
|
||
|
|
INTEGER NWISE_I, CWISE_I
|
||
|
|
PARAMETER (NWISE_I = 1, CWISE_I = 1)
|
||
|
|
INTEGER BND_I, COND_I
|
||
|
|
PARAMETER (BND_I = 2, COND_I = 3)
|
||
|
|
|
||
|
|
* Create the loop to test out the Hilbert matrices
|
||
|
|
|
||
|
|
FACT = 'E'
|
||
|
|
UPLO = 'U'
|
||
|
|
TRANS = 'N'
|
||
|
|
EQUED = 'N'
|
||
|
|
EPS = DLAMCH('Epsilon')
|
||
|
|
NFAIL = 0
|
||
|
|
N_AUX_TESTS = 0
|
||
|
|
LDA = NMAX
|
||
|
|
LDAB = (NMAX-1)+(NMAX-1)+1
|
||
|
|
LDAFB = 2*(NMAX-1)+(NMAX-1)+1
|
||
|
|
C2 = PATH( 2: 3 )
|
||
|
|
|
||
|
|
* Main loop to test the different Hilbert Matrices.
|
||
|
|
|
||
|
|
printed_guide = .false.
|
||
|
|
|
||
|
|
DO N = 1 , NMAX
|
||
|
|
PARAMS(1) = -1
|
||
|
|
PARAMS(2) = -1
|
||
|
|
|
||
|
|
KL = N-1
|
||
|
|
KU = N-1
|
||
|
|
NRHS = n
|
||
|
|
M = MAX(SQRT(DBLE(N)), 10.0D+0)
|
||
|
|
|
||
|
|
* Generate the Hilbert matrix, its inverse, and the
|
||
|
|
* right hand side, all scaled by the LCM(1,..,2N-1).
|
||
|
|
CALL ZLAHILB(N, N, A, LDA, INVHILB, LDA, B,
|
||
|
|
$ LDA, WORK, INFO, PATH)
|
||
|
|
|
||
|
|
* Copy A into ACOPY.
|
||
|
|
CALL ZLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)
|
||
|
|
|
||
|
|
* Store A in band format for GB tests
|
||
|
|
DO J = 1, N
|
||
|
|
DO I = 1, KL+KU+1
|
||
|
|
AB( I, J ) = (0.0D+0,0.0D+0)
|
||
|
|
END DO
|
||
|
|
END DO
|
||
|
|
DO J = 1, N
|
||
|
|
DO I = MAX( 1, J-KU ), MIN( N, J+KL )
|
||
|
|
AB( KU+1+I-J, J ) = A( I, J )
|
||
|
|
END DO
|
||
|
|
END DO
|
||
|
|
|
||
|
|
* Copy AB into ABCOPY.
|
||
|
|
DO J = 1, N
|
||
|
|
DO I = 1, KL+KU+1
|
||
|
|
ABCOPY( I, J ) = (0.0D+0,0.0D+0)
|
||
|
|
END DO
|
||
|
|
END DO
|
||
|
|
CALL ZLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)
|
||
|
|
|
||
|
|
* Call Z**SVXX with default PARAMS and N_ERR_BND = 3.
|
||
|
|
IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
|
||
|
|
CALL ZSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
|
||
|
|
$ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
|
||
|
|
$ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
|
||
|
|
$ PARAMS, WORK, RWORK, INFO)
|
||
|
|
ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN
|
||
|
|
CALL ZPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
|
||
|
|
$ EQUED, S, B, LDA, X, LDA, ORCOND,
|
||
|
|
$ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
|
||
|
|
$ PARAMS, WORK, RWORK, INFO)
|
||
|
|
ELSE IF ( LSAMEN( 2, C2, 'HE' ) ) THEN
|
||
|
|
CALL ZHESVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
|
||
|
|
$ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
|
||
|
|
$ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
|
||
|
|
$ PARAMS, WORK, RWORK, INFO)
|
||
|
|
ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN
|
||
|
|
CALL ZGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY,
|
||
|
|
$ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B,
|
||
|
|
$ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND,
|
||
|
|
$ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, RWORK,
|
||
|
|
$ INFO)
|
||
|
|
ELSE
|
||
|
|
CALL ZGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA,
|
||
|
|
$ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND,
|
||
|
|
$ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
|
||
|
|
$ PARAMS, WORK, RWORK, INFO)
|
||
|
|
END IF
|
||
|
|
|
||
|
|
N_AUX_TESTS = N_AUX_TESTS + 1
|
||
|
|
IF (ORCOND .LT. EPS) THEN
|
||
|
|
! Either factorization failed or the matrix is flagged, and 1 <=
|
||
|
|
! INFO <= N+1. We don't decide based on rcond anymore.
|
||
|
|
! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
|
||
|
|
! NFAIL = NFAIL + 1
|
||
|
|
! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
|
||
|
|
! END IF
|
||
|
|
ELSE
|
||
|
|
! Either everything succeeded (INFO == 0) or some solution failed
|
||
|
|
! to converge (INFO > N+1).
|
||
|
|
IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN
|
||
|
|
NFAIL = NFAIL + 1
|
||
|
|
WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND
|
||
|
|
END IF
|
||
|
|
END IF
|
||
|
|
|
||
|
|
* Calculating the difference between Z**SVXX's X and the true X.
|
||
|
|
DO I = 1,N
|
||
|
|
DO J =1,NRHS
|
||
|
|
DIFF(I,J) = X(I,J) - INVHILB(I,J)
|
||
|
|
END DO
|
||
|
|
END DO
|
||
|
|
|
||
|
|
* Calculating the RCOND
|
||
|
|
RNORM = 0
|
||
|
|
RINORM = 0
|
||
|
|
IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) .OR.
|
||
|
|
$ LSAMEN( 2, C2, 'HE' ) ) THEN
|
||
|
|
DO I = 1, N
|
||
|
|
SUMR = 0
|
||
|
|
SUMRI = 0
|
||
|
|
DO J = 1, N
|
||
|
|
SUMR = SUMR + S(I) * CABS1(A(I,J)) * S(J)
|
||
|
|
SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (S(J) * S(I))
|
||
|
|
END DO
|
||
|
|
RNORM = MAX(RNORM,SUMR)
|
||
|
|
RINORM = MAX(RINORM,SUMRI)
|
||
|
|
END DO
|
||
|
|
ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
|
||
|
|
$ THEN
|
||
|
|
DO I = 1, N
|
||
|
|
SUMR = 0
|
||
|
|
SUMRI = 0
|
||
|
|
DO J = 1, N
|
||
|
|
SUMR = SUMR + R(I) * CABS1(A(I,J)) * C(J)
|
||
|
|
SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (R(J) * C(I))
|
||
|
|
END DO
|
||
|
|
RNORM = MAX(RNORM,SUMR)
|
||
|
|
RINORM = MAX(RINORM,SUMRI)
|
||
|
|
END DO
|
||
|
|
END IF
|
||
|
|
|
||
|
|
RNORM = RNORM / CABS1(A(1, 1))
|
||
|
|
RCOND = 1.0D+0/(RNORM * RINORM)
|
||
|
|
|
||
|
|
* Calculating the R for normwise rcond.
|
||
|
|
DO I = 1, N
|
||
|
|
RINV(I) = 0.0D+0
|
||
|
|
END DO
|
||
|
|
DO J = 1, N
|
||
|
|
DO I = 1, N
|
||
|
|
RINV(I) = RINV(I) + CABS1(A(I,J))
|
||
|
|
END DO
|
||
|
|
END DO
|
||
|
|
|
||
|
|
* Calculating the Normwise rcond.
|
||
|
|
RINORM = 0.0D+0
|
||
|
|
DO I = 1, N
|
||
|
|
SUMRI = 0.0D+0
|
||
|
|
DO J = 1, N
|
||
|
|
SUMRI = SUMRI + CABS1(INVHILB(I,J) * RINV(J))
|
||
|
|
END DO
|
||
|
|
RINORM = MAX(RINORM, SUMRI)
|
||
|
|
END DO
|
||
|
|
|
||
|
|
! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
|
||
|
|
! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
|
||
|
|
NCOND = CABS1(A(1,1)) / RINORM
|
||
|
|
|
||
|
|
CONDTHRESH = M * EPS
|
||
|
|
ERRTHRESH = M * EPS
|
||
|
|
|
||
|
|
DO K = 1, NRHS
|
||
|
|
NORMT = 0.0D+0
|
||
|
|
NORMDIF = 0.0D+0
|
||
|
|
CWISE_ERR = 0.0D+0
|
||
|
|
DO I = 1, N
|
||
|
|
NORMT = MAX(CABS1(INVHILB(I, K)), NORMT)
|
||
|
|
NORMDIF = MAX(CABS1(X(I,K) - INVHILB(I,K)), NORMDIF)
|
||
|
|
IF (INVHILB(I,K) .NE. 0.0D+0) THEN
|
||
|
|
CWISE_ERR = MAX(CABS1(X(I,K) - INVHILB(I,K))
|
||
|
|
$ /CABS1(INVHILB(I,K)), CWISE_ERR)
|
||
|
|
ELSE IF (X(I, K) .NE. 0.0D+0) THEN
|
||
|
|
CWISE_ERR = DLAMCH('OVERFLOW')
|
||
|
|
END IF
|
||
|
|
END DO
|
||
|
|
IF (NORMT .NE. 0.0D+0) THEN
|
||
|
|
NWISE_ERR = NORMDIF / NORMT
|
||
|
|
ELSE IF (NORMDIF .NE. 0.0D+0) THEN
|
||
|
|
NWISE_ERR = DLAMCH('OVERFLOW')
|
||
|
|
ELSE
|
||
|
|
NWISE_ERR = 0.0D+0
|
||
|
|
ENDIF
|
||
|
|
|
||
|
|
DO I = 1, N
|
||
|
|
RINV(I) = 0.0D+0
|
||
|
|
END DO
|
||
|
|
DO J = 1, N
|
||
|
|
DO I = 1, N
|
||
|
|
RINV(I) = RINV(I) + CABS1(A(I, J) * INVHILB(J, K))
|
||
|
|
END DO
|
||
|
|
END DO
|
||
|
|
RINORM = 0.0D+0
|
||
|
|
DO I = 1, N
|
||
|
|
SUMRI = 0.0D+0
|
||
|
|
DO J = 1, N
|
||
|
|
SUMRI = SUMRI
|
||
|
|
$ + CABS1(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
|
||
|
|
END DO
|
||
|
|
RINORM = MAX(RINORM, SUMRI)
|
||
|
|
END DO
|
||
|
|
! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
|
||
|
|
! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
|
||
|
|
CCOND = CABS1(A(1,1))/RINORM
|
||
|
|
|
||
|
|
! Forward error bound tests
|
||
|
|
NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
|
||
|
|
CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
|
||
|
|
NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
|
||
|
|
CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
|
||
|
|
! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
|
||
|
|
! $ condthresh, ncond.ge.condthresh
|
||
|
|
! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
|
||
|
|
IF (NCOND .GE. CONDTHRESH) THEN
|
||
|
|
NGUAR = 'YES'
|
||
|
|
IF (NWISE_BND .GT. ERRTHRESH) THEN
|
||
|
|
TSTRAT(1) = 1/(2.0D+0*EPS)
|
||
|
|
ELSE
|
||
|
|
IF (NWISE_BND .NE. 0.0D+0) THEN
|
||
|
|
TSTRAT(1) = NWISE_ERR / NWISE_BND
|
||
|
|
ELSE IF (NWISE_ERR .NE. 0.0D+0) THEN
|
||
|
|
TSTRAT(1) = 1/(16.0*EPS)
|
||
|
|
ELSE
|
||
|
|
TSTRAT(1) = 0.0D+0
|
||
|
|
END IF
|
||
|
|
IF (TSTRAT(1) .GT. 1.0D+0) THEN
|
||
|
|
TSTRAT(1) = 1/(4.0D+0*EPS)
|
||
|
|
END IF
|
||
|
|
END IF
|
||
|
|
ELSE
|
||
|
|
NGUAR = 'NO'
|
||
|
|
IF (NWISE_BND .LT. 1.0D+0) THEN
|
||
|
|
TSTRAT(1) = 1/(8.0D+0*EPS)
|
||
|
|
ELSE
|
||
|
|
TSTRAT(1) = 1.0D+0
|
||
|
|
END IF
|
||
|
|
END IF
|
||
|
|
! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
|
||
|
|
! $ condthresh, ccond.ge.condthresh
|
||
|
|
! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
|
||
|
|
IF (CCOND .GE. CONDTHRESH) THEN
|
||
|
|
CGUAR = 'YES'
|
||
|
|
IF (CWISE_BND .GT. ERRTHRESH) THEN
|
||
|
|
TSTRAT(2) = 1/(2.0D+0*EPS)
|
||
|
|
ELSE
|
||
|
|
IF (CWISE_BND .NE. 0.0D+0) THEN
|
||
|
|
TSTRAT(2) = CWISE_ERR / CWISE_BND
|
||
|
|
ELSE IF (CWISE_ERR .NE. 0.0D+0) THEN
|
||
|
|
TSTRAT(2) = 1/(16.0D+0*EPS)
|
||
|
|
ELSE
|
||
|
|
TSTRAT(2) = 0.0D+0
|
||
|
|
END IF
|
||
|
|
IF (TSTRAT(2) .GT. 1.0D+0) TSTRAT(2) = 1/(4.0D+0*EPS)
|
||
|
|
END IF
|
||
|
|
ELSE
|
||
|
|
CGUAR = 'NO'
|
||
|
|
IF (CWISE_BND .LT. 1.0D+0) THEN
|
||
|
|
TSTRAT(2) = 1/(8.0D+0*EPS)
|
||
|
|
ELSE
|
||
|
|
TSTRAT(2) = 1.0D+0
|
||
|
|
END IF
|
||
|
|
END IF
|
||
|
|
|
||
|
|
! Backwards error test
|
||
|
|
TSTRAT(3) = BERR(K)/EPS
|
||
|
|
|
||
|
|
! Condition number tests
|
||
|
|
TSTRAT(4) = RCOND / ORCOND
|
||
|
|
IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0D+0)
|
||
|
|
$ TSTRAT(4) = 1.0D+0 / TSTRAT(4)
|
||
|
|
|
||
|
|
TSTRAT(5) = NCOND / NWISE_RCOND
|
||
|
|
IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0D+0)
|
||
|
|
$ TSTRAT(5) = 1.0D+0 / TSTRAT(5)
|
||
|
|
|
||
|
|
TSTRAT(6) = CCOND / NWISE_RCOND
|
||
|
|
IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0D+0)
|
||
|
|
$ TSTRAT(6) = 1.0D+0 / TSTRAT(6)
|
||
|
|
|
||
|
|
DO I = 1, NTESTS
|
||
|
|
IF (TSTRAT(I) .GT. THRESH) THEN
|
||
|
|
IF (.NOT.PRINTED_GUIDE) THEN
|
||
|
|
WRITE(*,*)
|
||
|
|
WRITE( *, 9996) 1
|
||
|
|
WRITE( *, 9995) 2
|
||
|
|
WRITE( *, 9994) 3
|
||
|
|
WRITE( *, 9993) 4
|
||
|
|
WRITE( *, 9992) 5
|
||
|
|
WRITE( *, 9991) 6
|
||
|
|
WRITE( *, 9990) 7
|
||
|
|
WRITE( *, 9989) 8
|
||
|
|
WRITE(*,*)
|
||
|
|
PRINTED_GUIDE = .TRUE.
|
||
|
|
END IF
|
||
|
|
WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
|
||
|
|
NFAIL = NFAIL + 1
|
||
|
|
END IF
|
||
|
|
END DO
|
||
|
|
END DO
|
||
|
|
|
||
|
|
c$$$ WRITE(*,*)
|
||
|
|
c$$$ WRITE(*,*) 'Normwise Error Bounds'
|
||
|
|
c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
|
||
|
|
c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
|
||
|
|
c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
|
||
|
|
c$$$ WRITE(*,*)
|
||
|
|
c$$$ WRITE(*,*) 'Componentwise Error Bounds'
|
||
|
|
c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
|
||
|
|
c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
|
||
|
|
c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
|
||
|
|
c$$$ print *, 'Info: ', info
|
||
|
|
c$$$ WRITE(*,*)
|
||
|
|
* WRITE(*,*) 'TSTRAT: ',TSTRAT
|
||
|
|
|
||
|
|
END DO
|
||
|
|
|
||
|
|
WRITE(*,*)
|
||
|
|
IF( NFAIL .GT. 0 ) THEN
|
||
|
|
WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
|
||
|
|
ELSE
|
||
|
|
WRITE(*,9997) C2
|
||
|
|
END IF
|
||
|
|
9999 FORMAT( ' Z', A2, 'SVXX: N =', I2, ', RHS = ', I2,
|
||
|
|
$ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
|
||
|
|
$ ' test(',I1,') =', G12.5 )
|
||
|
|
9998 FORMAT( ' Z', A2, 'SVXX: ', I6, ' out of ', I6,
|
||
|
|
$ ' tests failed to pass the threshold' )
|
||
|
|
9997 FORMAT( ' Z', A2, 'SVXX passed the tests of error bounds' )
|
||
|
|
* Test ratios.
|
||
|
|
9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
|
||
|
|
$ 'Guaranteed case: if norm ( abs( Xc - Xt )',
|
||
|
|
$ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
|
||
|
|
$ / 5X,
|
||
|
|
$ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
|
||
|
|
9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
|
||
|
|
9994 FORMAT( 3X, I2, ': Backwards error' )
|
||
|
|
9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
|
||
|
|
9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
|
||
|
|
9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
|
||
|
|
9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
|
||
|
|
9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )
|
||
|
|
|
||
|
|
8000 FORMAT( ' Z', A2, 'SVXX: N =', I2, ', INFO = ', I3,
|
||
|
|
$ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )
|
||
|
|
*
|
||
|
|
* End of ZEBCHVXX
|
||
|
|
*
|
||
|
|
END
|