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

[David Zaslavsky]

That's true, you don't have to pass in the accel object. But if you do
that, you wouldn't get the benefit of acceleration, which is a speedup
when you evaluate the spline at many nearby points in a row. (The accel
object basically makes the search for the evaluation point start from
the last point at which you evaluated the spline, instead of starting
from the middle of the interpolation range.)

If you want to get the benefit of acceleration, you can have one
globally-shared spline but a separate accel object for each thread. The
point Patrick was making was that the spline can be shared between
threads (after it gets initialized), and you know that because the
gsl_interp_spline_eval() functions take a const gsl_spline pointer.
Since the pointer is const, you know that the function will not change
anything in the spline object itself, and so it's safe to use in
multiple threads. But the evaluation functions *do* change something in
the accel object, so you can't share accel objects between threads. You
can tell because the evaluation functions take a non-const pointer to
the accel object.

:) David

On 3/23/16 12:32 AM, Leek, Jim wrote:
> OK, it appears that I don't actually need to pass in the accel object?  
> Because the gsl_spline* itself is const, it's only the accel object that 
> isn't.  And digging around a bit, I see code like this:
> 
>  if (acc != 0)
>     {
>       index_a = gsl_interp_accel_find (acc, x_array, size, a);
>       index_b = gsl_interp_accel_find (acc, x_array, size, b);
>     }
>   else
>     {
>       index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
>       index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
>     }
> 
> So, it looks like I can just pass in NULL as the accel object, and I'll be 
> able to use the same gsl_spline object with multiple threads, and GSL will 
> still produce that correct value (if a little more slowly.)  
> 
> That's correct? 
> Thanks,
> Jim
>   
> ________________________________________
> From: Patrick Alken address@hidden
> Sent: Monday, March 21, 2016 1:47 PM
> To: Leek, Jim; address@hidden
> Subject: Re: [Help-gsl] GSL 1.16 gsl_interp_accel_find not threadsafe?
> 
> Yes these functions are thread-safe in the sense that they do not set
> any global variables. So if you have 2 threads, each with its own accel,
> then each thread will update its own accel object just fine, without
> causing any problems due to global variables being set.
> 
> However, as you noted, when you interpolate a point, internal variables
> inside the accel object are updated, so if 2 threads are accessing the
> same accel object, they may try to update these internal variables at
> the same time, causing a data race condition.
> 
> I would say its only safe to access the same gsl object from different
> threads when you only need to pass a const * pointer to the object -
> this means that no internal state variables are updated by the function
> call, so read-only routines should be safe to call from any number of
> threads.
> 
> But for routines which update the internal state, each thread needs its
> own copy of the workspace.
> 
> Patrick
> 
> On 03/21/2016 02:37 PM, Jim Leek wrote:
>> 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
>>>
>>>
>>
> 
>