/* Integrate.c        Last Modified 10/07/98 */
/* Call Genz's DCUHRE to integrate a funcion */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

void main(int argc, char *argv[ ])
{
  int i;
  long ndim;      /* Number of variables in integrand(s) */
  long numfun;    /* Number of components in the vector integrand function */
  double *a;       /* Lower limit of integration */
  double *b;       /* Upper limit of integration */
  long minpts;    /* Minimum number of funsub calls */
  long maxpts;    /* Maximum number of funsub calls */
  double epsabs;  /* Requested absolute error */
  double epsrel;  /* Requested relative error */
  long key;       /* Selects integration rule */
  long nw;        /* Length of work array */
  long restar;    /* restar = 0, this is the first attempt
		      = 1, continuation of a previous attempt */
  double *result;  /* An array of the approximation to the integrals */
  double *abserr;  /* An array of estimates for absolute errors */
  long neval;     /* Number of function evaluations used */
  long ifail;     /* Exit condition indicator */
  double *work;    /* Work array */
  double *vecalloc(int nr);
  char   *charalloc(int nr);
  void func(long *nd, double *arg, long *nf, double *fv);
  void dcuhre_(long *ndim,long *numfun,double *a,double *b, 
	       long *minpts,long *maxpts,void *funsub,
	       double *epsabs,double *epsrel,long *key, 
	       long *nw,long *restar,double *result,
	       double *abserr,long *neval,long *ifail,
	       double *work);
  
  ndim=5;          /* Number of variables in integrand(s) */
  numfun=1;        /* Number of components in the vector integrand
		      function */
  a=vecalloc((int) ndim); /*allocate memory for ndim double vector*/
  b=vecalloc((int) ndim);
  for (i=0;i<ndim;i++){
    a[i]=0;          /* Lower limit of integration */
    b[i]=1.0;        /* Upper limit of integration */
  }
  minpts=0;        /* Minimum number of funsub calls */
  maxpts=10e6;     /* Maximum number of funsub calls */
  epsabs=0;        /* Requested absolute error */
  epsrel=10e-4;    /* Requested relative error */
  key=0;           /* Selects integration rule */
  nw=1000000;     /* Length of work array */
  restar=0;        /* restar = 0, this is the first attempt
		      = 1, continuation of a previous attempt */
  result=vecalloc((int) numfun);  /* resultant approximation to the integral(s) */
  abserr=vecalloc((int) numfun);  /* An array of estimates for
				     absolute errors */
  printf("Allocating %6.4eMB Memory\n",(((double)nw)*sizeof(double))/1000000);
  work=vecalloc((int) nw);        /* Work array */
  neval=0;
  dcuhre_(&ndim,&numfun,a,b,&minpts,&maxpts,func,
	  &epsabs,&epsrel,&key,&nw,&restar,result,abserr,
	  &neval,&ifail,work);
  printf("\n");
  printf("Integral=%6.4e\n",result[0]);
  printf(" w/ abs error = %6.4e\n",abserr[0]);
  printf(" %i function evaluations used\n",(int) neval);
  if (((int) ifail)==0)
    printf(" Note:  program exited with no errors.\n");
  else{
    printf(" ERROR:  program exited with error code %1i\n",((int) ifail));
  }
  printf("\n");
}