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