/* 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__1, eps, sigma, sum, temp;
  Integer exit_status, 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, *y=0;
  Integer *point=0;

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

  exit_status = 0;
  INIT_FAIL(fail);

  Vprintf( "nag_2d_spline_fit_panel (e02dac) Example Program Results\n");

  /* Skip heading in data file */
  Vscanf("%*[^\n] ");

  while (scanf("%lf", &eps) != EOF)
    {

      /* Read data, interchanging X and Y axes if PX.LT.PY */
      Vscanf("%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))
               )
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

        }
      else
        {
          Vprintf("Invalid m.\n");
          exit_status = 1;
          return exit_status;
        }
      Vscanf("%ld%ld%*[^\n] ",   &px,  &py);
      if (px >= 8 && py >= 8)
        {
          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)) )
            {
              Vprintf("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))
                   )
                {
                  Vprintf("Allocation failure\n");
                  exit_status = -1;
                  goto END;
                }

              for (i = 1; i <= m; ++i)
                {
                  Vscanf("%lf%lf%lf%lf", &y[i - 1], &x[i - 1], &f[i - 1],
                         &w[i - 1]);
                }
              Vscanf("%*[^\n] ");

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

              if (px > 8)
                {

                  for (j = 5; j <= px - 4; ++j)
                    {
                      Vscanf("%lf", &lamda[j - 1]);
                    }
                  Vscanf("%*[^\n] ");

                }
              if (py > 8)
                {
                  for (j = 5; j <= py - 4; ++j)
                    {
                      Vscanf("%lf", &mu[j - 1]);
                    }
                  Vscanf("%*[^\n] ");
                }
            }
          Vprintf("\n%s%1.1s %s\n",   "Interior ",   label + itemp,  "-knots");
          for (j = 5; j <= px - 4; ++j)
            {
              Vprintf("%11.4f\n",   lamda[j - 1]);
            }
          if (px == 8)
            {
              Vprintf("%s\n",   "None");
            }
          Vprintf("\n%s%1.1s %s\n", "Interior ", label + (2 - itemp - 1),
                  "-knots");
          for (j = 5; j <= py - 4; ++j)
            {
              Vprintf("%1s%11.4f\n",   "",  mu[j - 1]);
            }
          if (py == 8)
            {

              Vprintf("%s\n",   "None");
            }
          /* Sort points into panel order */

          /* 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)
            {
              Vprintf("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)) )
            {
              Vprintf("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)
            {
              Vprintf("Error from nag_2d_spline_fit_panel (e02dac).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          Vprintf("\n%s%16.3e\n",   "Sum of squares of residual RHS",  sigma);
          Vprintf("\n%s%5ld\n",   "Rank",  rank);

          /* Evaluate spline at the data points */
          /* 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)
            {
              Vprintf("Error from nag_2d_spline_eval (e02dec).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          sum = 0.;
          if (itemp == 1)
            {
              Vprintf("\n%s\n",   "X and Y have been interchanged");
            }

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

                  Vprintf("%11.4f%11.4f%11.4f%11.4f%11.3e\n", x[iadres - 1],
                          y[iadres - 1], f[iadres - 1], ff[iadres - 1], temp);
                  /* Computing 2nd power */
                  d__1 = temp * w[iadres - 1];
                  sum += d__1 * d__1;
                }
            }

          Vprintf("\n%s%16.3e\n",   "Sum of squared residuals",  sum);

          Vprintf("\n%s\n",   "Spline coefficients");
          for (i = 1; i <= px - 4; ++i)
            {
              for (j = 1; j <= py - 4; ++j)
                {
                  Vprintf("%11.4f",   spline.c[(i - 1) * (py - 4) + j - 1]);
                }
              Vprintf("\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;
}