я есть реализация матричного решателя на основе BiCCG (Conjugate Gradient), который также учитывает периодичность. Случается, что реализация требует значительных вычислительных ресурсов, и цикл не векторизован автоматически из-за проблемы с зависимостями. Я немного исследовал, и кажется, что красно-черный алгоритм Гаусса Зейделя более эффективно распараллеливается, чем ванильная версия (которая также имеет аналогичную проблему зависимости).

Может ли этот цикл / алгоритм быть изменен так, чтобы он мог эффективно векторизоваться?

 // FORWARD
        #pragma omp for schedule(static, nx/NTt)
        for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++)
        {


            dummy = res_sparse_s[i][j][k];

                                           dummy -= COEFF[i][j][k][7] * RLL[i-1][j][k];
            if (PeriodicBoundaryX && i==nx)dummy -= COEFF[i][j][k][8] * RLL[1  ][j][k];


                                            dummy -= COEFF[i][j][k][2] * RLL[i][j-1][k];
            if (PeriodicBoundaryY && j==ny) dummy -= COEFF[i][j][k][3] * RLL[i][1  ][k];


                                            dummy -= COEFF[i][j][k][4] * RLL[i][j][k-1];
            if (PeriodicBoundaryZ && k==nz) dummy -= COEFF[i][j][k][5] * RLL[i][j][1  ];


            RLL[i][j][k] = dummy / h_sparse_s[i][j][k];
        }

Постскриптум RLL для любой итерации i, j, k включает в себя обновленный «RLL» в i-1, j-1 и k-1 через переменную dummy. Также теперь цикл разрушается только в направлении х с помощью директивыschedule(static, nx/NTt) где NTt - это просто макрос для доступного количества потоков. Можно ли разбить его во всех направлениях, используя директиву?collapse?

------- ОСНОВНОЕ РЕДАКТИРОВАНИЕ -------------------------- следующий ответ Аджая - это минимальный рабочий пример

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<omp.h>

typedef double lr;

#define nx 4
#define ny 4
#define nz 4

void
print3dmatrix(double a[nx+2][ny+2][nz+2])
{
    for(int i=1; i<= nx; i++) {
        for(int j=1; j<= ny; j++) {
            for(int k=1; k<= nz; k++) {
                printf("%f ", a[i][j][k]);
            }
            printf("\n");
        }
        printf("\n");
    }
}

int 
main()
{

    double a[nx+2][ny+2][nz+2];
    double b[nx+2][ny+2][nz+2];

    srand(3461833726);


    // matrix filling 
    // b is just a copy of a
    for(int i=0; i< nx+2; i++) for(int j=0; j< ny+2; j++) for(int k=0; k< nz+2; k++)
    {
        a[i][j][k] = rand() % 5;
        b[i][j][k] = a[i][j][k];
    }

    // loop 1
    //#pragma omp parallel for num_threads(1)
    for(int i=1; i<= nx; i++) for(int j=1; j<= ny; j++) for(int k=1; k<= nz; k++)
    {
        a[i][j][k] = -1*a[i-1][j][k] - 1*a[i][j-1][k] -1 * a[i][j][k-1] + 4 * a[i][j][k];
    }

    print3dmatrix(a);
    printf("******************************\n");

    // loop 2
    //#pragma omp parallel for num_threads(1)
    for(int i=1; i<= nx; i++) 
        for(int j=1; j<= ny; j++)
            // #pragma omp simd
            for(int m=j+1; m<= j+nz; m++)
            {
                b[i][j][m-j] = -1*b[i-1][j][m-j] - 1*b[i][j-1][m-j] -1 * b[i][j][m-j-1] + 4 * b[i][j][m-j];
            }

    print3dmatrix(b);
    printf("=========================\n");

    return 0;
}

Ключевые наблюдения-

Матрица а заполнена случайными числами от 0 до 5, и цикл 1 является не преобразованным исходным циклом, тогда как цикл 2 является преобразованным цикломПреобразованный цикл был перекошен для удаления зависимости.После операции матрицы a и b одинаковы, если запустить без распараллеливания openmpесли развернут открытый mp, ответы меняются (возможно из-за расы) [цикл не может быть паралеллизуемым независимо от того, где размещена прагма]если#pragma omp simd используется для принудительной векторизации самого внутреннего цикла, в котором происходит сбой.