/* nag_2d_spline_fit_grid(e02dcc) Example Program
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 *
 * Mark 6 revised, 2000.
 */

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

#define MXMAX 11
#define MYMAX 9
#define NMAX 7
#define NXEST MXMAX + 4
#define NYEST MYMAX + 4

int main(void)
{
  Integer mx, my, nx, ny, npx, npy, i, j;
  double f[MXMAX*MYMAX], x[MXMAX], y[MYMAX], fg[NMAX*NMAX], px[NMAX], py[NMAX];
  double xhi, yhi, xlo, ylo, delta, s, fp;
  Nag_Start start;
  Nag_2dSpline spline;
  Nag_Comm warmstartinf;

  Vprintf("e02dcc Example Program Results\n");
  Vscanf("%*[^\n]"); /* Skip heading in data file */
  /* Input the number of x, y co-ordinates mx, my. */
  Vscanf("%ld%ld",&mx,&my);

  if (mx>0 && mx<=MXMAX && my>0 && my<=MYMAX)
    {
      /* Input the x co-ordinates followed by the y co-ordinates. */
      for (i=0; i<mx;  i++)
        Vscanf("%lf",&x[i]);
      for (i=0; i<my; i++)
        Vscanf("%lf",&y[i]);
      /* Input the mx*my function values f at the grid points. */
      for (i=0; i<mx*my; i++)
        Vscanf("%lf",&f[i]);
      start = Nag_Cold;
      Vscanf("%lf",&s);
      /* Determine the spline approximation. */

      e02dcc(start, mx, x, my, y, f, s, (Integer)(NXEST), (Integer)(NYEST),
             &fp, &warmstartinf, &spline, NAGERR_DEFAULT);
      nx = spline.nx;
      ny = spline.ny;

      Vprintf("\nCalling with smoothing factor s = %11.4e:"
              " spline.nx = %2ld, spline.ny = %2ld.\n\n",s,nx,ny);
      /* Print the knot sets, lamda and mu. */
      Vprintf("Distinct knots in x direction located  at\n");
      for (j=3; j<spline.nx-3; j++)
        Vprintf("%12.4f%s",spline.lamda[j],((j-3)%5==4 || j==spline.nx-4)
                ? "\n" : " ");
      Vprintf("\nDistinct knots in y direction located  at\n");
      for (j=3; j<spline.ny-3; j++)
        Vprintf("%12.4f%s",spline.mu[j],((j-3)%5==4 || j==spline.ny-4)
                ? "\n" : " ");
      Vprintf("\nThe B-spline coefficients:\n\n");
      for (i=0; i<ny-4; i++)
        {
          for (j=0; j<nx-4; j++)
            Vprintf("%9.4f",spline.c[i+j*(ny-4)]);
          Vprintf("\n");
        }
      Vprintf("\nSum of squared residuals fp = %11.4e\n",fp);
      if (fp==0.0)
        Vprintf("\nThe spline is an interpolating spline\n");
      else if (nx==8 && ny==8)
        Vprintf("\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).
       */
      Vscanf("%ld%lf%lf",&npx,&xlo,&xhi);
      Vscanf("%ld%lf%lf",&npy,&ylo,&yhi);
      if (npx<=NMAX && npy<=NMAX)
        {
          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);

          e02dfc(npx, npy, px, py, fg, &spline, NAGERR_DEFAULT);

          Vprintf("\nValues of computed spline:\n");
          Vprintf("\n          x");
          for (i=0; i<npx; i++)
            Vprintf("%7.2f  ",px[i]);
          Vprintf("\n     y\n");
          for (i=npy-1; i>=0; i--)
            {
              Vprintf("%7.2f   ",py[i]);
              for (j=0; j<npx; ++j)
                Vprintf("%8.2f ",fg[npy*j+i]);
              Vprintf("\n");
            }
          NAG_FREE(spline.lamda);
          NAG_FREE(spline.mu);
          NAG_FREE(spline.c);
          NAG_FREE(warmstartinf.nag_w);
          NAG_FREE(warmstartinf.nag_iw);
        }
      else
        {
          Vfprintf(stderr,"npx or npy is out of range: npx = %4ld,"
                   "npy = %4ld\n", npx,npy);
          NAG_FREE(spline.lamda);
          NAG_FREE(spline.mu);
          NAG_FREE(spline.c);
          NAG_FREE(warmstartinf.nag_w);
          NAG_FREE(warmstartinf.nag_iw);
          return EXIT_FAILURE;
        }
      return EXIT_SUCCESS;
    }
  else
    {
      Vfprintf(stderr,"mx or my is out of range: mx = %4ld, my = %4ld\n",mx,my);
      return EXIT_FAILURE;
    }
}