zebchvxx()

subroutine zebchvxx	(	double precision	thresh,
		character*3	path
	)

ZEBCHVXX

Purpose:

  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.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 95 of file zebchvxx.f.

      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 )',
     $     .LE.' / norm ( Xt )  ERRBND( *, nwise_i, bnd_i ), then',
     $     / 5x,
     $     .LE.'ERRBND( *, nwise_i, bnd_i )  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
*

Here is the call graph for this function:

Here is the caller graph for this function: