/* nag_quad_1d_gauss_vec (d01uac) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */
#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("%5ld   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;
    }
}