/* nag_quad_md_simplex (d01pac) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.1, 2024.
 */

#include <math.h>
#include <nag.h>
#include <stdio.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("%" NAG_IFMT " %" NAG_IFMT "%*[^\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("%4" NAG_IFMT "%13.5f%16.3e%15" NAG_IFMT "\n", maxord,
           finvls[maxord - 1], esterr, comm.iuser[0]);
  }

END:
  NAG_FREE(finvls);
  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]);
}