Re: parallel C code -- help needed

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Tue, May 31, 2022 at 09:44:49AM +0100, John Cremona wrote:
> I have a C++ program which uses libpari in just one function, and I
> want to parallelise the program using omp.    I have done this
> successfully with other programs not using libpari, so I think I have
> the omp basics right (while being a novice in parallel programming).
> My libpari was configured with --mt=pthread.

This is an important point, see below.

> A simplified version of the function looks like this before
> parallelising, and this works fine (the array *ai passed has length
> 5):
> 
> long pari_sturm(long *ai, int pos_only=0, int neg_only=0)
> // Return the number of real roots (default), number of positive real
> // roots (if pos_only==1) or number of negative real roots (if
> // neg_only==1)
> {
>   long res;
>   pari_sp av = avma;
>   GEN g0 = stoi(ai[0]);
>   GEN g1 = stoi(ai[1]);
>   GEN g2 = stoi(ai[2]);
>   GEN g3 = stoi(ai[3]);
>   GEN g4 = stoi(ai[4]);
>   GEN f = mkpoln(5,g0,g1,g2,g3,g4);
>   f = gdiv(f,ggcd(f,derivpol(f)));
> 
>   if (pos_only)
>     res = sturmpart(f,gen_0,NULL); // #roots >=0
>   else
>     {
>       if (neg_only)
>         res = sturmpart(f,NULL,gen_0); // #roots <=0
>       else
>         res = sturm(f);                  // #roots
>     }
>   gerepileupto(av, stoi(0));
>   return res;
> }
> 
> Now I run in 4 threads, set by inserting
> #pragma omp parallel num_threads(4)
> in a higher level function; each thread will call pari_sturm() many
> times.  Since pari_sturm()  uses at least one global variable (avma) I
> enclose the code in an omp critical block so it looks like
> 
> long pari_sturm(long *ai, int pos_only=0, int neg_only=0)
> {
>   long res;
> #pragma omp critical(pari)
> {
>   pari_sp av = avma;
>   GEN g0 = stoi(ai[0]);
>   <same as above>
>   gerepileupto(av, stoi(0));
> }
>  return res;
> }
> 
> but when I run it I get a segfault on the line "GEN g0 = stoi(ai[0]);"

> Can this be made to work?

Using a critical section to marshall access to the PARI stack could work
but only with the non-TLS version of PARI!

When using the TLS version of PARI (--enable-tls or --mt=pthread) each
thread automatically get its own private (TLS) avma and stack, which
need to be initialized when the thread starts.

This is done by using pari_thread_alloc in the main thread and 
pari_thread_start/pari_thread_close in each secondary threads.

The recommended way to use openmp with PARI is explained in the file
examples/openmp.c and the appendix D.

Cheers,
Bill.