I tried to call the Monte Carlo integration subroutine of GSL library to do some numerical calculation. Because my for loop is rather simple, meaning the results of different runs are independent, I expect it should be very straightforward to parallelize using OpenMP. However when I compiled it, it always said "Internal Compiler Error: Segmentation Fault", and generated nothing. Here is my code:
#include <stdlib.h>
#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <math.h>
double
Reg1 (double *x, double t, size_t dim, void *params)
{
return sin(x[1])*cos(t*t)*x[0]*x[0]*cos(x[0])*cos(3*x[0]);
}
void
display_results (char *title, double result, double error)
{
printf ("%s ==================\n", title);
printf ("result = % .10f\n", result);
printf ("sigma = % .10f\n\n", error);
}
void
VEGAS_integration_routine (double (*function)(double *, size_t, void *),
double xl[], double xu[], size_t calls, gsl_rng * r)
{
double res, err;
gsl_monte_function Function = { function, 2, 0 };
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (2);
gsl_monte_vegas_integrate (&Function, xl, xu, 2, 20000, r, s, &res, &err);
display_results ("vegas warm-up", res, err);
printf ("converging...\n");
do
{
gsl_monte_vegas_integrate (&Function, xl, xu, 2, calls, r, s, &res, &err);
printf ("result = % .10f; sigma = % .10f; "
"chisq/dof = %.2f\n", res, err, gsl_monte_vegas_chisq (s));
}
while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.05);
display_results ("vegas final", res, err);
gsl_monte_vegas_free (s);
}
int
main (void)
{
int seg = 200;
double xl[2] = { 195.0, -1000.0 };
double xu[2] = { 205.0, 1000.0 };
const gsl_rng_type *T;
gsl_rng *r;
size_t calls = 1000000;
gsl_rng_env_setup ();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
/* Calculate G1 */
int i;
double j=0;
#pragma omp parallel for shared(xl,xu,calls,r,seg) private(i,j)
for(i=0; i<=seg; i=i+1)
{
j=(double)i * 0.05;
printf("Now it's t = %.2f \n", j);
printf("Compute Re(G1)...\n");
double g(double * x, size_t dim, void *params)
{
return Reg1(x, j, dim, ¶ms);
}
VEGAS_integration_routine (&g, xl, xu, calls, r);
}
gsl_rng_free (r);
return 0;
}
This code is basically modified from the sample code on GSL webpage. Without using OpenMP, it works perfectly. But when I used gcc to compile with the following commands (with -fopenmp added) , which works on our server,
gcc -c -Wall -ansi -I/usr/include -L/usr/lib/gcc/x86_64-redhat-linux/4.4.4/ -lgomp -fopenmp test2.c -o test2.o
gcc -L/usr/lib64/ test2.o -L/usr/lib/gcc/x86_64-redhat-linux/4.4.4/ -lgomp -fopenmp -lgsl -lgslcblas -lm -o test2.out
I got the following error message:
test2.c: In function 'main':
test2.c:85: internal compiler error: Segmentation fault
Please submit a full bug report,
with preprocessed source if appropriate.
See <http://bugzilla.redhat.com/bugzilla> for instructions.
Preprocessed source stored into /tmp/ccAGFe3v.out file, please attach this to your bugreport.
make: *** [test2.o] Error 1
Since I couldn't compile it and the error message shown above was so limited, I'd really like to know what's wrong, so I decompose the subroutine VEGAS_integration_routine I called and then to run it line by line. What I found is that the compilation stopped at its second line
gsl_monte_function Function = { function, 2, 0 };
which made me so confused. Can't I declare a GSL function in a loop when using OpenMP to flatten the loop? Is there an internal conflict between GSL and OpenMP? I did search on Stack Overflow as well as Google, but it seems that no related post exists (so odd!).I would really appreciate if anyone could either explain what's going on here, or point out an alternative way to do parallel computing. I'm sure the way I wrote my for loop wasn't the best and the neatest...
There is a regression bug in the way the implementation of lexical closures in GCC interacts with the OpenMP code transformation engine. The bug appears to be fixed in later versions of the GCC 4.7 branch. If you cannot replace GCC 4.4.4 with GCC 4.7.1 or later, then you should rewrite your code to not use nested functions as in the GSL example that you refer in your post. Besides, nested functions are non-standard GCC C extension and you should refrain from using them wherever possible.
The Monte Carlo integration interface in GSL supports passing additional arguments to the integrand via the void *params argument as explained here. You would have to change Reg1 as follows:
double
Reg1 (double *x, size_t dim, void *params)
{
double t = *(double *)params;
return sin(x[1])*cos(t*t)*x[0]*x[0]*cos(x[0])*cos(3*x[0]);
}
and update the second line of VEGAS_integration_routine to read:
gsl_monte_function Function = { function, 2, &t };
You will also have to add t to the list of dummy arguments of VEGAS_integration_routine and pass its value from main. The relevant parts therefore become:
void
VEGAS_integration_routine (double (*function)(double *, size_t, void *),
double t,
double xl[], double xu[], size_t calls, gsl_rng * r)
{
double res, err;
gsl_monte_function Function = { function, 2, &t };
...
}
#pragma omp parallel for
for(i=0; i <= seg; i++)
{
double t = (double)i * 0.05;
...
VEGAS_integration_routine (Reg1, t, xl, xu, calls, r);
}