用OpenMP方程求解Gauss

Question

，我在Gauss上做了一个代码，用OpenMP来求解线性方程组。但我没有得到正确的输出。如果n设置为1，则输出是正确的。不过，对我来说，算法似乎是正确的。有人能帮忙吗？

代码：

    int seidel_p(double *A, double *b, double *x, int n)
{
    double *dx;
    dx = (double*) calloc(n,sizeof(double));
    int epsilon = 1.0e-4;
    int i,j,k,id;
    double dxi;
    double sum;
    for(int i =0; i<n ; i++)
    {
        x[i] = 0;
    }

    int maxit = 2*n*n;

    for(k=0; k<maxit; k++)
   {
      for(i=0; i<n; i++)
      {
         dx[i] = b[i];
         #pragma omp parallel shared(A,x)
         {
            for(j=0; j<n; j++)
            {
                if(i!=j)
                {
                    dxi += A[i*n+j]*x[j];
                }
            }

            #pragma omp critical
               dx[i] -= dxi;
         }
         dx[i] /= A[i*n +i];
         x[i] = dx[i];
         sum += ( (dx[i] >= 0.0) ? dx[i] : -dx[i]);
         if(sum<= epsilon) break;

      }
   }

   for( int i = 0 ; i<n ; i++ )
   {
        printf("\t\t%f\n", dx[i]);
   }

    free(dx);
}

Answer 1

除非你有很好的理由自己去做这件事，否则不要这样做。有很多库会这样做(还有更多的)。艾根是一个标头唯一的库，检查它。

http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html

您还可以使用BLAS / LAPACK，其中有许多实现。

https://github.com/Reference-LAPACK/lapack (这个实现不是特别优化的)。

如果您坚持要滚动自己的实现，那么首先应该检查的是是否只使用一个线程。您可以通过设置环境变量OMP_NUM_THREADS等于1或通过回调omp_set_num_threads(1)来尝试这一点。如果这样做有效，那么您就知道您的算法是正确的，只是您对OpenMP的使用失败了。