/* nag_4d_shep_eval (e01tlc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 23, 2010.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage01.h>
#include <nagg05.h>
#include <math.h>

#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL funct(double x[]);
#ifdef __cplusplus
}
#endif

#define X(I, J) x[I *4 + J]
#define XE(I, J) xe[I *4 + J]

int main(void)
{
  /* Scalars */
  Integer     exit_status, i, m, n, nq, nw, liq, lrq, lstate, subid;
  Integer     lseed = 1;
  double      fun;
  Nag_BaseRNG genid;
  NagError    fail;
  /* Arrays */
  double      *f = 0, *q = 0, *qx = 0, *rq = 0, *xe = 0, *x = 0;
  Integer     *iq = 0, *state = 0;
  Integer     seed[1], seed2[1];

  exit_status = 0;

  INIT_FAIL(fail);

  printf("nag_4d_shep_eval (e01tlc) Example Program Results\n");

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

  /* Input the seeds. */
  scanf("%ld%ld%*[^\n] ", &seed[0], &seed2[0]);

  /* Choose the base generator */
  genid = Nag_Basic;
  subid = 0;

  /* Get the length of the state array */
  lstate = -1;
  nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
            fail.message);
    exit_status = 1;
    goto END;
  }

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

  /* Allocate memory */
  lrq = 21 * m + 11;
  liq = 2 * m + 1;
  if (!(f = NAG_ALLOC(m, double)) ||
      !(x = NAG_ALLOC(m*4, double)) ||
      !(rq = NAG_ALLOC(lrq, double)) ||
      !(iq = NAG_ALLOC(liq, Integer)) ||
      !(state = NAG_ALLOC(lstate, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Initialise the generator to a repeatable sequence */
  nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",fail.message);
    exit_status = 1;
    goto END;
  }

  /* Generate the data points X */
  nag_rand_basic(m*4, state, x, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_basic (g05sac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Evaluate F */
  for (i = 0; i < m; ++i) {
    f[i] = funct(&X(i,0));
  }

  /* Generate the interpolant. */
  nq = 0;
  nw = 0;

  /* nag_4d_shep_interp (e01tkc).
   * Interpolating functions, modified Shepard's method, four
   * variables
   */
  nag_4d_shep_interp(m, x, f, nw, nq, iq, rq, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_4d_shep_interp (e01tkc).\n%s\n",fail.message);
    exit_status = 1;
    goto END;
  }

  /* Input the number of evaluation points. */
  scanf("%ld%*[^\n] ", &n);

  /* Allocate memory for nag_4d_shep_eval (e01tlc) */
  if (!(q = NAG_ALLOC(n, double)) ||
      !(qx = NAG_ALLOC(n*4, double)) ||
      !(xe = NAG_ALLOC(n*4, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Generate repeatable evaluation points. */
  nag_rand_init_repeatable(genid, subid, seed2, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",fail.message);
    exit_status = 1;
    goto END;
  }
  nag_rand_basic(n*4, state, xe, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_basic (g05sac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* nag_4d_shep_eval (e01tlc).
   * Interpolated values, evaluate interpolant and first derivatives
   * computed by nag_4d_shep_interp (e01tkc).
   */
  fail.print = Nag_TRUE;
  nag_4d_shep_eval(m, x, f, iq, rq, n, xe, q, qx, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_4d_shep_eval (e01tlc).\n%s\n",fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n     i     f(x)      q(x)   |f(x)-q(x)|\n");
  for (i = 0; i < n; ++i) {
    fun = funct(&XE(i,0));
    printf("%6ld%10.4f%10.4f%10.4f\n", i, fun, q[i], fabs(fun-q[i]));
  }

 END:
  NAG_FREE(f);
  NAG_FREE(q);
  NAG_FREE(qx);
  NAG_FREE(rq);
  NAG_FREE(xe);
  NAG_FREE(x);
  NAG_FREE(iq);
  NAG_FREE(state);

  return exit_status;
}

static double NAG_CALL funct(double x[])
{
  /* Scalars */
  double ret_val;

  ret_val = ((1.25+cos(5.4*x[3]))*cos(6.0*x[0])*cos(6.0*x[1]))/
    (6.0*(1.0+pow((3.0*x[2]-1.0),2.0)));
  return ret_val;
}