/* nag_2d_spline_fit_panel  (e02dac) Example Program.
 *
 * Copyright 2005 Numerical Algorithms Group.
 *
 * Mark 8, 2004.
 */

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

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

  printf("nag_2d_spline_fit_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("%ld%*[^\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("%ld%ld%*[^\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_2d_panel_sort (e02zac).
       * Sort two-dimensional data into panels for fitting bicubic
       * splines
       */
      nag_2d_panel_sort(px, py, lamda, mu, m, x, y, point, &fail);
      if (fail.code != NE_NOERROR)
	{
	  printf("Error from nag_2d_panel_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_2d_spline_fit_panel (e02dac).
       * Least-squares surface fit, bicubic splines
       */
      nag_2d_spline_fit_panel(m, x, y, f, w, point, dl, eps, &sigma, &rank,
			      &spline, &fail);
      if (fail.code != NE_NOERROR)
	{
	  printf("Error from nag_2d_spline_fit_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%5ld\n", rank);

      /* nag_2d_spline_eval (e02dec).
       * Evaluation of bicubic spline, at a set of points
       */
      nag_2d_spline_eval(m, x, y, ff, &spline, &fail);
      if (fail.code != NE_NOERROR)
	{
	  printf("Error from nag_2d_spline_eval (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:
      if (dl) NAG_FREE(dl);
      if (f) NAG_FREE(f);
      if (ff) NAG_FREE(ff);
      if (lamda) NAG_FREE(lamda);
      if (mu) NAG_FREE(mu);
      if (w) NAG_FREE(w);
      if (x) NAG_FREE(x);
      if (y) NAG_FREE(y);
      if (point) NAG_FREE(point);
      if (spline.lamda) NAG_FREE(spline.lamda);
      if (spline.mu) NAG_FREE(spline.mu);
      if (spline.c) NAG_FREE(spline.c);
    }
  return exit_status;
}