Cuba.jl icon indicating copy to clipboard operation
Cuba.jl copied to clipboard

suave not handling large maxevals

Open agerlach opened this issue 4 years ago • 2 comments

Hello, the docs state:

nvec, minevals, maxevals, nnew, nmin, neval are passed to the Cuba library as 64-bit integers, so they are limited to be at most typemax(Int64)

However, suave() seems to have issues returning NaN for large maxevals < typemax(Int64). I have not observed this issue with the other integrators in Cuba.jl

MWE:

suave((x,f)->f.=3.14, 1, 1, 
        rtol=1e-3, atol=1e-3, 
        nvec=1, maxevals = typemax(Int64)*.499 #works
suave((x,f)->f.=3.14, 1, 1, 
        rtol=1e-3, atol=1e-3, 
        nvec=1, maxevals = typemax(Int64)*.5) # returns NaN or random small number

agerlach avatar Jun 01 '20 21:06 agerlach

I believe this is an upstream issue and I can't do anything about it. When running the following program:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "cuba.h"

#define NDIM 1
#define NCOMP 1
#define USERDATA NULL
#define NVEC 1
#define EPSREL 1e-8
#define EPSABS 1e-8
#define VERBOSE 0
#define LAST 4
#define SEED 0
#define MINEVAL 0
#define MAXEVAL 4611686018427387904

#define STATEFILE NULL
#define SPIN NULL

#define NNEW 1000
#define NMIN 2

static int Integrand(const int *ndim, const double xx[],
		     const int *ncomp, double ff[], void *userdata) {
  ff[0]  = 3.14;
  return 0;
}

int main() {
  int comp, nregions, fail;
  long long int neval;
  double integral[NCOMP], error[NCOMP], prob[NCOMP];
  float start, end;
  llSuave(NDIM, NCOMP, Integrand, USERDATA, NVEC,
	  EPSREL, EPSABS, VERBOSE | LAST, SEED,
	  MINEVAL, MAXEVAL, NNEW, NMIN, 25.0,
	  STATEFILE, SPIN,
	  &nregions, &neval, &fail, integral, error, prob);
  printf("Result: %g\n", integral[0]);
  return 0;
}

I get

$ ./test-cuba
zsh: segmentation fault (core dumped)  ./test-cuba

Instead this works fine with

#define MAXEVAL 4611686018427387903

Please, report the issue to the author of the Cuba library, you can find his email in his website.

giordano avatar Jun 01 '20 22:06 giordano

Thanks for looking into this. I will send something to the author. In the mean time, maybe an update to the docs is in order since they explicitly state :

limited to be at most typemax(Int64)

agerlach avatar Jun 02 '20 13:06 agerlach