/* nag_2d_spline_deriv_rect (e02dhc) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <math.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage02.h>
#include <nagx04.h>

#ifdef __cplusplus
extern "C" {
#endif

static void NAG_CALL print_spline(Integer *ngx, double *gridx, Integer *ngy,
				  double *gridy, double *z, double *zder,
				  Integer *exit_status);

#ifdef __cplusplus
}
#endif

#define F(I,J) f[my*(I)+(J)]

int main(void)
{
  /* Scalars */
  Integer      exit_status = 0;
  Integer      i, j, mx, my, ngx, ngy, nux, nuy, nxest, nyest;
  double       delta, fp, s, xhi, xlo, yhi, ylo;
  /* Arrays */
  double       *f = 0, *gridx = 0, *gridy = 0, *x = 0, *y = 0, *z = 0,
               *zder = 0;
  /* NAG types */
  Nag_2dSpline spline;
  Nag_Comm     warmstartinf;
  Nag_Start    startc;
  NagError     fail;

  INIT_FAIL(fail);

  printf("nag_2d_spline_deriv_rect (e02dhc) Example Program Results\n");

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

  /* Input the number of X, Y co-ordinates MX, MY.*/
  scanf("%ld %ld%*[^\n]", &mx, &my);
  nxest = mx + 4;
  nyest = my + 4;
  spline.nx = 4;
  spline.ny = 4;

  /* Alocations for spline fitting */
  if (!(x = NAG_ALLOC((mx), double)) ||
      !(y = NAG_ALLOC((my), double)) ||
      !(f = NAG_ALLOC((mx * my), double))
      )
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  for (i = 0; i < mx; i++) scanf("%lf", &x[i]);
  scanf("%*[^\n]");
  for (i = 0; i < my; i++) scanf("%lf", &y[i]);
  scanf("%*[^\n]");

  /* Input the MX*MY function values F at grid points and smoothing factor.*/
  for (i = 0; i < mx; i++)
    for (j = 0; j < my; j++)
      scanf("%lf", &F(i,j));
  
  scanf("%*[^\n]");
  scanf("%lf%*[^\n]", &s);

  /* nag_2d_spline_fit_grid (e02dcc).
   * Least squares bicubic spline fit with automatic knot placement,
   * two variables (rectangular grid)
   */
  startc = Nag_Cold;
  nag_2d_spline_fit_grid(startc, mx, x, my, y, f, s, nxest, nyest, &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;
    }

  printf("\nSpline fit used smoothing factor s = %13.4e.\n", s);
  printf("Number of knots in each direction  = %5ld,%5ld.\n\n",
	 spline.nx, spline.ny);
  printf("Sum of squared residuals           = %13.4e.\n", fp);

  /* Spline and its derivative to be evaluated on rectangular grid with
   * ngx*ngy points on the domain [xlo,xhi] by [ylo,yhi].
   */
  scanf("%ld%lf%lf%*[^\n]", &ngx, &xlo, &xhi);
  scanf("%ld%lf%lf%*[^\n]", &ngy, &ylo, &yhi);

  /* Allocations for spline evaluation.*/
  if (!(gridx = NAG_ALLOC((ngx), double)) ||
      !(gridy = NAG_ALLOC((ngy), double)) ||
      !(z = NAG_ALLOC((ngx * ngy), double)) ||
      !(zder = NAG_ALLOC((ngx * ngy), double))
      )
    {
      printf("Allocation failure\n");
      exit_status = -2;
      goto END;
    }

  delta = (xhi - xlo)/(double) (ngx - 1);
  gridx[0] = xlo;
  for (i = 1; i < ngx - 1; i++) gridx[i] = gridx[i-1] + delta;
  gridx[ngx-1] = xhi;

  delta = (yhi - ylo)/(double) (ngy - 1);
  gridy[0] = ylo;
  for (i = 1; i < ngy - 1; i++) gridy[i] = gridy[i-1] + delta;
  gridy[ngy-1] = yhi;

  /* Evaluate spline (nux=nuy=0) using
   * nag_2d_spline_deriv_rect (e02dhc).
   * Evaluation of spline surface at mesh of points with derivatives
   */
  nux = 0;
  nuy = 0;
  nag_2d_spline_deriv_rect(ngx, ngy, gridx, gridy, nux, nuy, z, &spline, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_2d_spline_deriv_rect (e02dhc))\n%s\n",
	     fail.message);
      exit_status = 2;
      goto END;
    }


  /* Evaluate spline partial derivative of order (nux,nuy)*/
  scanf("%ld%ld%*[^\n]", &nux, &nuy);
  printf("\nDerivative of spline has order nux, nuy =%5ld, %5ld.\n",
	 nux, nuy);
  /* nag_2d_spline_deriv_rect (e02dhc).
   * Evaluation of spline surface at mesh of points with derivatives
   */
  nag_2d_spline_deriv_rect(ngx, ngy, gridx, gridy, nux, nuy, zder, &spline,
			   &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_2d_spline_deriv_rect (e02dhc)\n%s\n",
	     fail.message);
      exit_status = 3;
      goto END;
    }

  /* Print tabulated spline and derivative evaluations.*/
  print_spline(&ngx, gridx, &ngy, gridy, z, zder, &exit_status);

 END:
  if (f) NAG_FREE(f);
  if (gridx) NAG_FREE(gridx);
  if (gridy) NAG_FREE(gridy);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (z) NAG_FREE(z);
  if (zder) NAG_FREE(zder);
  if (spline.lamda) NAG_FREE(spline.lamda);
  if (spline.mu) NAG_FREE(spline.mu);
  if (spline.c) NAG_FREE(spline.c);
  if (warmstartinf.nag_w) NAG_FREE(warmstartinf.nag_w);
  if (warmstartinf.nag_iw) NAG_FREE(warmstartinf.nag_iw);
  return exit_status;
}

