Program d04bafe

!     Mark 24 Release. NAG Copyright 2012.

!     .. 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