/* nag_fit_dim2_spline_grid (e02dcc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 *
 */

#include <nag.h>
#include <stdio.h>

int main(void) {
  Integer exit_status = 0, i, j, mx, my, npx, npy, nx, ny;
  Nag_2dSpline spline;
  Nag_Comm warmstartinf;
  Nag_Start start;
  double delta, *f = 0, *fg = 0, fp, *px = 0, *py = 0, s, *x = 0, xhi, xlo,
                *y = 0, yhi;
  double ylo;
  NagError fail;

  INIT_FAIL(fail);

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

  warmstartinf.nag_w = 0;
  warmstartinf.nag_iw = 0;

  printf("nag_fit_dim2_spline_grid (e02dcc) Example Program Results\n");
  scanf("%*[^\n]"); /* Skip heading in data file */
  /* Input the number of x, y co-ordinates mx, my. */
  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;
  }
  /* Input the x co-ordinates followed by the y co-ordinates. */
  for (i = 0; i < mx; i++)
    scanf("%lf", &x[i]);
  for (i = 0; i < my; i++)
    scanf("%lf", &y[i]);
  /* Input the mx*my function values f at the grid points. */
  for (i = 0; i < mx * my; i++)
    scanf("%lf", &f[i]);
  start = Nag_Cold;
  scanf("%lf", &s);
  /* Determine the spline approximation. */

  /* nag_fit_dim2_spline_grid (e02dcc).
   * Least squares bicubic spline fit with automatic knot
   * placement, two variables (rectangular grid)
   */
  nag_fit_dim2_spline_grid(start, mx, x, my, y, f, s, mx + 4, my + 4, &fp,
                           &warmstartinf, &spline, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_dim2_spline_grid (e02dcc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  nx = spline.nx;
  ny = spline.ny;

  printf("\nCalling with smoothing factor s = %13.4e:"
         " spline.nx = %2" NAG_IFMT ", spline.ny = %2" NAG_IFMT ".\n\n",
         s, nx, ny);
  /* 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" : " ");
  printf("\nThe B-spline coefficients:\n\n");
  for (i = 0; i < ny - 4; i++) {
    for (j = 0; j < nx - 4; j++)
      printf("%9.4f", spline.c[i + j * (ny - 4)]);
    printf("\n");
  }
  printf("\nSum of squared residuals fp = %13.4e\n", fp);
  if (fp == 0.0)
    printf("\nThe spline is an interpolating spline\n");
  else if (nx == 8 && ny == 8)
    printf("\nThe spline is the least squares bi-cubic polynomial\n");

    /* Evaluate the spline on a 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)) ||
        !(px = NAG_ALLOC(npx, double)) || !(py = 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;
  }
  delta = (xhi - xlo) / (npx - 1);
  for (i = 0; i < npx; i++)
    px[i] = MIN(xlo + i * delta, xhi);
  for (i = 0; i < npy; i++)
    py[i] = MIN(ylo + i * delta, yhi);

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

  printf("\nValues of computed spline:\n");
  printf("\n          x");
  for (i = 0; i < npx; i++)
    printf("%7.2f  ", px[i]);
  printf("\n     y\n");
  for (i = npy - 1; i >= 0; i--) {
    printf("%7.2f   ", py[i]);
    for (j = 0; j < npx; ++j)
      printf("%8.2f ", fg[npy * j + i]);
    printf("\n");
  }
END:
  NAG_FREE(spline.lamda);
  NAG_FREE(spline.mu);
  NAG_FREE(spline.c);
  NAG_FREE(warmstartinf.nag_w);
  NAG_FREE(warmstartinf.nag_iw);
  NAG_FREE(f);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(fg);
  NAG_FREE(px);
  NAG_FREE(py);
  return exit_status;
}