3

我尝试调用 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 网页上的示例代码修改而来的。在不使用 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 上进行了搜索，但似乎没有相关的帖子存在（太奇怪了！）。如果有人能解释这里发生了什么，或者指出另一种进行并行计算的方法，我将不胜感激。我确定我编写 for 循环的方式不是最好和最整洁的......

4

1 回答 1

4

GCC 中的词法闭包实现与 OpenMP 代码转换引擎交互的方式存在回归错误。该错误似乎已在 GCC 4.7 分支的更高版本中得到修复。如果您无法将 GCC 4.4.4 替换为 GCC 4.7.1 或更高版本，那么您应该重写代码以不使用您在帖子中引用的 GSL 示例中的嵌套函数。此外，嵌套函数是非标准的 GCC C 扩展，您应该尽可能避免使用它们。

GSL 中的 Monte Carlo 积分接口支持通过参数将附加参数传递给被积函数，如此处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);
}
于 2013-05-20T12:00:12.477 回答