/* nag_fit_dim2_spline_sctr (e02ddc) 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, m, npx, npy, nx, nxest, ny, nyest, rank;
  NagError fail;
  Nag_2dSpline spline;
  Nag_Start start;
  double delta, *f = 0, *fg = 0, fp, *px = 0, *py = 0, s, warmstartinf;
  double *weights = 0, *x = 0, xhi, xlo, *y = 0, yhi, ylo;

  INIT_FAIL(fail);

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

  nxest = 14;
  nyest = 14;
  printf("nag_fit_dim2_spline_sctr (e02ddc) Example Program Results\n");
  scanf("%*[^\n]"); /* Skip heading in data file */
  /* Input the number of data-points m. */
  scanf("%" NAG_IFMT "", &m);
  if (m >= 16) {
    if (!(f = NAG_ALLOC(m, double)) || !(weights = NAG_ALLOC(m, double)) ||
        !(x = NAG_ALLOC(m, double)) || !(y = NAG_ALLOC(m, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  } else {
    printf("Invalid m.\n");
    exit_status = 1;
    return exit_status;
  }
  /* Input the data-points and the weights. */
  for (i = 0; i < m; i++)
    scanf("%lf%lf%lf%lf", &x[i], &y[i], &f[i], &weights[i]);
  start = Nag_Cold;
  if (scanf("%lf", &s) != EOF) {
    /* Determine the spline approximation. */

    /* nag_fit_dim2_spline_sctr (e02ddc).
     * Least squares bicubic spline fit with automatic knot
     * placement, two variables (scattered data)
     */
    nag_fit_dim2_spline_sctr(start, m, x, y, f, weights, s, nxest, nyest, &fp,
                             &rank, &warmstartinf, &spline, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_dim2_spline_sctr (e02ddc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
    nx = spline.nx;
    ny = spline.ny;
    printf("\nCalling with smoothing factor s = %13.4e, nx = %1" NAG_IFMT ","
           " ny = %1" NAG_IFMT "\n",
           s, nx, ny);

    printf("rank deficiency = %1" NAG_IFMT "\n\n", (nx - 4) * (ny - 4) - rank);

    /* 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.2f", spline.c[i + j * (ny - 4)]);
      printf("\n");
    }

    printf("\n Sum of squared residuals fp = %13.4e\n", fp);
    if (nx == 8 && ny == 8)
      printf("The 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(xhi, xlo + i * delta);
    for (i = 0; i < npy; i++)
      py[i] = MIN(yhi, ylo + i * delta);

    /* 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\n");
    printf("          x");
    for (i = 0; i < npx; i++)
      printf("%8.2f ", px[i]);
    printf("\n     y\n");
    for (i = npy - 1; i >= 0; i--) {
      printf("%8.2f   ", py[i]);
      for (j = 0; j < npx; j++)
        printf("%8.2f ", fg[npy * j + i]);
      printf("\n");
    }
    /* Free memory used by spline */
    NAG_FREE(spline.lamda);
    NAG_FREE(spline.mu);
    NAG_FREE(spline.c);
    NAG_FREE(fg);
    NAG_FREE(px);
    NAG_FREE(py);
  }
END:
  NAG_FREE(f);
  NAG_FREE(weights);
  NAG_FREE(x);
  NAG_FREE(y);
  return exit_status;
}