1

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).

Drade
  • 157
  • 10

2 Answers2

2

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);
mhrsalehi
  • 1,124
  • 1
  • 12
  • 29
0

All you need to do is to set your own error handler via

https://www.gnu.org/software/gsl/doc/html/err.html#c.gsl_set_error_handler

See e.g. Callbacks with ctypes (How to call a python function from C) how to declare an appropriate callback and pass it.

The default handler aborts (as lined out in the documentation link above), but it's your prerogative to not do that of course :)

deets
  • 6,285
  • 29
  • 28