Thread safety

[aafkegros]

Currently the SHTOOLS SHExpandGLQ function i regularly use fails in multithreaded applications. I think this comes from the modulation of fftw3.3.3 plans as only the fftw execute calls are thread-safe: https://www.fftw.org/fftw3_doc/Thread-safety.html
This probably affects (almost) all fftw calls in SHTOOLS.

Would it be an idea to implement the brute fftw_make_planner_thread_safe call and force fft expansions to be thread safe?

[MarkWieczorek]

I have used the fortran routines with openmp in the past without any problems. Could you give me more info about when this fails and when it doesn't? Or does it fail for even the most simply tasks?

[mreineck]

If the problem is indeed FFTW-related, it's most likely not connected with ducc, since ducc doesn't use FFTW internally.

[aafkegros]

Oh interesting! I had it fail when multiple threads were trying to execute SHExpandGLQ calls on different datasets (with the same dimensions) at the same time. The failure was then a segmentation fault without clear error codes. I found the FFTW thread safety information, and saw that the variables they talk about there were edited in the SHExpandGLQ code.
I tried now changing backend to ducc as @mreineck suggested, but my segmentation faults stay.

The threading implementation i'm using is a python-rewrite of Java threading that is implemented in the bigger software i'm writing in, so this may also have specific vulnerabilities?

I managed to fix it in my code by encapsulating the SHExpandGLQ in an RLock - so it's not very pressing for me, but weird that it happens especially if it's stable with openMP

[MarkWieczorek]

I don't use openmp very often, so I'll need to investigate...

[mreineck]

OK, I have an idea why the Fortran code cannot be called concurrently: in the function you linked there are arrays with the "save" attribute, which means that they are equivalent to static variables in C/C++, i.e. global state. Calling into such a function concurrently will almost certainly lead to undesired effects.

Not sure why you also encounter this with the ducc back-end as well though ... it shouldn't have any global state.

Semi-related question: how do you manage to execute concurrent calls from Python? Doesn't the GIL prevent this?

[mreineck]

Sorry, I thake that back: all variables marked save are also threadprivate, so this could in principle work.

On the other hand, threadprivate is an OpenMP-specific keyword, so I have no idea what happens if this is called from multiple non-OpenMP threads...

[aafkegros]

I had to look it up to figure out how the software im using is doing this: I'm never releasing the GIL, I call SHExpandGLQ in multiple consecutive threads (that exist for other non-GIL-locked reasons outside of my code) and it seems like the entering into SHExpandGLQ unlocks the GIL somewhere, causing another of my threads to progress until it collides inside of pyshtools.

[aafkegros]

The GIL unlocking was proposed in #304 and is done here it seems

[mreineck]

Sure! I was just wondering how you can call the SHTOOLS Python interface from two concurrently running Python threads. I most likely have an oversimplified picture in my mind of what you are actually doing.

[mreineck]

To get a clearer picture: which function from the SHTOOLS Python interface are you calling exactly?

[aafkegros]

the only function i have isssues with is SHExpandGLQ, if i lock this function to wait for thread execution, my code works fine

[aafkegros]

And i have to say i dont fully understand your question on having multiple concurrent python threads, this is something that does not require the release of the GIL (such as the multiprocessing.ThreadPool implementation)

[aafkegros]

I'll try to make a minimal working example soon 😄

[mreineck]

Thanks, that will certainly be very helpful!

[aafkegros]

from concurrent import futures
import pyshtools as pysh
import numpy as np

your_patches = np.random.randint(1,1e6, size=(100, 251,503))

zero, w = pysh.expand.SHGLQ(251)
with futures.ThreadPoolExecutor(max_workers=10) as executor:
    jobs = [executor.submit(pysh.expand.SHExpandGLQ, patch, w, zero) for patch in your_patches]
    [fut.result() for fut in futures.as_completed(jobs)]

This breaks with a segmentation fault for me, but works with max_workers=1 . The 251 lmax is just because i use this as a default grid size (~80*π), but it also breaks with other lmax that i tested