/* nag_quad_1d_gen_vec_multi_rcomm (d01rac) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>
#include <nagx01.h>

/* Print information on splitting and evaluations over subregions? */
Nag_Boolean disp_integration_info = Nag_TRUE;

static void display_integration_details(const Integer ni,
                                        const Integer iopts[],
                                        const double opts[],
                                        const Integer icom[],
                                        const double com[]);
static void display_option(const char *optstr, const Nag_VariableType optype,
                           const Integer ivalue, const double rvalue,
                           const char *cvalue);

int main(void)
{
#define FM(J,I) fm[(I-1)*ldfm + J-1]

  /* Scalars */
  int exit_status = 0;
  Integer len_cvalue;
  double a, b, rvalue;
  Integer irevcm, ivalue, i, j, lcmax, lcmin, lcom, ldfm, ldfmrq,
         lenx, lenxrq, licmax, licmin, licom, liopts, lopts, ni, nx,
         sdfm, sdfmrq, sid;
  /* Arrays */
  char cvalue[17];
  double *com = 0, *dinest = 0, *errest = 0, *fm = 0, *opts = 0, *x = 0;
  Integer *icom = 0, *iopts = 0, *needi = 0;

  /* NAG types */
  Nag_VariableType optype;
  NagError fail;

  printf("nag_quad_1d_gen_vec_multi_rcomm (d01rac) Example Program Results"
         "\n\n");

  /* Setup phase. */
  /* Set problem parameters. */
  ni = 2;
  a = 0.0;
  b = nag_pi;

  liopts = 100;
  lopts = 100;
  if (!(opts = NAG_ALLOC((lopts), double)) ||
      !(iopts = NAG_ALLOC((liopts), Integer))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  INIT_FAIL(fail);
  /* Initialize option arrays using nag_quad_opt_set (d01zkc). */
  nag_quad_opt_set("Initialize = nag_quad_1d_gen_vec_multi_rcomm",
                   iopts, liopts, opts, lopts, &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("Quadrature Rule = gk41", iopts, liopts, opts, lopts,
                   &fail);
  nag_quad_opt_set("Absolute Tolerance = 1.0e-7", iopts, liopts, opts, lopts,
                   &fail);
  nag_quad_opt_set("Relative Tolerance = 1.0e-7", iopts, liopts, opts, lopts,
                   &fail);

  /* Determine required array dimensions for
   * nag_quad_1d_gen_vec_multi_rcomm (d01rac) using
   * nag_quad_1d_gen_vec_multi_dimreq (d01rcc).
   */
  nag_quad_1d_gen_vec_multi_dimreq(ni, &lenxrq, &ldfmrq, &sdfmrq,
                                   &licmin, &licmax, &lcmin, &lcmax,
                                   iopts, opts, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_quad_1d_gen_vec_multi_dimreq (d01rcc).\n%s\n",
           fail.message);
    exit_status = 2;
    goto END;
  }
  ldfm = ldfmrq;
  sdfm = sdfmrq;
  lenx = lenxrq;
  licom = licmax;
  lcom = lcmax;

  /* Allocate remaining arrays. */
  if (!(x = NAG_ALLOC((lenx), double)) ||
      !(needi = NAG_ALLOC((ni), Integer)) ||
      !(fm = NAG_ALLOC((ldfm) * (sdfm), double)) ||
      !(dinest = NAG_ALLOC((ni), double)) ||
      !(errest = NAG_ALLOC((ni), double)) ||
      !(com = NAG_ALLOC((lcom), double)) ||
      !(icom = NAG_ALLOC((licom), Integer))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Solve phase. */
  /* Use nag_quad_1d_gen_vec_multi_rcomm (d01rac) to evaluate the
   * definite integrals of:
   *  f_1 = (x*sin(2*x))*cos(15*x)
   *  f_2 = (x*sin(2*x))*(x*cos(50*x))
   */

  INIT_FAIL(fail);
  /* Set initial irevcm. */
  irevcm = 1;
  while (irevcm) {
    /* nag_quad_1d_gen_vec_multi_rcomm (d01rac).
     * One-dimensional quadrature, adaptive, vectorized, multi-integral,
     * reverse communication.
     */
    nag_quad_1d_gen_vec_multi_rcomm(&irevcm, ni, a, b,
                                    &sid, needi, x, lenx, &nx, fm, ldfm,
                                    dinest, errest,
                                    iopts, opts, icom, licom, com, lcom,
                                    &fail);
    switch (irevcm) {
    case 11:
      /* Initial returns.
       * These will occur during the non-adaptive phase.
       * All values must be supplied.
       * dinest and errest do not contain approximations over the complete
       * interval at this stage.
       */

      for (i = 1; i <= nx; i++) {
        FM(2, i) = x[i - 1] * sin(2.0 * x[i - 1]);
        FM(1, i) = FM(2, i) * cos(15.0 * x[i - 1]);
        /* Complete f_2 calculation. */
        FM(2, i) = FM(2, i) * x[i - 1] * cos(50.0 * x[i - 1]);
      }

      break;
    case 12:
      /* Intermediate returns.
       * These will occur during the adaptive phase.
       * All requested values must be supplied.
       * dinest and errest contain approximations over the complete
       * interval at this stage.
       */
      if ((needi[0] == 1) && (needi[1] == 1)) {
        for (i = 1; i <= nx; i++) {
          FM(2, i) = x[i - 1] * sin(2.0 * x[i - 1]);
          FM(1, i) = FM(2, i) * cos(15.0 * x[i - 1]);
          /* Complete f_2 calculation. */
          FM(2, i) = FM(2, i) * x[i - 1] * cos(50.0 * x[i - 1]);
        }
      }
      else if (needi[0] == 1) {
        /* Only calculation of f_1 is requried. */
        for (i = 1; i <= nx; i++)
          FM(1, i) =
                 (x[i - 1] * sin(2.0 * x[i - 1])) * (cos(15.0 * x[i - 1]));
      }
      else if (needi[1] == 1) {
        /* Only calculation of f_2 is requried. */
        for (i = 1; i <= nx; i++)
          FM(2, i) =
                 (x[i - 1] * sin(2.0 * x[i - 1])) * (x[i - 1] *
                                                     cos(50.0 * x[i - 1]));
      }
      break;
    case 0:
      /* Final return. Test fail.code. */
      switch (fail.code) {
      case NE_NOERROR:
        break;
      case NE_ACCURACY:
        printf("Warning: nag_quad_1d_gen_vec_multi_rcomm (d01rac) has "
               "returned at \n least one error estimate exceeding the"
               " requested tolerances\n");
        break;
      case NE_QUAD_BAD_SUBDIV_INT:
        /* Useful information has been returned. */
        printf("Warning: nag_quad_1d_gen_vec_multi_rcomm (d01rac) has "
               "detected \n extremely bad behaviour for at least"
               " one integral\n");
        break;
      default:;
        /* An unrecoverable error has been detected. */
        printf("Error from nag_quad_1d_gen_vec_multi_rcomm (d01rac).\n%s\n",
               fail.message);
        exit_status = 3;
        goto END;
      }
    }
  }

  /* Query some currently set options using nag_quad_opt_get (d01zlc). */
  len_cvalue = 17;
  nag_quad_opt_get("Quadrature rule", &ivalue, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  display_option("Quadrature rule", optype, ivalue, rvalue, cvalue);
  nag_quad_opt_get("Maximum Subdivisions", &ivalue, &rvalue, cvalue,
                   len_cvalue, &optype, iopts, opts, &fail);
  display_option("Maximum Subdivisions", optype, ivalue, rvalue, cvalue);
  nag_quad_opt_get("Extrapolation", &ivalue, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  display_option("Extrapolation", optype, ivalue, rvalue, cvalue);
  nag_quad_opt_get("Extrapolation Safeguard", &ivalue, &rvalue,
                   cvalue, len_cvalue, &optype, iopts, opts, &fail);
  display_option("Extrapolation safeguard", optype, ivalue, rvalue, cvalue);

  /* Print solution. */
  printf("\n Integral |  needi  |   dinest   |   errest   \n");
  for (j = 1; j <= ni; j++)
    printf(" %9" NAG_IFMT " %9" NAG_IFMT " %12.4e %12.4e\n",
           j, needi[j - 1], dinest[j - 1], errest[j - 1]);

  /* Investigate integration strategy. */
  if (disp_integration_info)
    display_integration_details(ni, iopts, opts, icom, com);

END:
  NAG_FREE(com);
  NAG_FREE(dinest);
  NAG_FREE(errest);
  NAG_FREE(fm);
  NAG_FREE(opts);
  NAG_FREE(x);
  NAG_FREE(icom);
  NAG_FREE(iopts);
  NAG_FREE(needi);
  return exit_status;
}

static void display_integration_details(const Integer ni,
                                        const Integer iopts[],
                                        const double opts[],
                                        const Integer icom[],
                                        const double com[])
{
#define FCOUNT(J) icom[fcp + J-2]
#define EVALS(J,K) icom[evalsp +(K-1)*ldi + J-2]
#define SINFOI(L,K) icom[sinfoip +(K-1)*ldi + L-2]
#define DS(J,K) com[dsp + (K-1)*ldr + J-2]
#define ES(J,K) com[esp + (K-1)*ldr + J-2]
#define SINFOR(L,K) com[sinforp + (K-1)*ldr + L-2]
  double lbnd, ubnd, rvalue;
  Integer ldi, ldr, sinfoip, sinforp, evalsp, fcp, dsp, esp, nseg, nsdiv;
  Integer child1, child2, j, k, level, parent, sid, index, len_cvalue;
  char cvalue[17];
  NagError fail;
  Nag_VariableType optype;

  /* Request communication array indices */
  INIT_FAIL(fail);
  len_cvalue = 17;
  nag_quad_opt_get("Index nseg", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  nseg = icom[index - 1];
  nag_quad_opt_get("Index nsdiv", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  nsdiv = icom[index - 1];
  nag_quad_opt_get("Index ldi", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  ldi = icom[index - 1];
  nag_quad_opt_get("Index ldr", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  ldr = icom[index - 1];
  nag_quad_opt_get("Index fcp", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  fcp = icom[index - 1];
  nag_quad_opt_get("Index evalsp", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  evalsp = icom[index - 1];
  nag_quad_opt_get("Index sinfoip", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  sinfoip = icom[index - 1];
  nag_quad_opt_get("Index dsp", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  dsp = icom[index - 1];
  nag_quad_opt_get("Index esp", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  esp = icom[index - 1];
  nag_quad_opt_get("Index sinforp", &index, &rvalue, cvalue, len_cvalue,
                   &optype, iopts, opts, &fail);
  sinforp = icom[index - 1];

  printf("\n Information on integration:\n ni = %3" NAG_IFMT
         ", nseg = %3" NAG_IFMT ", nsdiv = %3" NAG_IFMT ".", ni, nseg, nsdiv);
  for (j = 1; j <= ni; j++)
    printf("\n Integral %2" NAG_IFMT " total approximations: %3" NAG_IFMT ".",
           j, FCOUNT(j));
  printf("\n\n Information on subdivision and evaluations over segments.\n");
  for (k = 1; k <= nseg; k++) {
    printf("\n");
    sid = SINFOI(1, k);
    parent = SINFOI(2, k);
    child1 = SINFOI(3, k);
    child2 = SINFOI(4, k);
    level = SINFOI(5, k);
    lbnd = SINFOR(1, k);
    ubnd = SINFOR(2, k);
    printf(" Segment %3" NAG_IFMT ".\n Sid = %3" NAG_IFMT ", Parent = "
           "%3" NAG_IFMT ", Level = %3" NAG_IFMT ".\n", k, sid, parent,
           level);
    if (child1 > 0)
      printf(" Children = (%3" NAG_IFMT ",%3" NAG_IFMT ").\n", child1,
             child2);
    printf(" Bounds (%11.4e,%11.4e).\n", lbnd, ubnd);
    for (j = 1; j <= ni; j++) {
      if (EVALS(j, k) > 0 && EVALS(j, k) < 5) {
        printf(" Integral %2" NAG_IFMT " approximation : %11.4e.\n", j,
               DS(j, k));
        printf(" Integral %2" NAG_IFMT " error estimate: %11.4e.\n", j,
               ES(j, k));
        if (EVALS(j, k) == 3)
          printf(" Integral %2" NAG_IFMT " evaluation has been superseded by "
                 "descendants.\n", j);
      }
    }
  }
  fflush(stdout);
}

static void display_option(const char *optstr, const Nag_VariableType optype,
                           const Integer ivalue, const double rvalue,
                           const char *cvalue)
{
  /* Query optype and print the appropriate option values. */
  switch (optype) {
  case Nag_Integer:
    printf("   %30s :    %13" NAG_IFMT "\n", optstr, ivalue);
    break;
  case Nag_Real:
    printf("   %30s :    %13.4e\n", optstr, rvalue);
    break;
  case Nag_Character:
    printf("   %30s : %16s\n", optstr, cvalue);
    break;
  case Nag_Integer_Additional:
    printf("   %30s :    %3" NAG_IFMT " %16s\n", optstr, ivalue, cvalue);
    break;
  case Nag_Real_Additional:
    printf("   %30s :    %13.4e %16s\n", optstr, rvalue, cvalue);
    break;
  default:;
  }
  fflush(stdout);
}