!   E02DDF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.
    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
!      .. 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,7F9.2)
       END SUBROUTINE cprint
    END MODULE e02ddfe_mod
    PROGRAM e02ddfe

!      E02DDF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : e02ddf, e02dff, nag_wp
       USE e02ddfe_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, 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 Functions ..
       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