PROGRAM d04bafe

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : d04baf, d04bbf, nag_wp, s14aef, x04cbf
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER   :: x_0 = 0.05_nag_wp
       INTEGER, PARAMETER              :: indent = 0, ncols = 80, nout = 6,    &
                                          n_der_comp = 3, n_display = 3,       &
                                          n_hbase = 4, zeroth = 0
       CHARACTER (1), PARAMETER        :: chlabel =   'C',  diag =   'N',      &
                                          form =   '',  matrix =   'G',        &
                                          nolabel =   'N'
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: hbase
       INTEGER                         :: ifail, j, k
       CHARACTER (50)                  :: title
!      .. Local Arrays ..
       REAL (KIND=nag_wp)              :: actder(n_display), der(14),          &
                                          der_comp(n_hbase,n_der_comp,14),     &
                                          erest(14), fval(21), xval(21)
       CHARACTER (10)                  :: clabs(n_der_comp), rlabs(1)
!      .. Executable Statements ..
       WRITE (nout,*) 'D04BAF Example Program Results'
       WRITE (nout,*)
       WRITE (nout,*) ' Find the derivatives of the polygamma (psi) function'
       WRITE (nout,*) ' using function values generated by S14AEF.'
       WRITE (nout,*)
       WRITE (nout,*) &
          ' Demonstrate the effect of successively reducing HBASE.'
       WRITE (nout,*)
!      Select an initial separation distance HBASE. 
       hbase = 0.0025_nag_wp

!      Compute the actual derivatives at target location x_0 using s14aef for
!      comparison.
       DO j = 1, n_display
          ifail = 0
          actder(j) = s14aef(x_0,j,ifail)
       END DO

!      Attempt N_HBASE approximations, reducing HBASE by factor 0.1 each time.
       DO j = 1, n_hbase

!         Generate the abscissa XVAL using D04BBF
          CALL d04bbf(x_0,hbase,xval)

!         Calculate the corresponding objective function values.
          DO k = 1, 21
             ifail = 0
             fval(k) = s14aef(xval(k),zeroth,ifail)
          END DO

!         Call D04BAF to calculate the derivative estimates
          ifail = 0
          CALL d04baf(xval,fval,der,erest,ifail)

!         Store results in DER_COMP
          der_comp(j,1,1:14) = hbase
          der_comp(j,2,1:14) = der(1:14)
          der_comp(j,3,1:14) = erest(1:14)

!         Decrease hbase for next loop
          hbase = hbase*0.1_nag_wp
       END DO

!      Display Results for first N_DISPLAY derivatives

       DO j = 1, n_display
          WRITE (nout,99996) j, actder(j)
          WRITE (clabs(1),99997) 'hbase     '
          WRITE (clabs(2),99998) 'DER', j, '   '
          WRITE (clabs(3),99999) 'EREST', j
          WRITE (title,99999) &
             ' Derivative and error estimates for derivative ', j
          FLUSH (nout)

!         Use X04CBF to display the matrix
          ifail = 0
          CALL x04cbf(matrix,diag,n_hbase,n_der_comp,der_comp(1,1,j),n_hbase, &
             form,title,nolabel,rlabs,chlabel,clabs,ncols,indent,ifail)
          WRITE (nout,*)
       END DO

99999  FORMAT (A,'(',I1,')')
99998  FORMAT (A,'(',I1,')',A)
99997  FORMAT (A)
99996  FORMAT (1X,' Derivative (',I1,') calculated using S14AEF :',1X,ES11.4)
    END PROGRAM d04bafe