static void print_spline(Integer *ngx, double *gridx, Integer *ngy,
			 double *gridy, double *z, double *zder,
			 Integer *exit_status)
{
  /* Print spline function and spline derivative evaluation*/
  Integer        indent = 0, ncols = 80;
  char           formc[] = "%8.3f";
  Integer        i;
  char           title[49];
  char           *outfile = 0, *clabs = 0, *rlabs = 0;
  char           **clabsc = 0, **rlabsc = 0;
  Nag_OrderType  order;
  Nag_MatrixType matrixc;
  Nag_DiagType   diagc;
  Nag_LabelType  chlabelc;
  NagError       fail;

  INIT_FAIL(fail);

  /* Allocate for row and column label*/
  if (
      !(clabs = NAG_ALLOC((*ngx) * 10+1, char)) ||
      !(rlabs = NAG_ALLOC((*ngy) * 10+1, char)) ||
      !(clabsc = NAG_ALLOC(*ngx, char *)) ||
      !(rlabsc = NAG_ALLOC(*ngy, char *))
      )
    {
      printf("Allocation failure\n");
      *exit_status = -3;
      goto END;
    }
  for (i = 0; i < *ngx; i++)
    {
      clabsc[i] = NAG_ALLOC(11, char);
      sprintf(clabsc[i], "%5.2f%5s", gridx[i], "");
    }
  for (i = 0; i < *ngy; i++)
    {
      rlabsc[i] = NAG_ALLOC(11, char);
      sprintf(rlabsc[i], "%5.2f%5s", gridy[i], "");
    }
  /* Generate column and row labels to print the results with.*/
  memset(clabs, ' ', *ngx * 10 + 1);
  memset(rlabs, ' ', *ngy * 10 + 1);

  for (i = 0; i < *ngx; i++) sprintf(&clabs[i*10], "%5.2f%5s", gridx[i], "");
  for (i = 0; i < *ngy; i++) sprintf(&rlabs[i*10], "%5.2f%5s", gridy[i], "");

  order = Nag_ColMajor;
  matrixc = Nag_GeneralMatrix;
  diagc = Nag_NonUnitDiag;
  chlabelc = Nag_CharacterLabels;

  /* Print the spline evaluations, z. */
  strcpy(title, "Spline evaluated on X-Y grid (X across, Y down):");
  printf("\n");
  /* nag_gen_real_mat_print_comp (x04cbc).
   * Print real general matrix (comprehensive)
   */
  nag_gen_real_mat_print_comp(order, matrixc, diagc, *ngy, *ngx, z, *ngy, formc,
                              title, chlabelc, (const char **)rlabsc, chlabelc, 
			      (const char **)clabsc, ncols,
			      indent, outfile, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_gen_real_mat_print_comp (x04cbc)\n%s\n",
	     fail.message);
      *exit_status = 4;
      goto END;
    }

  /* Print the spline derivative evaluations, zder. */
  strcpy(title, "Spline derivative evaluated on X-Y grid:");
  printf("\n");

  nag_gen_real_mat_print_comp(order, matrixc, diagc, *ngy, *ngx, zder, *ngy,
			      formc, title, chlabelc, (const char **)rlabsc,
			      chlabelc, (const char **)clabsc, ncols, indent,
			      outfile, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_gen_real_mat_print_comp (x04cbc)\n%s\n",
	     fail.message);
      *exit_status = 5;
      goto END;
    }

 END:
  for (i = 0; i < *ngy; i++) if (rlabsc[i]) NAG_FREE(rlabsc[i]);
  if (rlabsc) NAG_FREE(rlabsc);
  if (rlabs) NAG_FREE(rlabs);
  for (i = 0; i < *ngx; i++) if (clabsc[i]) NAG_FREE(clabsc[i]);
  if (clabsc) NAG_FREE(clabsc);
  if (clabs) NAG_FREE(clabs);
}