debchvxx.f

Go to the documentation of this file.

*> \brief \b DEBCHVXX
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*     SUBROUTINE DEBCHVXX( THRESH, PATH )
*
*     .. Scalar Arguments ..
*      DOUBLE PRECISION  THRESH
*      CHARACTER*3       PATH
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>  DEBCHVXX will run D**SVXX on a series of Hilbert matrices and then
*>  compare the error bounds returned by D**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 D**SVXX
*>          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. 
*
*> \date November 2011
*
*> \ingroup double_lin
*
*  =====================================================================
      SUBROUTINE debchvxx( 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

*     .. Local Arrays ..
      DOUBLE PRECISION   TSTRAT(ntests), RINV(nmax), PARAMS(nparams),
     $                   s(nmax),r(nmax),c(nmax), diff(nmax, nmax),
     $                   errbnd_n(nmax*3), errbnd_c(nmax*3),
     $                   a(nmax,nmax),invhilb(nmax,nmax),x(nmax,nmax),
     $                   ab( (nmax-1)+(nmax-1)+1, nmax ),
     $                   abcopy( (nmax-1)+(nmax-1)+1, nmax ),
     $                   afb( 2*(nmax-1)+(nmax-1)+1, nmax ),
     $                   work(nmax*3*5), af(nmax, nmax),b(nmax, nmax),
     $                   acopy(nmax, nmax)
      INTEGER            IPIV(nmax), IWORK(3*nmax)

*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH

*     .. External Subroutines ..
      EXTERNAL           dlahilb, dgesvxx, dposvxx, dsysvxx,
     $                   dgbsvxx, dlacpy, lsamen
      LOGICAL            LSAMEN

*     .. Intrinsic Functions ..
      INTRINSIC          sqrt, max, abs, dble

*     .. 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 dlahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)

*        Copy A into ACOPY.
         CALL dlacpy('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
            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
            END DO
         END DO
         CALL dlacpy('ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)

*        Call D**SVXX with default PARAMS and N_ERR_BND = 3.
         IF ( lsamen( 2, c2, 'SY' ) ) THEN
            CALL dsysvxx(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, iwork, info)
         ELSE IF ( lsamen( 2, c2, 'PO' ) ) THEN
            CALL dposvxx(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, iwork, info)
         ELSE IF ( lsamen( 2, c2, 'GB' ) ) THEN
            CALL dgbsvxx(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, iwork,
     $           info)
         ELSE
            CALL dgesvxx(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, iwork, 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 D**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.0d+0
         rinorm = 0.0d+0
         IF ( lsamen( 2, c2, 'PO' ) .OR. lsamen( 2, c2, 'SY' ) ) THEN
            DO i = 1, n
               sumr = 0.0d+0
               sumri = 0.0d+0
               DO j = 1, n
                  sumr = sumr + s(i) * abs(a(i,j)) * s(j)
                  sumri = sumri + abs(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.0d+0
               sumri = 0.0d+0
               DO j = 1, n
                  sumr = sumr + r(i) * abs(a(i,j)) * c(j)
                  sumri = sumri + abs(invhilb(i, j)) / (r(j) * c(i))
               END DO
               rnorm = max(rnorm,sumr)
               rinorm = max(rinorm,sumri)
            END DO
         END IF

         rnorm = rnorm / abs(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) + abs(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 + abs(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 = abs(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(abs(invhilb(i, k)), normt)
               normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
               IF (invhilb(i,k) .NE. 0.0d+0) THEN
                  cwise_err = max(abs(x(i,k) - invhilb(i,k))
     $                            /abs(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) + abs(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
     $                 + abs(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 = abs(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( ' D', a2, 'SVXX: N =', i2, ', RHS = ', i2,
     $     ', NWISE GUAR. = ', a, ', CWISE GUAR. = ', a,
     $     ' test(',i1,') =', g12.5 )
 9998 FORMAT( ' D', a2, 'SVXX: ', i6, ' out of ', i6,
     $     ' tests failed to pass the threshold' )
 9997 FORMAT( ' D', 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( ' D', a2, 'SVXX: N =', i2, ', INFO = ', i3,
     $     ', ORCOND = ', g12.5, ', real RCOND = ', g12.5 )

      END