如何处理 GSL 中的集成错误

Question

我正在编写 python 代码，需要非常快地计算大量积分，所以我没有使用 scipy.integrate 使用 c 函数，而是使用 ctypes 来计算我在 c.

中的所有积分

我正在使用以下代码导入 c 函数：

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

在 c 中，我使用以下方法计算积分：

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

问题是有时（约 5% 的时间），由于我的积分收敛缓慢，它会终止我的程序并打印：

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

如果我用1e-10代替1e-6，就不会出现这个问题，但它会使积分计算速度变慢。

我想要一种方法来做类似 try, catch 的事情，这样大多数时候它使用 1e-6 并且运行速度很快，当它失败时它使用 1e-10。

我试过：

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

但这只打印 0，因为错误在返回任何值之前中止。

我的猜测是我需要创建自己的错误处理程序方法，但我不知道该怎么做。

非常感谢！（如果有帮助，我可以向您展示更多代码，我尽量让一切都简洁明了）。

Answer 1

您需要做的就是通过

设置您自己的错误处理程序

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

参见例如如何声明适当的回调并传递它。

默认处理程序中止（如上文 link 中所述），但您当然有权不这样做:)

Answer 2

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

您可以调用 gsl_set_error_handler_off() 来禁用默认错误处理程序（此函数通过定义一个不执行任何操作的错误处理程序来关闭错误处理程序）。

我认为以下内容可以完成工作

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

如何处理 GSL 中的集成错误

How to handle error in integration in GSL

c

python

integration

ctypes

gsl