!   E02DHF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.
    MODULE e02dhfe_mod

!      E02DHF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: nin = 5, nout = 6
    CONTAINS
       SUBROUTINE print_spline(ngx,gridx,ngy,gridy,z,zder)

!         Print spline function and spline derivative evaluation

!         .. Use Statements ..
          USE nag_library, ONLY : x04cbf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          INTEGER, PARAMETER                  :: indent = 0, ncols = 80
          CHARACTER (1), PARAMETER            :: chlabel =   'C',  diag =      &
                                                 'N',  matrix =   'G'
          CHARACTER (4), PARAMETER            :: form = 'F8.3'
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: ngx, ngy
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: gridx(ngx), gridy(ngy),       &
                                                 z(ngx*ngy), zder(ngx*ngy)
!         .. Local Scalars ..
          INTEGER                             :: i, ifail
          CHARACTER (48)                      :: title
!         .. Local Arrays ..
          CHARACTER (10), ALLOCATABLE         :: clabs(:), rlabs(:)
!         .. Executable Statements ..
!         Allocate for row and column label
          ALLOCATE (clabs(ngx),rlabs(ngy))
!         Generate column and row labels to print the results with.
          DO i = 1, ngx
             WRITE (clabs(i),99999) gridx(i)
          END DO
          DO i = 1, ngy
             WRITE (rlabs(i),99999) gridy(i)
          END DO

!         Print the spline evaluations.
          title = 'Spline evaluated on X-Y grid (X across, Y down):'
          WRITE (nout,*)
          FLUSH (nout)
          ifail = 0
          CALL x04cbf(matrix,diag,ngy,ngx,z,ngy,form,title,chlabel,rlabs, &
             chlabel,clabs,ncols,indent,ifail)

!         Print the spline derivative evaluations.
          title = 'Spline derivative evaluated on X-Y grid:'
          WRITE (nout,*)
          FLUSH (nout)
          ifail = 0
          CALL x04cbf(matrix,diag,ngy,ngx,zder,ngy,form,title(1:40),chlabel, &
             rlabs,chlabel,clabs,ncols,indent,ifail)

          DEALLOCATE (clabs,rlabs)

99999     FORMAT (F5.2)
       END SUBROUTINE print_spline
    END MODULE e02dhfe_mod
    PROGRAM e02dhfe

!      E02DHF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : e02dcf, e02dhf, nag_wp
       USE e02dhfe_mod, ONLY : nin, nout, print_spline
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: delta, fp, s, xhi, xlo, yhi, ylo
       INTEGER                             :: i, ifail, liwrk, lwrk, mx, my,   &
                                              nc, ngx, ngy, nux, nuy, nx,      &
                                              nxest, ny, nyest
       CHARACTER (1)                       :: start
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: c(:), f(:), gridx(:), gridy(:),  &
                                              lamda(:), mu(:), wrk(:), x(:),   &
                                              y(:), z(:), zder(:)
       INTEGER, ALLOCATABLE                :: iwrk(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              max, real
!      .. Executable Statements ..
       WRITE (nout,*) 'E02DHF Example Program Results'

!      Skip heading in data file
       READ (nin,*)

!      Input the number of X, Y co-ordinates MX, MY.

       READ (nin,*) mx, my
       nxest = mx + 4
       nyest = my + 4
       nc = (nxest-4)*(nyest-4)

!      Allocations for spline fit
       ALLOCATE (lamda(nxest),mu(nyest),c(nc))

!      Allocations for e02dcf only
       lwrk = 4*(mx+my) + 11*(nxest+nyest) + nxest*my + max(my,nxest) + 54
       liwrk = 3 + mx + my + nxest + nyest
       ALLOCATE (x(mx),y(my),f(mx*my),wrk(lwrk),iwrk(liwrk))

       READ (nin,*) x(1:mx)
       READ (nin,*) y(1:my)

!      Input the MX*MY function values F at grid points nad smoothing factor.

       READ (nin,*) f(1:mx*my)
       READ (nin,*) s

!      Determine the spline approximation.
       start = 'C'
       ifail = 0
       CALL e02dcf(start,mx,x,my,y,f,s,nxest,nyest,nx,lamda,ny,mu,c,fp,wrk, &
          lwrk,iwrk,liwrk,ifail)

       DEALLOCATE (x,y,f,wrk,iwrk)

       WRITE (nout,*)
       WRITE (nout,99999) 'Spline fit used smoothing factor S =', s, '.'
       WRITE (nout,99998) 'Number of knots in each direction  =', nx, ny
       WRITE (nout,*)
       WRITE (nout,99999) 'Sum of squared residuals           =', fp, '.'

!      Spline and its derivative to be evaluated on rectangular grid with
!      ngx*ngy points on the domain [xlo,xhi]x[ylo,yhi].

       READ (nin,*) ngx, xlo, xhi
       READ (nin,*) ngy, ylo, yhi

!      Allocations for e02dhf (spline evaluation).
       ALLOCATE (gridx(ngx),gridy(ngy),z(ngx*ngy),zder(ngx*ngy))

       delta = (xhi-xlo)/real(ngx-1,kind=nag_wp)
       gridx(1) = xlo
       DO i = 2, ngx - 1
          gridx(i) = gridx(i-1) + delta
       END DO
       gridx(ngx) = xhi

       delta = (yhi-ylo)/real(ngy-1,kind=nag_wp)
       gridy(1) = ylo
       DO i = 2, ngy - 1
          gridy(i) = gridy(i-1) + delta
       END DO
       gridy(ngy) = yhi

!      Evaluate spline (nux=nuy=0)
       nux = 0
       nuy = 0
       ifail = 0
       CALL e02dhf(ngx,ngy,nx,ny,gridx,gridy,lamda,mu,c,nux,nuy,z,ifail)
!      Evaluate spline partial derivative of order (nux,nuy)
       READ (nin,*) nux, nuy
       WRITE (nout,*)
       WRITE (nout,99998) 'Derivative of spline has order nux, nuy =', nux, &
          nuy
       ifail = 0
       CALL e02dhf(ngx,ngy,nx,ny,gridx,gridy,lamda,mu,c,nux,nuy,zder,ifail)

!      Print tabulated spline and derivative evaluations.
       CALL print_spline(ngx,gridx,ngy,gridy,z,zder)

99999  FORMAT (1X,A,1P,E13.4,A)
99998  FORMAT (1X,A,I5,',',I5,'.')
    END PROGRAM e02dhfe