GSL Error Handler status with Wrapper function in C++

Question

I'm using this wrapper function in c++ for gsl that I saw on other solutions

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

In documentation of GSL, https://www.gnu.org/software/gsl/manual/html_node/Error-Reporting-Examples.html#Error-Reporting-Examples

It say int status = gsl_function... will give the error codes. However i couldnt figure out how to get status with wrapper function above.

**************************UPDATE 1

So I'm using this function for solving fzero problems

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

How can I get the status output for above function, I tried by iterate and interval functions but didnt work, already I have no idea and would like learn why & how.

For example when GSl has no solution I get error such as:

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

So that's what I want actually. I just want to handle this cases. Sorry for confusion, hope it's more more clear question now.

Hello. First: I have written some answers where I use this wrapper but I am not the original creator of it. This is a relatively well known solution. So, wrapper from "Vinicius Miranda" is not fair (also not helpful as citation). Can you please just say "wrapper that I saw on those solutions" and actually specify where you saw it and how you use it in your code.. — Vivian Miranda, Jun 10 '14 at 18:18
Second: There are functions in gsl that return error status like, for example: "int gsl_integration_qag(...)". gsl_function is not one of them. Actually gsl_function is just a c-struct that holds a function and a void pointer. See gsl_math.h in source code lines 123-130 — Vivian Miranda, Jun 10 '14 at 18:21
:( no i didnt understand or asked in awrong way.. Now let me tell again; in my case above, I'm calling a member function with wrapper.. Which variable will hold the status basically? it's okay gsl_func is a struct but how can i use that gsl_inteegration_qag() in my case? — questioner, Jul 01 '14 at 21:09

Vivian Miranda · Accepted Answer · 2014-07-02T16:47:48.887

I think my comment solves the problem, but I will summarize it here to close the question. GSL function is just a typedef to the following c-struct

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

typedef struct gsl_function_struct gsl_function ;

Then, it is clear that gsl_function is not a routine, like "int gsl_fft_complex_radix2_forward(...)" for example, that performs some operation and returns an integer with the error status.

Update 1: I still don't understand what is the relation between error status and the wrapper. So let me show a nearly complete example (gsl integration)

 // 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

With the fft functions you should do something very similar int status = ....; if(status) {cout << std::string( gsl_strerror (status) ) << std::endl; exit(1); // or throw an exception}

Update 2: I hope it is now clear that error handle has nothing to do with the C++ wrapper.

First: when you set gsl_set_error_handler_off() you are the responsible to check the error on gsl calls. Then you must replace your current code to something like this

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;

Second: if you check the GSL source code you will see that the message

gsl: brent.c:74: ERROR: endpoints do not straddle y=0

is written on function

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

Furthermore, the function set is just a pointer to brent_init in your particular case

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

Then, you couldn't detect the error because you missed the status on the following call

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

Last, it is important to check the GSL source code when something like this happens. The error message shows the file name and line number where the error occurs. So enjoy that GSL is an open source code and check what is happening there! Besides, GSL source code is very well organized, clean and easy to understand.

thanks for your efforts vinicius but so how can i use gsl_fft_complex_radix2_forward(...) with my case, with wrapper while calling member function..?? — questioner, Jul 01 '14 at 21:11
cool, so definitely now my question is more clear(also for me, sorry about that.. ) So, now I use something like::::---> 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); — questioner, Jul 02 '14 at 13:02
so here i tried to output status with gsl_strerror but it didnt come out. so what do u suggest for getting gsl errors correctly? — questioner, Jul 02 '14 at 13:08
aha got an output as GSL FAIL: the iteration has not converged yet. This is after gsl_root_test_interval but this is not what i want. because i need output when ""gsl: brent.c:74: ERROR: endpoints do not straddle y=0 Default GSL error handler invoked. "" — questioner, Jul 02 '14 at 13:18
PERFECT VINICIUS! Okay that status = ..._set was really important, i wouldnt figure out without u. and You re right it is not about wrapper, but i was really confused a lot, I'm trying to continue on a project written by another engineer, and I had no sleep around these days :/ sorry to getting too much time.+ Anywayy I learned a lot now over all those points. Hopefully I can help to other guys. — questioner, Jul 02 '14 at 17:11
!!! !!! !!Just a quick question: now i see 2 messages: ""the iteration has not converged yet"" and ""invalid argument supplied by user"" normally with defaul GSL error handler, it was tellioing that start is bigger than end, or endpoints do not straddle, so it was more detailed but now it seems more abstract. If u have any idea about that also i''ll bbe gladd. Nevertheless if there is no way, never mind, just asking. Thanks again thanks ..Gracias. campeón de brasil :))) — questioner, Jul 02 '14 at 17:15
Once you have one error, you should not proceed (just exit or throw an exception...I haven't done that in my example but you should because the second error is jut consequence of the first). The problem seems that you set 2 boundaries [xlower, xhigher] but f(slower) and f(higher) are both > 0 (or < 0). There is no way the root finder can converge in that situation because there is no root in your interval (or there are 2 (or even number of roots)..in any case - no convergence is possible)! — Vivian Miranda, Jul 02 '14 at 18:40
Yess so how can I get some error outputs like ""endpoints do not straddle y=0 Default GSL error handler invoked"" ?? because now I get hundreds of outputs like iteration has not converged yet. But I need to see only the cases when there is a certain error — questioner, Jul 03 '14 at 13:38
so can you clarify that "Once you have one error, you should not proceed " . Which cases are error? I see that status is 0, -2 or 2, and i guess 4 and may be some others.. so -2 is not an error it is continue, in which cases I should exit? or break? there re a lot of flags :/ — questioner, Jul 03 '14 at 15:55
if (status) returns true means it is an error so you should exit or throw an exception. Again, you should go to the source code. /gsl/roots/brent.c line 72 shows that this error message occurs (f_lower < 0.0 && f_upper < 0.0) or (f_lower > 0.0 && f_upper > 0.0)...this was exactly what I told you before. The function values at the end point are both positive or negative and this implies no root or an even number of roots. — Vivian Miranda, Jul 03 '14 at 16:51

GSL Error Handler status with Wrapper function in C++

1 Answers1