Program e01dafe

!     E01DAF Example Program Text

!     Mark 29.2 Release. NAG Copyright 2023.

!     .. 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 (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