/* nag_quad_md_sgq_multi_vec (d01esc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

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

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL f(Integer ni, Integer ndim, Integer nx, double xtr,
                         Integer nntr, const Integer *icolzp,
                         const Integer *irowix, const double *xs,
                         const Integer *qs, double *fm,
                         Integer *iflag, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

/* We define some structures to serve as a demonstration of safely operating
 * with the NAG communication structure comm when running in parallel.
 */

/* par_sh: any information to be shared between threads in the function f. */
typedef struct
{
  double s_tr;
} par_sh;

/* par_pr: any private workspace that a single thread will require in the
 * execution of the function f.
 */
typedef struct
{
  double *s;
  double *logs;
} par_pr;

/* parallel_comm: a container for par_sh and par_pr. */
typedef struct
{
  par_sh shared;
  par_pr *tprivate;
} parallel_comm;

int main(void)
{
  Integer exit_status = 0;
  Integer ndim, ni;
  Integer maxnx, smpthd, lcvalue;
  double rvalue;
  char cvalue[16];
  Integer *ivalid, *iopts, *maxdlv;
  double *dinest, *errest, *opts;
  parallel_comm parcom;
  int i, thdnum;
  /* Nag Types */
  Nag_VariableType optype;
  Nag_Comm comm;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_quad_md_sgq_multi_vec (d01esc) Example Program Results\n");

  ni = 10;
  ndim = 4;
  lcvalue = 16;
  if (!(iopts = NAG_ALLOC(100, Integer)) ||
      !(opts = NAG_ALLOC(100, double)) ||
      !(maxdlv = NAG_ALLOC(ndim, Integer)) ||
      !(dinest = NAG_ALLOC(ni, double)) ||
      !(errest = NAG_ALLOC(ni, double)) || !(ivalid = NAG_ALLOC(ni, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize option arrays. */
  nag_quad_opt_set("Initialize = d01esc", iopts, 100, opts, 100, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Set any required options. */
  nag_quad_opt_set("Absolute Tolerance = 0.0", iopts, 100, opts, 100, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_quad_opt_set("Relative Tolerance = 1.0e-3", iopts, 100, opts, 100,
                   &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_quad_opt_set("Maximum Level = 6", iopts, 100, opts, 100, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  nag_quad_opt_set("Index Level = 5", iopts, 100, opts, 100, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Set any required maximum dimension levels. */
  for (i = 0; i < ndim; ++i)
    maxdlv[i] = 0;

  /* As a demonstration of safely operating with the communication structure
   * comm when running in parallel, we will create an instance of our
   * parallel_comm structure with fields indexed by the current thread number.
   * The size of the array subfields in this structure are a function of
   * Maximum Nx.
   */
  nag_quad_opt_get("Maximum Nx", &maxnx, &rvalue, cvalue, lcvalue, &optype,
                   iopts, opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_opt_get (d01zlc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  smpthd = nag_omp_get_max_threads();

  /* Allocate an array of smpthd pointers to private structures. */
  if (!(parcom.tprivate = NAG_ALLOC(smpthd, par_pr))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* For each thread, allocate a double array of size maxnx in the 
   * private component of the parallel structure.
   */
  for (thdnum = 0; thdnum < smpthd; thdnum++) {
    if (!(parcom.tprivate[thdnum].s = NAG_ALLOC(maxnx, double)) ||
        !(parcom.tprivate[thdnum].logs = NAG_ALLOC(maxnx, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }

  /* Finally, store the parallel structure in comm for communication to f. */
  comm.p = &parcom;

  /* Approximate the integrals. */
  nag_quad_md_sgq_multi_vec(ni, ndim, f, maxdlv, dinest, errest, ivalid,
                            iopts, opts, &comm, &fail);
  if (fail.code != NE_NOERROR && fail.code != NE_ACCURACY
      && fail.code != NE_TOTAL_PRECISION_LOSS && fail.code != NE_USER_STOP) {
    /* If internal memory allocation failed consider reducing the options
     * 'Maximum Nx' and 'Index Level', or run with fewer threads.
     */
    printf("Error from nag_quad_md_sgq_multi_vec (d01esc).\n%s\n",
           fail.message);
    exit_status = 2;
  }
  else {
    /* NE_NOERROR:
     *   The result returned satisfies the requested accuracy requirements.
     * NE_ACCURACY, NE_TOTAL_PRECISION_LOSS:
     *   The result returned is inaccurate for at least one integral.
     * NE_USER_STOP:
     *   Exit was requested by setting iflag negative in f.
     *   A result will be returned if at least one call to f was successful.
     */
    printf("Integral # | Estimated value | Error estimate |"
           " Final state of integral\n");
    for (i = 0; i < ni; ++i)
      printf("%11d|%17.5e|%16.5e|%8" NAG_IFMT "\n",
             i, dinest[i], errest[i], ivalid[i]);
  }
  for (thdnum = 0; thdnum < smpthd; ++thdnum) {
    NAG_FREE(parcom.tprivate[thdnum].s);
    NAG_FREE(parcom.tprivate[thdnum].logs);
  }

END:
  NAG_FREE(maxdlv);
  NAG_FREE(dinest);
  NAG_FREE(errest);
  NAG_FREE(ivalid);
  NAG_FREE(iopts);
  NAG_FREE(opts);
  NAG_FREE(parcom.tprivate);
  return exit_status;
}

static void NAG_CALL f(Integer ni, Integer ndim, Integer nx, double xtr,
                       Integer nntr, const Integer *icolzp,
                       const Integer *irowix, const double *xs,
                       const Integer *qs, double *fm,
                       Integer *iflag, Nag_Comm *comm)
{
  Integer i, j, k, tid;
  double s_tr, s_ntr;
  double *s, *logs;
  parallel_comm *parcom;

  /* For each evaluation point x_i, i = 1, ..., nx, return in fm the computed
   * values of the ni integrals f_j, j = 1, ..., ni defined by
   *
   *   fm(j,i) = f_j(x_i) 
   *                                                   ndim
   *           = sin(j + S(i))*log(S(i)), where S(i) =  Sum  k*x_i(k).
   *                                                    k=1
   *
   * Split the S expression into two components, one involving only the
   * 'trivial' value xtr:
   *
   *          ndim             ndim
   *   S(i) =  Sum  (k*xtr)  +  Sum  (k*(x_i(k)-xtr))
   *           k=1              k=1
   *
   *                ndim*(ndim+1)   ndim
   *        = xtr * ------------- +  Sum  (k*(x_i(k)-xtr))
   *                     2           k=1
   *
   *       := s_tr                + s_ntr(i)
   *
   * By definition the summands in the s_ntr(i) term on the right-hand side
   * are zero for those k outside the range of indices defined in irowix.
   */

  /* As a demonstration of safely operating with the communication structure
   * comm when running in parallel, the p field of comm is itself (a pointer
   * to) an instance of our parallel_comm structure 'partitioned' by the current
   * thread number. Store some of the s_tr and s_ntr computations in these.
   */

  /* The thread number. */
  tid = nag_omp_get_thread_num();

  /* The S and log(S) terms from above, extracted from comm. */
  parcom = (parallel_comm *) comm->p;
  s = (*parcom).tprivate[tid].s;
  logs = (*parcom).tprivate[tid].logs;

  if (*iflag == 0) {
    /* First call: nx=1, no non-trivial dimensions.
     * The constant s_tr can be reused by all subsequent calculations.
     */
    s_tr = 0.5 * xtr * ((double) (ndim * (ndim + 1)));
    parcom->shared.s_tr = s_tr;
    s[0] = s_tr;
    logs[0] = log(s_tr);
  }
  else {
    /* Calculate S(i) = s_tr + s_ntr(i). */
    s_tr = parcom->shared.s_tr;
    for (i = 0; i < nx; ++i) {
      s_ntr = 0.0;
      for (k = icolzp[i]; k < icolzp[i + 1]; ++k)
        s_ntr = s_ntr + ((double) irowix[k - 1]) * (xs[k - 1] - xtr);
      s[i] = s_tr + s_ntr;
      logs[i] = log(s[i]);
    }
  }
  /* Finally we obtain fm(j,:) = sin(j+S(:))*log(S(:)). */
  for (j = 0; j < ni; ++j)
    for (i = 0; i < nx; ++i)
      fm[i * ni + j] = sin(((double) j + 1) + s[i]) * logs[i];
}