/* nag_quad_1d_gauss_vec (d01uac) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL f(const double x[], const Integer nx, double fv[],
                         Integer *iflag, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  Integer exit_status = 0;
  double a, b, dinest;
  Integer funid, i, nstor;
  /* Arrays */
  Integer iuser[1];
  /* Nag Types */
  Nag_Comm comm;
  NagError fail;
  Nag_QuadType quad_type;

  printf("nag_quad_1d_gauss_vec (d01uac) Example Program Results\n");

  INIT_FAIL(fail);

  /* Use comm to pass information to f. */
  comm.iuser = iuser;

  for (funid = 1; funid <= 6; funid++) {
    switch (funid) {
    case 1:
      {
        printf("\nGauss-Legendre example\n");
        a = 0.0;
        b = 1.0;
        quad_type = Nag_Quad_Gauss_Legendre;
        break;
      }
    case 2:
      {
        printf("\nRational Gauss example\n");
        a = 2.0;
        b = 0.0;
        quad_type = Nag_Quad_Gauss_Rational_Adjusted;
        break;
      }
    case 3:
      {
        printf("\nGauss-Laguerre example (adjusted weights)\n");
        a = 2.0;
        b = 1.0;
        quad_type = Nag_Quad_Gauss_Laguerre_Adjusted;
        break;
      }
    case 4:
      {
        printf("\nGauss-Laguerre example (normal weights)\n");
        a = 2.0;
        b = 1.0;
        quad_type = Nag_Quad_Gauss_Laguerre;
        break;
      }
    case 5:
      {
        printf("\nGauss-Hermite example (adjusted weights)\n");
        a = -1.0;
        b = 3.0;
        quad_type = Nag_Quad_Gauss_Hermite_Adjusted;
        break;
      }
    case 6:
      {
        printf("\nGauss-Hermite example (normal weights)\n");
        a = -1.0;
        b = 3.0;
        quad_type = Nag_Quad_Gauss_Hermite;
        break;
      }
    }

    iuser[0] = funid;
    for (i = 0; i < 6; i++) {
      nstor = pow(2, i + 1);
      /* Compute the one-dimensional integral employing Gaussian quadrature,
       * with quadrature type and weights specified in quad_type, using
       * nag_quad_1d_gauss_vec (d01uac).
       */
      nag_quad_1d_gauss_vec(quad_type, a, b, nstor, f, &dinest, &comm, &fail);
      switch (fail.code) {
      case NE_NOERROR:
      case NE_QUAD_GAUSS_NPTS_RULE:
      case NE_UNDERFLOW:
      case NE_WEIGHT_ZERO:
        {
          /* The definite integral has been estimated. */
          printf("%5" NAG_IFMT "   Points  Answer = %10.5f\n", nstor, dinest);
          break;
        }
      default:
        {
          /* A solution could not be calculated due to an illegal parameter
           * or a requested exit.
           */
          printf("%s\n", fail.message);
          exit_status++;
        }
      }
    }
  }
  return exit_status;
}

static void NAG_CALL f(const double x[], const Integer nx, double fv[],
                       Integer *iflag, Nag_Comm *comm)
{
  Integer i, funid;

  funid = comm->iuser[0];
  switch (funid) {
  case 1:
    for (i = 0; i < nx; i++)
      fv[i] = 4.0 / (1.0 + x[i] * x[i]);
    break;
  case 2:
    for (i = 0; i < nx; i++)
      fv[i] = 1.0 / (x[i] * x[i] * log(x[i]));
    break;
  case 3:
    for (i = 0; i < nx; i++)
      fv[i] = exp(-x[i]) / x[i];
    break;
  case 4:
    for (i = 0; i < nx; i++)
      fv[i] = 1.0 / x[i];
    break;
  case 5:
    for (i = 0; i < nx; i++)
      fv[i] = exp(-3.0 * x[i] * x[i] - 4.0 * x[i] - 1.0);
    break;
  case 6:
    for (i = 0; i < nx; i++)
      fv[i] = exp(2.0 * x[i] + 2.0);
    break;
  default:
    *iflag = -1;
  }
}