/* nag_interp_dimn_grid (e01zac) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.1, 2024.
 */
#include <math.h>
#include <nag.h>

static double prod(Integer *x, Integer a, Integer b);

int main(void) {
  /* Scalars */
  int i, j, m;
  Integer d, lx, exit_status = 0;
  double ans, wf, tv;
  /* Arrays */
  double *v = 0, *point = 0, *axis = 0;
  Integer *narr = 0;
  const char *methodname[3];
  /* Nag Types */
  NagError fail;
  Nag_Interp method;
  Nag_Boolean uniform = Nag_TRUE;

  methodname[0] = "Weighted Average";
  methodname[1] = "Linear Method   ";
  methodname[2] = "Cubic Method    ";

  INIT_FAIL(fail);

  printf("nag_interp_dimn_grid (e01zac) Example Program Results\n\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Input the number of nodes. */
  scanf("%" NAG_IFMT "%*[^\n] ", &d);

  lx = 2 * d;
  if (!(point = NAG_ALLOC(d, double)) || !(axis = NAG_ALLOC(lx, double)) ||
      !(narr = NAG_ALLOC(d, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Input narr and point. */
  for (i = 0; i < d; i++) {
    scanf("%" NAG_IFMT, &narr[i]);
  }
  scanf("%*[^\n] ");
  for (i = 0; i < d; i++) {
    scanf("%lf", &point[i]);
  }
  scanf("%*[^\n] ");
  for (i = 0; i < d; i++) {
    scanf("%lf", &axis[i * 2]);
    scanf("%lf%*[^\n]", &axis[i * 2 + 1]);
  }
  /* Allocate v. */
  if (!(v = NAG_ALLOC(prod(narr, 0, d), double))) {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }

  for (i = 0; i < prod(narr, 1, d); i++) {
    for (j = 0; j < narr[0]; j++) {
      scanf("%lf", &v[narr[0] * i + j]);
    }
    scanf("%*[^\n]");
  }
  scanf("%lf", &wf);

  /* get function value at point */
  tv = point[0] * point[0] * point[0] - point[1] * point[1] + point[2];
  printf("True Value: %10.4f\n\n", tv);

  printf("    method      |   answer   |    error \n");
  printf("-----------------------------------------\n");
  for (m = 1; m <= 3; m++) {
    /* nag_interp_dimn_grid (e01zac).
     * Interpolated values, gridded data, n dimensions.
     */
    if (m == 1) {
      method = Nag_WeightedAverage;
    } else if (m == 2) {
      method = Nag_LinearInterp;
    } else if (m == 3) {
      method = Nag_CubicInterp;
    }

    nag_interp_dimn_grid(d, narr, uniform, axis, lx, v, point, method, 1, wf,
                         &ans, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_interp_dimn_grid (e01zac) using %s.\n%s\n",
             methodname[m - 1], fail.message);
      exit_status = 1;
      goto END;
    } else {
      printf("%s| %10.4f | %10.4f\n", methodname[m - 1], ans, fabs(ans - tv));
    }
  }

END:

  NAG_FREE(v);
  NAG_FREE(axis);
  NAG_FREE(point);
  NAG_FREE(narr);

  return exit_status;
}

static double prod(Integer *x, Integer a, Integer b) {
  int i;
  Integer ans = 1;
  for (i = a; i < b; i++) {
    ans *= x[i];
  }
  return ans;
}