Внутренняя ошибка компилятора с вложенными функциями в параллельных регионах OpenMP

Question

Я попытался вызвать подпрограмму интеграции Монте-Карло библиотеки GSL, чтобы провести какой-то численный расчет. Поскольку цикл for довольно прост, то есть результаты разных запусков независимы, я ожидаю, что распараллеливание с использованием OpenMP должно быть очень простым. Однако, когда я компилировал его, он всегда говорил "Внутренняя ошибка компилятора: ошибка сегментации" и ничего не генерировал. Вот мой код:

#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, &params);
        }
        VEGAS_integration_routine (&g, xl, xu, calls, r);
   }
   gsl_rng_free (r);

   return 0;
 }

Этот код в основном изменен из примера кода на веб-странице GSL. Без использования OpenMP все работает отлично. Но когда я использовал gcc для компиляции с помощью следующих команд (с -fopenmp добавил), который работает на нашем сервере,

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

Я получил следующее сообщение об ошибке:

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

Так как я не смог скомпилировать его, а сообщение об ошибке, показанное выше, было настолько ограниченным, я бы очень хотел знать, что не так, поэтому я разлагаю подпрограмму VEGAS_integration_routine Я позвонил, а затем запустить его построчно. Я обнаружил, что компиляция остановилась на второй строке

gsl_monte_function Function = { function, 2, 0 };

что меня так смутило Не могу ли я объявить функцию GSL в цикле при использовании OpenMP для выравнивания цикла? Есть ли внутренний конфликт между GSL и OpenMP? Я выполнил поиск по переполнению стека так же, как и по Google, но похоже, что не существует связанных постов (так странно!). Я был бы очень признателен, если бы кто-нибудь мог либо объяснить, что здесь происходит, либо указать альтернативный способ параллельных вычислений. Я уверен, что способ, которым я написал свой цикл for, был не лучшим и не самым лучшим...

Answer 1

Существует ошибка регрессии в том, как реализация лексических замыканий в GCC взаимодействует с механизмом преобразования кода OpenMP. Эта ошибка исправлена ​​в более поздних версиях GCC 4.7. Если вы не можете заменить GCC 4.4.4 на GCC 4.7.1 или новее, вам следует переписать свой код, чтобы не использовать вложенные функции, как в примере GSL, который вы ссылаетесь в своем посте. Кроме того, вложенные функции являются нестандартным расширением GCC C, и вы должны воздерживаться от их использования везде, где это возможно.

Интерфейс интеграции Монте-Карло в GSL поддерживает передачу дополнительных аргументов в интегратор через void *params аргумент, как объяснено здесь. Вы должны были бы изменить Reg1 следующее:

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]);
}

и обновить вторую строку VEGAS_integration_routine читать:

gsl_monte_function Function = { function, 2, &t };

Вы также должны будете добавить t к списку фиктивных аргументов VEGAS_integration_routine и передать его значение из main, Поэтому соответствующие части становятся:

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);
}