/* nag_quad_md_simplex (d01pac) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>

#ifdef __cplusplus
extern "C" {
#endif
  static double NAG_CALL functn(Integer ndim, const double x[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
#define VERT(I, J) vert[(J-1)* (ndim+1) + I-1]
  /* Scalars */
  Integer  exit_status = 0;
  Integer  i, j, maxord, minord, mxord, ndim;
  double   esterr;
  /* Arrays */
  Integer  iuser[1];
  double   *finvls = 0, *vert = 0;
  /* Nag types */
  NagError fail;
  Nag_Comm comm;

  INIT_FAIL(fail);

  printf("nag_quad_md_simplex (d01pac) Example Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Input mxord and ndim */
  scanf("%ld %ld%*[^\n] ", &mxord, &ndim);

  if (!(finvls = NAG_ALLOC(mxord, double)) ||
      !(vert = NAG_ALLOC(2*(ndim+1)*(ndim+1), double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  for (i = 1; i <= ndim+1; i++)
    for (j = 1; j <= ndim; j++) VERT(i, j) = 0.0;
  for (j = 2; j <= ndim+1; j++) VERT(j, j - 1) = 1.0;

  minord = 0;
  iuser[0] = 0; /* Function counter */
  comm.iuser = iuser;
  for (maxord = 1; maxord <= mxord; maxord++)
    {
      /* nag_quad_md_simplex (d01pac).
       * Multidimensional quadrature over an n-simplex.
       */
      nag_quad_md_simplex(ndim, vert, functn, &minord, maxord, finvls,
                          &esterr, &comm, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_quad_md_simplex (d01pac).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

      if (maxord == 1)
        printf("maxord   Estimated      Estimated         Integrand\n"
               "           value         accuracy        evaluations\n");
      printf("%4ld%13.5f%16.3e%15ld\n",
             maxord, finvls[maxord-1], esterr, comm.iuser[0]);
    }

 END:
  if (finvls) NAG_FREE(finvls);
  if (vert) NAG_FREE(vert);

  return exit_status;
}

static double  NAG_CALL functn(Integer ndim, const double x[], Nag_Comm *comm)
{
  comm->iuser[0]++;
  return exp(x[0] + x[1] + x[2]) * cos(x[0] + x[1] + x[2]);
}