PROGRAM e01dafe

!      E01DAF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : e01daf, e02dff, nag_wp, x04cbf
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: indent = 0, ncols = 80, nin = 5,     &
                                          nout = 6
       CHARACTER (1), PARAMETER        :: chlabel =   'C',  diag =   'N',      &
                                          matrix =   'G'
       CHARACTER (4), PARAMETER        :: form = 'F8.3'
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: step, xhi, xlo, yhi, ylo
       INTEGER                         :: i, ifail, j, liwrk, lwrk, mx, my,    &
                                          nx, ny, px, py
       CHARACTER (54)                  :: title
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: c(:), f(:,:), fg(:), lamda(:),       &
                                          mu(:), tx(:), ty(:), wrk(:), x(:),   &
                                          y(:)
       INTEGER, ALLOCATABLE            :: iwrk(:)
       CHARACTER (10), ALLOCATABLE     :: clabs(:), rlabs(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          min, real
!      .. Executable Statements ..
       WRITE (nout,*) 'E01DAF Example Program Results'

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

!      Read the number of X points, MX, and the values of the
!      X co-ordinates.

       READ (nin,*) mx
       ALLOCATE (x(mx),lamda(mx+4))

       READ (nin,*) x(1:mx)

!      Read the number of Y points, MY, and the values of the
!      Y co-ordinates.

       READ (nin,*) my
       ALLOCATE (y(my),mu(my+4),c(mx*my),f(my,mx),wrk((mx+6)*(my+6)))

       READ (nin,*) y(1:my)

!      Read the function values at the grid points.

       DO j = 1, my
          READ (nin,*) f(j,1:mx)
       END DO

!      Generate the (X,Y,F) interpolating bicubic B-spline.

       ifail = 0
       CALL e01daf(mx,my,x,y,f,px,py,lamda,mu,c,wrk,ifail)

!      Print the knot sets, LAMDA and MU.

       WRITE (nout,*)
       WRITE (nout,*) '               I   Knot LAMDA(I)      J     Knot MU(J)'
       WRITE (nout,99997) (j,lamda(j),j,mu(j),j=4,min(px,py)-3)
       IF (px>py) THEN
          WRITE (nout,99997) (j,lamda(j),j=py-2,px-3)
       ELSE IF (px<py) THEN
          WRITE (nout,99996) (j,mu(j),j=px-2,py-3)
       END IF

!      Print the spline coefficients.

       WRITE (nout,*)
       WRITE (nout,*) 'The B-Spline coefficients:'
       WRITE (nout,99999) (c(i),i=1,mx*my)
       WRITE (nout,*)
       FLUSH (nout)

!      Evaluate the spline on a regular rectangular grid at nx*ny points
!      over the domain [xlo,xhi] x [ylo,yhi].

       DEALLOCATE (wrk)

       READ (nin,*) nx, xlo, xhi
       READ (nin,*) ny, ylo, yhi
       lwrk = min(4*nx+px,4*ny+py)

       IF (4*nx+px>4*ny+py) THEN
          liwrk = ny + py - 4
       ELSE
          liwrk = nx + px - 4
       END IF

       ALLOCATE (tx(nx),ty(ny),fg(nx*ny),wrk(lwrk),iwrk(liwrk))

!      Generate nx/ny equispaced x/y co-ordinates.

       step = (xhi-xlo)/real(nx-1,kind=nag_wp)
       tx(1) = xlo
       DO i = 2, nx - 1
          tx(i) = tx(i-1) + step
       END DO
       tx(nx) = xhi

       step = (yhi-ylo)/real(ny-1,kind=nag_wp)
       ty(1) = ylo
       DO i = 2, ny - 1
          ty(i) = ty(i-1) + step
       END DO
       ty(ny) = yhi

!      Evaluate the spline.
       ifail = 0
       CALL e02dff(nx,ny,px,py,tx,ty,lamda,mu,c,fg,wrk,lwrk,iwrk,liwrk,ifail)

!      Generate row and column labels and title for printing results.
       ALLOCATE (clabs(nx),rlabs(ny))
       DO i = 1, nx
          WRITE (clabs(i),99998) tx(i)
       END DO
       DO i = 1, ny
          WRITE (rlabs(i),99998) ty(i)
          FLUSH (nout)
       END DO
       title = 'Spline evaluated on a regular mesh (X across, Y down):'

!      Print the results.
       ifail = 0
       CALL x04cbf(matrix,diag,ny,nx,fg,ny,form,title,chlabel,rlabs,chlabel, &
          clabs,ncols,indent,ifail)

99999  FORMAT (1X,8F9.4)
99998  FORMAT (F5.2)
99997  FORMAT (1X,I16,F12.4,I11,F12.4)
99996  FORMAT (1X,I39,F12.4)
    END PROGRAM e01dafe