OpenMP並行區域中嵌套函數的內部編譯器錯誤

Question

我試圖調用GSL庫的Monte Carlo集成子例程來進行一些數值計算。因爲我的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 webpage。不使用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 }; 

這讓我很困惑。在使用OpenMP壓扁循環時，我無法在循環中聲明GSL函數嗎？ GSL和OpenMP之間是否存在內部衝突？我曾在Stack Overflow和Google上搜索過，但似乎沒有相關的文章存在（太奇怪了！）。如果有人能夠解釋這裏發生了什麼，或者指出另一種做並行計算的方法，我將非常感激。我敢肯定，我寫循環的方式並不是最好的，最好的...

Answer 1

在GCC的詞法關閉實現與OpenMP代碼轉換引擎交互的方式中有a regression bug。該錯誤似乎在GCC 4.7分支的更高版本中得到修復。如果你不能用GCC 4.7.1或更高版本替換GCC 4.4.4，那麼你應該重寫你的代碼，使其不像你在文章中提到的GSL示例那樣使用嵌套函數。此外，嵌套函數是非標準的GCC C擴展，你應該儘可能避免使用它們。

GSL中的Monte Carlo集成接口支持通過void *params參數將其他參數傳遞給被積函數，如here所解釋的。你將不得不改變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); 
}