Hello,
I converted your main.cpp to C (I couldn't make it compile with my
g++). Also I linked it against GSL v2.5, but the CQUAD code hasn't
changed since 2.3 as far as I know. The program works as expected (C
code attached).
$ ./main
result = 4.390157623791e-164
abserr = 0.000000000000e+00
nevals = 33
Mathematica agrees with the result to about 11 significant digits.
I am unable to diagnose problems with fgsl since I have never used it.
Patrick
On 11/20/18 9:36 AM, Rico Buchholz wrote:
>
> Hello!
>
> I am working with a fortran project which uses fgsl (v1.2.0, together
> with gsl v2.3.0).
>
> The code reports an arithmetic exception in an integration routine,
> when trying to integrate exp(-x^2) from x=19.3 to 20.0.
> I would expect the code to work without problems (maybe underflow or
> denormal).
> Trying to create a minimal working example i came up with the attached
> file (the cmake input i used for compilation is also attached, sorry I
> could not boil this further down).
> As a corresponding version for c++ works (also attached), this might
> relate to a bug in fgsl, but for me (a person without knowledge of the
> internals), that seems questionable according to what causes the
> problem and what not.
>
>
> The output is as follows:
>
> Program received signal SIGFPE: Floating-point exception - erroneous
> arithmetic operation.
>
> Backtrace for this error:
> #0 0x7fa39acc5d1d in ???
> #1 0x7fa39acc4f7d in ???
> #2 0x7fa39a1e405f in ???
> #3 0x560ac54c9eff in ???
> #4 0x560ac54860cf in neo2
> at
> /afs/itp.tugraz.at/user/buchholz/Programs/neo2-minimal/NEO-2-QL/neo2.f90:31
> #5 0x560ac54861d8 in main
> at
> /afs/itp.tugraz.at/user/buchholz/Programs/neo2-minimal/NEO-2-QL/neo2.f90:4
>
> It does not occur with the alternative integration routine commented
> in the code.
> This does not occur, if one reduces the lower bound to 19.2 or smaller
> (which is a reason why I suspect that this might be a gsl bug - a bug
> in the interface should have no influence on this).
>
>
> With kind regards
> Rico Buchholz
>
>
> PS: I did also send the bug report to fgsl.
>
/**
*
* g++ main.cpp -o main.x -L`gsl-config --libs`
*
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
double f(double x, void * params) {
return exp(-x*x);
}
int main (void)
{
size_t const n = 100;
double const a = 19.3;
double const b = 20.0;
double const epsabs = 1.0e-13;
double const epsrel = 1.0e-13;
double result = 0.0;
double abserr = 0.0;
size_t nevals = 0;
gsl_integration_cquad_workspace * workspace;
gsl_function F;
F.function = &f;
F.params = NULL;
workspace = gsl_integration_cquad_workspace_alloc(n);
gsl_integration_cquad(&F, a, b, epsabs, epsrel, workspace, &result, &abserr, &nevals);
fprintf(stderr, "result = %.12e\n", result);
fprintf(stderr, "abserr = %.12e\n", abserr);
fprintf(stderr, "nevals = %zu\n", nevals);
gsl_integration_cquad_workspace_free(workspace);
return 0;
}
