Блог

Я делаю код на Gauss Seidel с использованием 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("tt%fn", dx[i]);
   }

    free(dx);
}
  

Если у вас нет ОЧЕНЬ веской причины делать это самостоятельно, не делайте этого. Существует множество библиотек, которые делают это (и многое другое). Eigen — это библиотека только для заголовков, проверьте это.

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, которое нарушено.