Re: parallel C code -- help needed

[John Cremona <john.cremona@gmail.com>]

On Tue, 31 May 2022 at 13:34, Bill Allombert
<Bill.Allombert@math.u-bordeaux.fr> wrote:
>
> 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.

Thanks a lot Bill, I will try your suggestions.

John

>
> Cheers,
> Bill.