Re: [Help-gsl] GSL 1.16 gsl_interp_accel_find not threadsafe?

[Jim Leek]

OK, that's the problem then.  Yes, I am using the same object with 
multiple threads.

I had assumed that once the object was constructed, it was safe to read 
it from multiple threads.  So a very common thing I have done is:

//This constructor is only run on 1 thread of course
GSL1DFunction::GSL1DFunction(OBJECT* obj) {

  ....[snip].....

  acc = gsl_interp_accel_alloc ();
  spline = gsl_spline_alloc (gsl_meth, XX.size());

  gsl_spline_init (spline, &XX[0], &YY[0], XX.size());
}

//This eval function is called by many threads
double GSL1DFunction::eval(double X) const {

  return gsl_spline_eval (spline, X, acc);
}

//In this case I was calling this integration function from multiple 
threads:
double GSL1DFunction::integ(double AA, double BB) const {
  return gsl_spline_eval_integ (spline, AA, BB, acc);
}


So, I've done this a lot.  This is the first time it has caused an 
issue, but I think that is just because the loop was unusually tight and 
the integration call unusually long.  It it correct that both 
gsl_spline_eval_integ AND gsl_spline_eval are not thread-safe in this 
manner?  It would be difficult for me to make an individual object for 
each thread with OpenMP.  I will probably have to use a mutex.

Thanks,
Jim

On 03/21/2016 12:45 PM, Patrick Alken wrote:
> Hard to say without any example code, but are you trying to use the 
> same accel object from all the threads? Each thread will need its own 
> interpolator and accel workspace. From the two lines you posted, I 
> don't see any global variables or anything that would cause problems 
> as long as each thread has its own accel workspace.
>
> Patrick
>
> On 03/21/2016 12:34 PM, Leek, Jim wrote:
> > I am using GSL in an openMP program.  It usually goes pretty well, 
> > but I recently ran a problem where I started getting random results 
> > (the deterministic problem returned wildly different results on each 
> > run.)
> >
> > I was finally able to get Intel Inspector to give me a clue.  It says 
> > there is a data race in gsl_intep_accel_find at:
> >
> > Line 210,211:
> >    a->miss_count++;
> >    a->cache - gsl_interp_bsearch(ca, x, x_intex, len-1)
> >
> > The call stack is:
> > gsl_spline_eval_integ:190
> > gsl_intep_eval_integ:273
> > cspline_eval_integ:403
> > gsl_interp_accel_find:211
> >
> > Avoiding calling this in the threading region does seem to solve the 
> > issue.
> >
> > 1) Is this a known issue?
> > 2) Is there a known way to get around it?
> > 3) Is it already fixed in a newer version of GSL?
> >
> > Thanks,
> > Jim