GSL: Error reporting

I want to use the GSL for integration http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html However, I find no convenient way how the integrated function

(the function f in the example http://www.gnu.org/software/gsl/manual/html_node/Numerical-integration-examples.html)

can report an error to the integrator. I want to integrate a function which itself results from an integration that could fail. This is my sample program

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errn开发者_JAVA技巧o.h>

double f (double x, void * params) {
    GSL_ERROR("test error",GSL_FAILURE);
    return 0.0;
}



int main (void)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

double result, error;

gsl_function F;
F.function = &f;

gsl_set_error_handler_off();
int status = gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                     w, &result, &error); 

printf ("status          = %d\n", status);
status  = GSL_FAILURE;
printf ("status          = %d\n", status);


gsl_integration_workspace_free (w);

return 0;
}

resulting in the output status = 0 status = -1

I think the integrator should rather stop and return my error code. How can I achieve this?

Thank you very much for your help!!!

2011-04-27: I also tried this variant, after Brian Gough told me,

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

double f (double x, void * params) {
    GSL_ERROR("test error",GSL_FAILURE);
    return GSL_NAN;
}



int main (void)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

double result, error;

gsl_function F;
F.function = &f;

gsl_set_error_handler_off();
int status = gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                     w, &result, &error); 

printf ("status          = %d\n", status);
status  = GSL_FAILURE;
printf ("status          = %d\n", status);


gsl_integration_workspace_free (w);

return 0;
}

it did not help either. I will now fill out a bug report.

Thanks to Xuebin Wu from the GSL Mailing list the problem is solved:

Hi,

GSL_ERROR itself is a macro, it looks like

      gsl_error (reason, __FILE__, __LINE__, gsl_errno);
      return gsl_errno;

The function already returns before you return NAN, because GSL_ERROR has been called. Turning the handler off just let the first line do nothing. The default error handler abort the program after printing error message.

I do not think it is a bug. Maybe you can write your own error handler to solve your problem. For example, you can use "goto" to jump out of gsl_integration_qags, or set some global variable to indicate the integration result is incorrect.

PS: I believe this macro is what you need,

Macro: GSL_ERROR_VAL (reason, gsl_errno, value) This macro is the same as GSL_ERROR but returns a user-defined value of value instead of an error code. It can be used for mathematical functions that return a floating point value.

The following example shows how to return a NaN at a mathematical singularity using the GSL_ERROR_VAL macro,

 if (x == 0)
   {
     GSL_ERROR_VAL("argument lies on singularity",
                   GSL_ERANGE, GSL_NAN);
   }

So I adjusted the code according to

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

double f (double x, void * params) {
//  return GSL_NAN;
    GSL_ERROR_VAL ("argument lies on singularity", GSL_ERANGE, GSL_NAN);
}



int main (void)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

double result, error;

gsl_function F;
F.function = &f;

gsl_set_error_handler_off();
int status = gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
                     w, &result, &error); 

printf ("status          = %d\n", status);
status  = GSL_FAILURE;
printf ("status          = %d\n", status);


gsl_integration_workspace_free (w);

return 0;
}

and everything works as expected...

A bit hackish, but I'd probably have your function store some flag. When it encounters an error it sets the flag and returns zero for all subsequent evaluations. Then, after you've integrated it you can check this flag to see if the result is valid.

What about to write a wrapper for the function which returns pointer to a structure, containing function results and error status ? Or if you use c++, this encapsulation can be made with use of objects ....