/* nag_fit_dim2_spline_panel  (e02dac) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.5, 2022.
 */

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

int main(void) {
  /* Initialized data */
  char label[] = "xy";

  /* Scalars */
  double d, eps, sigma, sum, temp;
  Integer exit_status = 0, i, iadres, itemp, j, m, nc, np, npoint, px, py, rank;

  /* Arrays */
  double *dl = 0, *f = 0, *ff = 0, *lamda = 0, *mu = 0, *w = 0, *x = 0;
  double *y = 0;
  Integer *point = 0;

  /* Nag Types */
  Nag_2dSpline spline;
  NagError fail;

  exit_status = 0;
  INIT_FAIL(fail);

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

  printf("nag_fit_dim2_spline_panel (e02dac) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  while (scanf("%lf", &eps) != EOF && exit_status == 0)
  {
    /* Read data, interchanging X and Y axes if PX.LT.PY */
    scanf("%" NAG_IFMT "%*[^\n] ", &m);
    if (m > 1) {
      /* Allocate memory */
      if (!(f = NAG_ALLOC(m, double)) || !(ff = NAG_ALLOC(m, double)) ||
          !(w = 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;
      goto END;
    }
    scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &px, &py);
    if (px < 8 && py < 8) {
      printf("px or py is too small.\n");
      exit_status = 1;
      goto END;
    }
    nc = (px - 4) * (py - 4);
    np = (px - 7) * (py - 7);
    npoint = m + (px - 7) * (py - 7);

    /* Allocate memory */
    if (!(dl = NAG_ALLOC(nc, double)) ||
        !(point = NAG_ALLOC(npoint, Integer))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    if (px < py) {
      itemp = px;
      px = py;
      py = itemp;
      itemp = 1;
      /* Allocate memory */
      if (!(lamda = NAG_ALLOC(px, double)) || !(mu = NAG_ALLOC(py, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      for (i = 0; i < m; ++i)
        scanf("%lf%lf%lf%lf", &y[i], &x[i], &f[i], &w[i]);
      scanf("%*[^\n] ");

      if (py > 8) {
        for (j = 4; j < py - 4; ++j)
          scanf("%lf", &mu[j]);
        scanf("%*[^\n] ");
      }
      if (px > 8) {
        for (j = 4; j < px - 4; ++j)
          scanf("%lf", &lamda[j]);
        scanf("%*[^\n] ");
      }
    } else {
      /* Allocate memory */
      if (!(lamda = NAG_ALLOC(px, double)) || !(mu = NAG_ALLOC(py, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
      itemp = 0;
      for (i = 0; i < m; ++i)
        scanf("%lf%lf%lf%lf", &x[i], &y[i], &f[i], &w[i]);
      scanf("%*[^\n] ");

      if (px > 8) {
        for (j = 4; j < px - 4; ++j)
          scanf("%lf", &lamda[j]);
        scanf("%*[^\n] ");
      }
      if (py > 8) {
        for (j = 4; j < py - 4; ++j)
          scanf("%lf", &mu[j]);
        scanf("%*[^\n] ");
      }
    }

    printf("\nInterior %1.1s -knots\n", label + itemp);
    for (j = 4; j < px - 4; ++j)
      printf("%11.4f\n", lamda[j]);
    if (px == 8)
      printf("None\n");

    printf("\nInterior %1.1s -knots\n", label + (2 - itemp - 1));
    for (j = 4; j < py - 4; ++j)
      printf("%1s%11.4f\n", "", mu[j]);
    if (py == 8)
      printf("None\n");

    /* nag_fit_dim2_spline_sort (e02zac).
     * Sort two-dimensional data into panels for fitting bicubic
     * splines
     */
    nag_fit_dim2_spline_sort(px, py, lamda, mu, m, x, y, point, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_dim2_spline_sort (e02zac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    /* Fit bicubic spline to data points */
    spline.nx = px;
    spline.ny = py;

    if (!(spline.c = NAG_ALLOC((spline.nx - 4) * (spline.ny - 4), double)) ||
        !(spline.lamda = NAG_ALLOC(spline.nx, double)) ||
        !(spline.mu = NAG_ALLOC(spline.ny, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    for (i = 0; i < spline.nx; i++)
      spline.lamda[i] = lamda[i];
    for (i = 0; i < spline.ny; i++)
      spline.mu[i] = mu[i];

    /* nag_fit_dim2_spline_panel (e02dac).
     * Least squares surface fit, bicubic splines
     */
    nag_fit_dim2_spline_panel(m, x, y, f, w, point, dl, eps, &sigma, &rank,
                              &spline, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_dim2_spline_panel (e02dac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("\nSum of squares of residual RHS%16.3e\n", sigma);
    printf("\nRank%5" NAG_IFMT "\n", rank);

    /* nag_fit_dim2_spline_evalv (e02dec).
     * Evaluation of bicubic spline, at a set of points
     */
    nag_fit_dim2_spline_evalv(m, x, y, ff, &spline, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_fit_dim2_spline_evalv (e02dec).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    sum = 0.;
    if (itemp == 1)
      printf("\nx and y have been interchanged\n\n");

    /*Output data points, fitted values and residuals */
    printf("      X          Y          Data       Fit      Residual\n");
    for (i = 0; i < np; ++i) {
      iadres = i + m;
      while ((iadres = point[iadres] - 1) >= 0) {
        temp = ff[iadres] - f[iadres];

        printf("%11.4f%11.4f%11.4f%11.4f%12.3e\n", x[iadres], y[iadres],
               f[iadres], ff[iadres], temp);
        /* Computing 2nd power */
        d = temp * w[iadres];
        sum += d * d;
      }
    }

    printf("\nSum of squared residuals%16.3e\n", sum);
    printf("\nSpline coefficients\n");
    for (i = 0; i < px - 4; ++i) {
      for (j = 0; j < py - 4; ++j)
        printf("%11.4f", spline.c[i * (py - 4) + j]);
      printf("\n");
    }
  END:
    NAG_FREE(dl);
    NAG_FREE(f);
    NAG_FREE(ff);
    NAG_FREE(lamda);
    NAG_FREE(mu);
    NAG_FREE(w);
    NAG_FREE(x);
    NAG_FREE(y);
    NAG_FREE(point);
    NAG_FREE(spline.lamda);
    NAG_FREE(spline.mu);
    NAG_FREE(spline.c);
  }
  return exit_status;
}