!   E02DCF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.
    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
!      .. Parameters ..
       INTEGER, PARAMETER                  :: nin = 5, nout = 6
    CONTAINS
       SUBROUTINE cprint(c,ny,nx,nout)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 nag_library, ONLY : e02dcf, e02dff, nag_wp
       USE e02dcfe_mod, ONLY : cprint, nin, nout
!      .. 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 Functions ..
       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