Reputation: 11
gsl_integration_qag works with 1 core (with/without openMP), but fails with multi-threads (i.e. >1).
Some information that may help...
gsl-2.5
#define _OPENMP 201107
Depending on the number of cores, I can get error reports of:
gsl: qag.c:248: ERROR: roundoff error prevents tolerance from being achieved (comment: usually with a small number of cores)
gsl: qag.c:257: ERROR: maximum number of subdivisions reached (comment: usually with a large number of cores)
A large max iteration number given to gsl_integration_qag only delays the code to crash.
The integration function is (can be more specific if needed):
double Func(double Param1, ..., double ParamN){
double result, error;
gsl_function F;
gsl_integration_workspace * w
= gsl_integration_workspace_alloc (1000);
struct parameters_gsl_int_ parameters_gsl = {
.Param1 = Param1,
...
.ParamN = ParamN,};
F.function = &func_integrand;
F.params = ¶meters_gsl;
gsl_integration_qag (&F, LOWER_LIMIT, UPPER_LIMIT, 0, 0.001,
1000, GSL_INTEG_GAUSS61, w, &result, &error);
gsl_integration_workspace_free (w);
return result;
}
The OpenMP part that calls the integration is:
void call_Func(int Nbin, double array[], double Param1[], double Param2, ... double ParamN){
int i;
...
#pragma omp parallel shared(Nbin, array, Param1, ..., ParamN) private(i)
{
#pragma omp for
for (i=0; i<Nbin; i++)
array[i] = Func(Param1[i], Param2, ..., ParamN);
}
...
}
I'm new to both GSL and openMP. I hope I am using gsl_integration_qag correctly and the definition of shared or private variables makes sense.
btw, it's the same question as this 2014 one (gsl openmp failed integration), but I couldn't find the solution in this post.
Upvotes: 0
Views: 954
Reputation: 11
Problem solved...
It is actually due to func_integrand having also a term which is estimated using gsl_integration_qag. There were some global variables adopted in this calculation, which I didn't capture before.
Upvotes: 1