/* nag_interp_dimn_grid (e01zac) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */
#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;
}