```    Program e01dafe

!     E01DAF Example Program Text

!     Mark 26.2 Release. NAG Copyright 2017.

!     .. 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 Procedures ..
Intrinsic                        :: min, real
!     .. Executable Statements ..
Write (nout,*) 'E01DAF Example Program Results'

!     Skip heading in data file

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

Allocate (x(mx),lamda(mx+4))

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

Allocate (y(my),mu(my+4),c(mx*my),f(my,mx),wrk((mx+6)*(my+6)))

!     Read the function values at the grid points.

Do j = 1, my
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
```