!   E02DHF Example Program Text
!   Mark 28.5 Release. NAG Copyright 2022.
    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
!     .. Accessibility Statements ..
      Private
      Public                           :: print_spline
!     .. Parameters ..
      Integer, Parameter, Public       :: 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
!       .. 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 e02dhfe_mod, Only: nin, nout, print_spline
      Use nag_library, Only: e02dcf, e02dhf, nag_wp
!     .. 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 Procedures ..
      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 and 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