C++ 中带有 Wrapper 函数的 GSL 错误处理程序状态

Question

我在 c++ 中将这个 wrapper 函数用于我在其他解决方案中看到的 gsl

gsl_function_pp Fp( std::bind(&Class::member_function, &(*this),  std::placeholders::_1) );
gsl_function *F = static_cast<gsl_function*>(&Fp);

在 GSL 的文档中， https://www.gnu.org/software/gsl/manual/html_node/Error-Reporting-Examples.html#Error-Reporting-Examples

它说 int status = gsl_function..。 将给出错误代码。但是我不知道如何使用上面的包装函数获取状态。

******************************更新1

所以我用这个函数来解决 fzero 问题

double gsl_root(gsl_function *F, double x_lo, double x_hi) {
    const gsl_root_fsolver_type *T;
    gsl_root_fsolver *s;

    T = gsl_root_fsolver_brent;
    s = gsl_root_fsolver_alloc (T);

    gsl_root_fsolver_set (s, F, x_lo, x_hi);

    int status;
    do {
        status = gsl_root_fsolver_iterate(s);

        x_lo = gsl_root_fsolver_x_lower (s);
        x_hi = gsl_root_fsolver_x_upper (s);

        status = gsl_root_test_interval (x_lo, x_hi, 0, 1e-12);

    } while(status == GSL_CONTINUE);

    gsl_root_fsolver_free(s);

    return 0.5*(x_lo + x_hi);   
}

如何获得上述函数的状态输出，我尝试了迭代和区间函数但没有奏效，我已经不知道为什么和如何。

例如，当 GSl 没有解决方案时，我会收到如下错误：

gsl: brent.c:74: ERROR: endpoints do not straddle y=0
Default GSL error handler invoked.

所以这就是我真正想要的。我只想处理这种情况。很抱歉造成混淆，希望现在问题更清楚。

Answer 1

我认为我的评论解决了问题，但我会在这里总结一下以结束问题。 GSL 函数只是以下 c-struct 的 typedef

struct gsl_function_struct 
{
  double (* function) (double x, void * params);
  void * params;
};

typedef struct gsl_function_struct gsl_function ;

那么，很明显 gsl_function 不是一个例程，例如“int gsl_fft_complex_radix2_forward(...)”，它执行一些操作并返回一个带有错误状态的整数。

更新 1：我仍然不明白错误状态和包装器之间的关系是什么。所以让我展示一个几乎完整的例子（gsl集成）

 // use smart pointers in c++
 template< typename T > class deleter;

 template<> class deleter< gsl_integration_workspace > {
    public:
    void operator()( gsl_integration_workspace* ptr ) { 
        gsl_integration_workspace_free(ptr);
        return; 
    }
  };

  typedef std::unique_ptr< gsl_integration_workspace, 
   deleter< gsl_integration_workspace > > cpp_gsl_integration_workspace;


 // main program
 // TURN OFF GSL ERROR ERROR HANDLER (GSL JUST PRINT ERROR MESSAGE AND KILL THE PROGRAM IF FLAG IN ON)
 gsl_set_error_handler_off();

 auto ptr = [](double x)->double{ return x; };
 std::function<double(const double)> F1( std::cref(ptr) );
 gsl_function_pp F2(F1);
 gsl_function *F = static_cast<gsl_function*>(&F2); 

 double result;
 double error;
 double lower_limit = 0;
 double upper_limit = 1;
 double abs_eps = 0;
 double rel_eps= 1e3;
 int size = 1000;

 cpp_gsl_integration_workspace w( gsl_integration_workspace_alloc( size ) );

 int status = gsl_integration_qag ( F, lower_limit, upper_limit, abs_eps, rel_eps, 
    2000, GSL_INTEG_GAUSS15, w.get(), result, error );

 // check status and get error message if integration failed.
 // Nothing to do with the use of the C++ wrapper
 if ( status )
 {
   // gsl_strerror prints a message explaining the error 
   std::cout << "GSL FAIL: " << std::string( gsl_strerror (status) ) << std::endl;
   exit(1);
 }

 std::cout << result << "   " << error << std::endl

使用 fft 函数，您应该做一些非常相似的事情 int status = ....; if(status) {cout << std::string( gsl_strerror (status) ) << std::endl; exit(1); // or throw an exception}

更新 2：我希望现在清楚错误句柄与 C++ 包装器无关。

首先：当您设置gsl_set_error_handler_off() 时，您有责任检查 gsl 调用的错误。然后你必须将你当前的代码替换为这样的东西

int status;

// this will be the crucial line for my the answer to your problem, but it is important
// to stress all places where you must check the status of the gsl function calls.
status = gsl_root_fsolver_set (s, F, x_lo, x_hi);
if(status) std::cout << std::string( gsl_strerror (status) ) << std::endl;

do {
    status = gsl_root_fsolver_iterate(s);

    if(status) std::cout << std::string( gsl_strerror (status) ) << std::endl;

    x_lo = gsl_root_fsolver_x_lower (s);
    x_hi = gsl_root_fsolver_x_upper (s);

    status = gsl_root_test_interval (x_lo, x_hi, 0, 1e-12);

    if(status) std::cout << std::string( gsl_strerror (status) ) << std::endl;

} while(status == GSL_CONTINUE);

if(status != GSL_SUCCESS ) std::cout << std::string( gsl_strerror (status) ) << std::endl;

第二：如果您检查 GSL 源代码，您将看到该消息

gsl: brent.c:74: 错误：端点不跨越 y=0

写在函数上

 static int brent_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper)

此外，函数 set 在您的特定情况下只是指向 brent_init 的指针

typedef struct
{
  const char *name;
  size_t size;
  int (*set) (void *state, gsl_function * f, double * root, double x_lower, double x_upper);
  int (*iterate) (void *state, gsl_function * f, double * root, double * x_lower, double * x_upper);
}
gsl_root_fsolver_type;

static const gsl_root_fsolver_type brent_type =
{"brent",                               /* name */
 sizeof (brent_state_t),
 &brent_init,
 &brent_iterate};

然后，您无法检测到错误，因为您错过了以下调用的状态

 status = gsl_root_fsolver_set (s, F, x_lo, x_hi);
 if(status) std::cout << std::string( gsl_strerror (status) ) << std::endl; 

最后，当发生这种情况时，检查 GSL 源代码很重要。错误消息显示发生错误的文件名和行号。所以享受 GSL 是一个开源代码并检查那里发生了什么！此外，GSL 源代码组织良好、干净且易于理解。