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...