Yuxiang Qin
Yuxiang Qin

Reputation: 11

gsl_integration_qag failed with gsl openmp

gsl_integration_qag works with 1 core (with/without openMP), but fails with multi-threads (i.e. >1).

Some information that may help...

  1. gsl-2.5

  2. #define _OPENMP 201107

  3. 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)
    
  4. A large max iteration number given to gsl_integration_qag only delays the code to crash.

  5. 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 = &parameters_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;
    }
    
  6. 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

Answers (1)

Yuxiang Qin
Yuxiang Qin

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

Related Questions