/* nag_2d_spline_fit_scat (e02ddc) Example Program.
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 *
 * Mark 6 revised, 2000.
 * Mark 8 revised, 2004.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage02.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);

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

  nxest = 14;
  nyest = 14;
  printf("nag_2d_spline_fit_scat (e02ddc) Example Program Results\n");
  scanf("%*[^\n]"); /* Skip heading in data file */
  /* Input the number of data-points m. */
  scanf("%ld", &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_2d_spline_fit_scat (e02ddc).
       * Least-squares bicubic spline fit with automatic knot
       * placement, two variables (scattered data)
       */
      nag_2d_spline_fit_scat(start, m, x, y, f, weights, s, nxest, nyest, &fp,
                             &rank, &warmstartinf, &spline, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_2d_spline_fit_scat (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 = %1ld,"
              " ny = %1ld\n", s, nx, ny);

      printf("rank deficiency = %1ld\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("%ld%lf%lf", &npx, &xlo, &xhi);
      scanf("%ld%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_2d_spline_eval_rect (e02dfc).
       * Evaluation of bicubic spline, at a mesh of points
       */
      nag_2d_spline_eval_rect(npx, npy, px, py, 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;
        }
      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);
      if (fg) NAG_FREE(fg);
      if (px) NAG_FREE(px);
      if (py) NAG_FREE(py);
    }
 END:
  if (f) NAG_FREE(f);
  if (weights) NAG_FREE(weights);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  return exit_status;
}