I am the author of CQUAD in the GSL. The interface is almost identical to that of QAGS, so if you've used the latter, it should not be difficult at all to try the former. Just remember not to convert your NaNs and Infs to zeros in the integrand -- the code will deal with these itself.
The routine is also available in Octave as quadcc, and in Matlab here.
Could you provide an example of the integrands you are dealing with?
Update
Here's an example of using CQUAD to integrate a function with a singularity at one of the endpoints:
#include <stdio.h>
#include <gsl/gsl_integration.h>
/* Our test integrand. */
double thefunction ( double x , void *param ) {
return sin(x) / x;
}
/* Driver function. */
int main ( int argc , char *argv[] ) {
gsl_function f;
gsl_integration_cquad_workspace *ws = NULL;
double res, abserr;
size_t neval;
/* Prepare the function. */
f.function = &thefunction;
f.params = NULL;
/* Initialize the workspace. */
if ( ( ws = gsl_integration_cquad_workspace_alloc( 200 ) ) == NULL ) {
printf( "main: call to gsl_integration_cquad_workspace_alloc failed.\n" );
abort();
}
/* Call the integrator. */
if ( gsl_integration_cquad( &f, 0.0 , 1.0 , 1.0e-10 , 1.0e-10 , ws , &res , &abserr , &neval ) != 0 ) {
printf( "main: call to gsl_integration_cquad failed.\n" );
abort();
}
/* Print the result. */
printf( "main: int of sin(x)/x in [0,1] is %.16e +/- %e (%i evals).\n" ,
res , abserr , neval );
/* Free the workspace. */
gsl_integration_cquad_workspace_free( ws );
/* Bye. */
return 0;
}
which I compiled with gcc -g -Wall cquad_test.c -lgsl -lcblas. The output is
main: int of sin(x)/x in [0,1] is 9.4608307036718275e-01 +/- 4.263988e-13 (63 evals).
Which, given the result computed in Maple to 20 digits, $0.94608307036718301494$, is correct to 14 digits.
Note that there is nothing special here, neither to tell CQUAD where the singularity is, or any special treatment within the integrand itself. I just let it return NaNs, and the integrator takes care of them automatically.
Note also that there is a bug in the latest GSL version 1.15 which can affect the treatment of singularities. It has been fixed, but has not made it to the official distribution. I used the most recent source, downloaded with bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.
