/* nag_2d_spline_fit_grid (e02dcc) Example Program.
 *
 * Copyright 2014 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, 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);

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

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

  printf("nag_2d_spline_fit_grid (e02dcc) Example Program Results\n");
  scanf("%*[^\n]"); /* Skip heading in data file */
  /* Input the number of x, y co-ordinates mx, my. */
  scanf("%ld%ld", &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_2d_spline_fit_grid (e02dcc).
   * Least-squares bicubic spline fit with automatic knot
   * placement, two variables (rectangular grid)
   */
  nag_2d_spline_fit_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_2d_spline_fit_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 = %2ld, spline.ny = %2ld.\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("%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(xlo+i*delta, xhi);
  for (i = 0; i < npy; i++)
    py[i] = MIN(ylo+i*delta, yhi);

  /* 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");
  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;
}