!   E02DDF Example Program Text
!   Mark 28.6 Release. NAG Copyright 2022.
    Module e02ddfe_mod

!     E02DDF 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,7F9.2)
      End Subroutine cprint
    End Module e02ddfe_mod
    Program e02ddfe

!     E02DDF Example Main Program

!     .. Use Statements ..
      Use e02ddfe_mod, Only: cprint, nin, nout
      Use nag_library, Only: e02ddf, 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, m, npx,    &
                                          npy, nx, nxest, ny, nyest, rank, u,  &
                                          v, ww
      Character (1)                    :: start
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: c(:), f(:), fg(:), lamda(:), mu(:),  &
                                          px(:), py(:), w(:), wrk(:), x(:),    &
                                          y(:)
      Integer, Allocatable             :: iwrk(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, min, real
!     .. Executable Statements ..
      Write (nout,*) 'E02DDF Example Program Results'

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

      Read (nin,*) m
      nxest = m
      nyest = nxest
      liwrk = m + 2*(nxest-7)*(nyest-7)
      u = nxest - 4
      v = nyest - 4
      ww = max(u,v)
      lwrk = (7*u*v+25*ww)*(ww+1) + 2*(u+v+4*m) + 23*ww + 56
      Allocate (x(m),y(m),f(m),w(m),lamda(nxest),mu(nyest),c((nxest-4)*(nyest- &
        4)),iwrk(liwrk),wrk(lwrk))

!     Input the data-points and the weights.

      Do i = 1, m
        Read (nin,*) x(i), y(i), f(i), w(i)
      End Do

      start = 'C'

      Read (nin,*) s

!     Determine the spline approximation.

      ifail = 0
      Call e02ddf(start,m,x,y,f,w,s,nxest,nyest,nx,lamda,ny,mu,c,fp,rank,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,99998) 'rank deficiency =', (nx-4)*(ny-4) - rank
      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,99996) j, lamda(j), j, mu(j)
        Else If (j<=nx-3) Then
          Write (nout,99996) j, lamda(j)
        Else If (j<=ny-3) Then
          Write (nout,99995) j, mu(j)
        End If

      End Do

      Call cprint(c,ny,nx,nout)

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

      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

      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,99994) '          X', (px(i),i=1,npx)
      Write (nout,*) '     Y'

      Do i = npy, 1, -1
        Write (nout,99993) 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,I5)
99997 Format (1X,A,1P,E13.4,A)
99996 Format (1X,I16,F12.4,I11,F12.4)
99995 Format (1X,I39,F12.4)
99994 Format (1X,A,7F8.2)
99993 Format (1X,F8.2,3X,7F8.2)
    End Program e02ddfe