!   E02DCF Example Program Text
!   Mark 28.4 Release. NAG Copyright 2022.
    Module e02dcfe_mod

!     E02DCF 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                           :: cprint
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
    Contains
      Subroutine cprint(c,ny,nx,nout)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: nout, nx, ny
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: c(ny-4,nx-4)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*) 'The B-spline coefficients:'
        Write (nout,*)

        Do i = 1, ny - 4
          Write (nout,99999) c(i,1:(nx-4))
        End Do

        Return

99999   Format (1X,8F9.4)
      End Subroutine cprint
    End Module e02dcfe_mod
    Program e02dcfe

!     E02DCF Example Main Program

!     .. Use Statements ..
      Use e02dcfe_mod, Only: cprint, nin, nout
      Use nag_library, Only: e02dcf, e02dff, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: delta, fp, s, xhi, xlo, yhi, ylo
      Integer                          :: i, ifail, j, liwrk, lwrk, mx, my,    &
                                          npx, npy, nx, nxest, ny, nyest
      Character (1)                    :: start
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: c(:), f(:), fg(:), lamda(:), mu(:),  &
                                          px(:), py(:), wrk(:), x(:), y(:)
      Integer, Allocatable             :: iwrk(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, min, real
!     .. Executable Statements ..
      Write (nout,*) 'E02DCF 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
      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),lamda(nxest),mu(nyest),wrk(lwrk),         &
        iwrk(liwrk),c((nxest-4)*(nyest-4)))

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

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

      Read (nin,*) f(1:mx*my)

      start = 'C'

      Read (nin,*) s

!     Determine the spline approximation.

      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 (wrk,iwrk)

      Write (nout,*)
      Write (nout,99999) 'Calling with smoothing factor S =', s, ': NX =', nx, &
        ', NY =', ny, '.'
      Write (nout,*)
      Write (nout,*) '               I   Knot LAMDA(I)      J     Knot MU(J)'
      Write (nout,*)

      Do j = 4, max(nx,ny) - 3

        If (j<=nx-3 .And. j<=ny-3) Then
          Write (nout,99997) j, lamda(j), j, mu(j)
        Else If (j<=nx-3) Then
          Write (nout,99997) j, lamda(j)
        Else If (j<=ny-3) Then
          Write (nout,99996) j, mu(j)
        End If

      End Do

      Call cprint(c,ny,nx,nout)

      Write (nout,*)
      Write (nout,99998) 'Sum of squared residuals FP =', fp

      If (fp==0.0E+0_nag_wp) Then
        Write (nout,*) '(The spline is an interpolating spline)'
      Else If (nx==8 .And. ny==8) Then
        Write (nout,*) '(The spline is the least squares bi-cubic polynomial)'
      End If

!     Evaluate the spline on a rectangular grid at NPX*NPY points
!     over the domain (XLO to XHI) x (YLO to YHI).

      Read (nin,*) npx, xlo, xhi
      Read (nin,*) npy, ylo, yhi

      lwrk = min(4*npx+nx,4*npy+ny)

      If (4*npx+nx>4*npy+ny) Then
        liwrk = npy + ny - 4
      Else
        liwrk = npx + nx - 4
      End If

      Allocate (px(npx),py(npy),fg(npx*npy),wrk(lwrk),iwrk(liwrk))

      delta = (xhi-xlo)/real(npx-1,kind=nag_wp)

      Do i = 1, npx
        px(i) = min(xlo+real(i-1,kind=nag_wp)*delta,xhi)
      End Do

      delta = (yhi-ylo)/real(npy-1,kind=nag_wp)

      Do i = 1, npy
        py(i) = min(ylo+real(i-1,kind=nag_wp)*delta,yhi)
      End Do

      ifail = 0
      Call e02dff(npx,npy,nx,ny,px,py,lamda,mu,c,fg,wrk,lwrk,iwrk,liwrk,ifail)

      Write (nout,*)
      Write (nout,*) 'Values of computed spline:'
      Write (nout,*)
      Write (nout,99995) '          X', (px(i),i=1,npx)
      Write (nout,*) '     Y'

      Do i = npy, 1, -1
        Write (nout,99994) py(i), (fg(npy*(j-1)+i),j=1,npx)
      End Do

99999 Format (1X,A,1P,E13.4,A,I5,A,I5,A)
99998 Format (1X,A,1P,E13.4)
99997 Format (1X,I16,F12.4,I11,F12.4)
99996 Format (1X,I39,F12.4)
99995 Format (1X,A,7F8.2)
99994 Format (1X,F8.2,3X,7F8.2)
    End Program e02dcfe