/* nag_fit_dim2_spline_derivm (e02dhc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

#include <math.h>
#include <nag.h>
#include <string.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_fit_dim2_spline_derivm (e02dhc) Example Program Results\n");
  fflush(stdout);

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

  /* Input the number of X, Y co-ordinates MX, MY. */
  scanf("%" NAG_IFMT " %" NAG_IFMT "%*[^\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_fit_dim2_spline_grid (e02dcc).
   * Least squares bicubic spline fit with automatic knot placement,
   * two variables (rectangular grid)
   */
  startc = Nag_Cold;
  nag_fit_dim2_spline_grid(startc, mx, x, my, y, f, s, nxest, nyest, &fp,
                           &warmstartinf, &spline, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_dim2_spline_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  = %5" NAG_IFMT ",%5" NAG_IFMT
         ".\n\n",
         spline.nx, spline.ny);
  printf("Sum of squared residuals           = %13.4e.\n", fp);
  fflush(stdout);

  /* Spline and its derivative to be evaluated on rectangular grid with
   * ngx*ngy points on the domain [xlo,xhi] by [ylo,yhi].
   */
  scanf("%" NAG_IFMT "%lf%lf%*[^\n]", &ngx, &xlo, &xhi);
  scanf("%" NAG_IFMT "%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_fit_dim2_spline_derivm (e02dhc).
   * Evaluation of spline surface at mesh of points with derivatives
   */
  nux = 0;
  nuy = 0;
  nag_fit_dim2_spline_derivm(ngx, ngy, gridx, gridy, nux, nuy, z, &spline,
                             &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_fit_dim2_spline_derivm (e02dhc))\n%s\n",
           fail.message);
    exit_status = 2;
    goto END;
  }

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

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

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

static void NAG_CALL 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;
  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 (!(clabsc = NAG_ALLOC(*ngx, char *)) ||
      !(rlabsc = NAG_ALLOC(*ngy, char *))) {
    printf("Allocation failure\n");
    *exit_status = -3;
    goto END;
  }

  /* Allocate memory to clabsc and rlabsc elements and generate
   * column and row labels to print the results with.
   */
  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], "");
  }

  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");
  fflush(stdout);

  /* nag_file_print_matrix_real_gen_comp (x04cbc).
   * Print real general matrix (comprehensive)
   */
  nag_file_print_matrix_real_gen_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_file_print_matrix_real_gen_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");
  fflush(stdout);

  nag_file_print_matrix_real_gen_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_file_print_matrix_real_gen_comp (x04cbc)\n%s\n",
           fail.message);
    *exit_status = 5;
    goto END;
  }

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