/* nag_multid_quad_monte_carlo (d01gbc) Example Program.

 *

 * Copyright 1991 Numerical Algorithms Group.

 *

 * Mark 2, 1991.

 * Mark 6 revised, 2000.

 * Mark 7 revised, 2001.

 * Mark 8 revised, 2004.

 *

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <math.h>

#include <nagd01.h>



#ifdef __cplusplus

extern "C" {

#endif

  static double f(Integer  ndim, double x[]);

#ifdef __cplusplus

}

#endif



#define MAXCLS 20000



int main(void)

{



  Integer exit_status=0, k, maxcls=MAXCLS, mincls, ndim=4;

  NagError fail;

  Nag_MCMethod method;

  Nag_Start cont;

  double *a=0, acc, *b=0, *comm_arr=0, eps, finest;



  INIT_FAIL(fail);

  Vprintf("nag_multid_quad_monte_carlo (d01gbc) Example Program Results\n");

  if (ndim>=1)

    {

      if ( !( a = NAG_ALLOC(ndim, double)) ||

           !( b = NAG_ALLOC(ndim, double)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

    }

  else

    {

      Vprintf("Invalid ndim.\n");

      exit_status = 1;

    }

  for (k=0; k<ndim; ++k)

    {

      a[k] = 0.0;

      b[k] = 1.0;

    }

  eps = 0.01;

  mincls = 1000;

  method = Nag_ManyIterations;

  cont = Nag_Cold;

  /* nag_multid_quad_monte_carlo (d01gbc).

   * Multi-dimensional quadrature, using Monte Carlo method

   */

  nag_multid_quad_monte_carlo(ndim, f, method, cont, a, b, &mincls, maxcls, eps,

                              &finest, &acc, &comm_arr, &fail);

  if (fail.code == NE_NOERROR || fail.code == NE_QUAD_MAX_INTEGRAND_EVAL)

    {

      if (fail.code == NE_QUAD_MAX_INTEGRAND_EVAL)

        {

          Vprintf("Error from nag_multid_quad_monte_carlo (d01gbc).\n%s\n",

                  fail.message);

          exit_status = 2;

        }

      Vprintf("Requested accuracy    = %10.2e\n",eps);

      Vprintf("Estimated value       = %10.5f\n", finest);

      Vprintf("Estimated accuracy    = %10.2e\n", acc);

      Vprintf("Number of evaluations = %5ld\n", mincls);

    }

  else

    {

      Vprintf("Error from nag_multid_quad_monte_carlo (d01gbc).\n%s\n",

              fail.message);

      Vprintf("%s\n", fail.message);

      exit_status = 1;

    }

 END:

  if (a) NAG_FREE(a);

  if (b) NAG_FREE(b);

  /* Free memory allocated internally */

  if (comm_arr) NAG_FREE(comm_arr);

  return exit_status;

}

static double f(Integer ndim, double x[])

{

  return x[0]*4.0*(x[2]*x[2])*exp(x[0]*2.0*x[2])/

    ((x[1]+1.0+x[3])*(x[1]+1.0+x[3]));

}