John Cremona on Tue, 31 May 2022 10:45:19 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

parallel C code -- help needed


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.

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?

John