| 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