How to handle error in integration in GSL

Question

I am writing a python code that needs to calculate a lot of integrals really fast, so instead of using scipy.integrate using a c function, I am using ctypes to calculate all of my integrals in c.

I am importing the c function with this code:

ctypes.CDLL('/usr/lib/i386-linux-gnu/libgslcblas.so', mode=ctypes.RTLD_GLOBAL)
ctypes.CDLL('/usr/lib/i386-linux-gnu/libgsl.so', mode=ctypes.RTLD_GLOBAL)
lib = ctypes.CDLL('/home/aurelien/Desktop/Project/within_opt_model_1/integral_test.so')

In c, I am calculating the integral using:

gsl_integration_workspace *w = gsl_integration_workspace_alloc(200);
double result, err;

gsl_function F;
F.function = &integral;
F.params = &params;

gsl_integration_qags(&F, z, 1, 1e-6, 1e-6, 200, w, &result, &err);
gsl_integration_workspace_free(w);

The problem is that sometimes (~5% of the time), since my integral is slowly convergent, it will kill my program and print:

gsl: qags.c:563: ERROR: interval is divergent, or slowly convergent
Default GSL error handler invoked
Aborted (core dumped)

This problem does not occur if I put 1e-10 instead of 1e-6, but it makes the integral calculation go a lot slower.

I would like a way to do something like a try, catch so that most of the time it uses 1e-6 and goes fast, and when it fails it uses 1e-10.

I have tried doing:

int status = gsl_integration_qags(&F, z, 1, 1e-6, 1e-6, 200, w, &result, &err);
printf("%i\n", status);

But this only prints 0, since the error aborts before returning any value.

My guess is that I need create my own error handler method, but I don't know how to do that.

Thank you very much! (I can show you more code if that's helpful, I tried to keep everything concise).

Answer 1

from GCC GSL :

Error Handlers

The default behavior of the GSL error handler is to print a short message and call abort(). When this default is in use programs will stop with a core-dump whenever a library routine reports an error. This is intended as a fail-safe default for programs which do not check the return status of library routines (we don’t encourage you to write programs this way).

you can call gsl_set_error_handler_off() to disable defualt error handeler (This function turns off the error handler by defining an error handler which does nothing).

I think the following will do the job

gsl_set_error_handler_off()
gsl_integration_workspace *w = gsl_integration_workspace_alloc(200);
double result, err;

gsl_function F;
F.function = &integral;
F.params = &params;

int status = gsl_integration_qags(&F, z, 1, 1e-6, 1e-6, 200, w, &result, &err);
if(status == GSL_EDIVERGE){
status = gsl_integration_qags(&F, z, 1, 1e-10, 1e-10, 200, w, &result, &err);

/*handle other errors here...*/
}
gsl_integration_workspace_free(w);