```/* nag_2d_spline_interpolant (e01dac) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.2, 2017.
*
*/

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage01.h>
#include <nage02.h>

#define F(I, J)  f[my*(I)+(J)]
#define FG(I, J) fg[npy*(I)+(J)]
#define C(I, J)  spline.c[my*(I)+(J)]
int main(void)
{
Integer exit_status = 0, i, j, mx, my, npx, npy;
NagError fail;
Nag_2dSpline spline;
double *f = 0, *fg = 0, step, *tx = 0, *ty = 0, *x = 0, xhi, xlo;
double *y = 0, yhi, ylo;

INIT_FAIL(fail);

/* Initialize spline */
spline.lamda = 0;
spline.mu = 0;
spline.c = 0;

printf("nag_2d_spline_interpolant (e01dac) Example Program Results\n");
scanf("%*[^\n]"); /* Skip heading in data file */
/* Read the number of x points, mx, and the values of the
* x co-ordinates.
*/
scanf("%" NAG_IFMT "%" NAG_IFMT "", &mx, &my);
if (mx >= 4 && my >= 4) {
if (!(f = NAG_ALLOC(mx * my, double)) ||
!(x = NAG_ALLOC(mx, double)) || !(y = NAG_ALLOC(my, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else {
printf("Invalid mx or my.\n");
exit_status = 1;
return exit_status;
}
for (i = 0; i < mx; i++)
scanf("%lf", &x[i]);
/* Read the number of y points, my, and the values of the
* y co-ordinates.
*/
for (i = 0; i < my; i++)
scanf("%lf", &y[i]);
/* Read the function values at the grid points. */
for (j = 0; j < my; j++)
for (i = 0; i < mx; i++)
scanf("%lf", &F(i, j));
/* Generate the (x,y,f) interpolating bicubic B-spline. */
/* nag_2d_spline_interpolant (e01dac).
* Interpolating function, bicubic spline interpolant, two
* variables
*/
nag_2d_spline_interpolant(mx, my, x, y, f, &spline, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_interpolant (e01dac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* Print the knot sets, lamda and mu. */
printf("Distinct knots in x direction located  at\n");
for (j = 3; j < spline.nx - 3; j++)
printf("%12.4f%s", spline.lamda[j],
((j - 3) % 5 == 4 || j == spline.nx - 4) ? "\n" : " ");
printf("\nDistinct knots in y direction located  at\n");
for (j = 3; j < spline.ny - 3; j++)
printf("%12.4f%s", spline.mu[j],
((j - 3) % 5 == 4 || j == spline.ny - 4) ? "\n" : " ");
/* Print the spline coefficients. */
printf("\nThe B-Spline coefficients:\n");
for (i = 0; i < mx; i++) {
for (j = 0; j < my; j++)
printf("%9.4f", C(i, j));
printf("\n");
}

/* Evaluate the spline on a regular rectangular grid at npx*npy
* points over the domain (xlo to xhi) x (ylo to yhi).
*/
scanf("%" NAG_IFMT "%lf%lf", &npx, &xlo, &xhi);
scanf("%" NAG_IFMT "%lf%lf", &npy, &ylo, &yhi);
if (npx >= 1 && npy >= 1) {
if (!(fg = NAG_ALLOC(npx * npy, double)) ||
!(tx = NAG_ALLOC(npx, double)) || !(ty = NAG_ALLOC(npy, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else {
printf("Invalid npx or npy.\n");
exit_status = 1;
return exit_status;
}
step = (xhi - xlo) / (double) (npx - 1);
printf("\nSpline evaluated on a regular mesh "
"     (x across, y down): \n     ");
/* Generate nx equispaced x co-ordinates. */
for (i = 0; i < npx; i++) {
tx[i] = MIN(xlo + i * step, xhi);
printf("   %5.2f ", tx[i]);
}
step = (yhi - ylo) / (npy - 1);
for (i = 0; i < npy; i++)
ty[i] = MIN(ylo + i * step, yhi);

/* Evaluate the spline. */
/* nag_2d_spline_eval_rect (e02dfc).
* Evaluation of bicubic spline, at a mesh of points
*/
nag_2d_spline_eval_rect(npx, npy, tx, ty, fg, &spline, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_2d_spline_eval_rect (e02dfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* Print the results. */
printf("\n");
for (j = 0; j < npy; j++) {
printf("%5.2f", ty[j]);
for (i = 0; i < npx; i++)
printf("%8.3f ", FG(i, j));
printf("\n");
}
/* Free memory allocated by nag_2d_spline_interpolant (e01dac) */
END:
NAG_FREE(spline.lamda);
NAG_FREE(spline.mu);
NAG_FREE(spline.c);
NAG_FREE(f);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(fg);
NAG_FREE(tx);
NAG_FREE(ty);
return exit_status;
}
```