Re: parallel C code -- help needed

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

On Tue, 31 May 2022 at 14:14, John Cremona <john.cremona@gmail.com> wrote:
>
> 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.

Fantastic, it now works perfectly.  I just added 2 lines outside the
parallel block to define the pari_thread object pth, and call
pari_thread_alloc on each.  Inside the parallel block I call
pari_thread_start() at the beginning and pari_thread_close() at the
end.  And I removed the #pragma omp critical in the function which
does the main pari work.

Thanks, Bill!

John

>
> John
>
> >
> > Cheers,
> > Bill.