openmp - продукт с параллельными векторными матрицами

Question

Я вычисляю векторную матрицу (матрица разрежена и хранится в CSR) продукта (не совсем продукт, небольшое отклонение для вычисления кратчайшего расстояния) с использованием метода внешнего продукта. Я новичок в параллельном программировании и, по сути, пытаюсь понять разницу между использованием параллели для секции с критической секцией для обновления VS с использованием задач и выполнением редукции. Какой подход лучше и почему?

Примечание. Этот вызов функции заключен в параллельный omp и одиночный omp.

Используя параллель для подхода, я делаю это:

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq) {
    initialize_T(tReq);
    int index;
    #pragma omp parallel for schedule(static, BLOCK_SIZE)
    for(int i=0;i<N;i++) {
        int num_edges = A->row_ptr[i+1] - A->row_ptr[i];
        index = 0;
        if(num_edges) {
            if(T[i] != INFINITY && tB[i] != INFINITY) {
                for(int j=0;j<num_edges;j++) {
                    index = A->col_ind[A->row_ptr[i] + j];
                    #pragma omp critical 
                    tReq[index] = min(tReq[index], T[i]+A->val[A->row_ptr[i]+j]);      
                }
            }
        }
    }
    return tReq;
}

Используя основанный на задачах подход с сокращением, это по сути моя идея:

int size = N/BLOCK_SIZE + 1;
double C[size][N];

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq, int size, double C[][N]) {

    int index;

    for(int i=0;i<size;i++) {
        for(int j=0;j<N;j++) {
            C[i][j] = INFINITY;
            tReq[j] = INFINITY;
        }
    }

    int k;

    for(k=0;k<size-1; k++) {
        #pragma omp task firstprivate(k) depend(inout: C[k]) 
        {
            int index;
            int delim;   
            delim = (k==size-1) ? N : BLOCK_SIZE;
            printf("kk is %d\n", k*BLOCK_SIZE);
            // printf("k is %d Delim is %d for thread %d\n", k, delim, omp_get_thread_num());
            for(int i=0;i<delim; i++) {
                int num_edges = A->row_ptr[k*BLOCK_SIZE + i+1] - A->row_ptr[k*BLOCK_SIZE + i];
                index = 0;
                if(num_edges) {
                    if(T[k*BLOCK_SIZE + i] != INFINITY && tB[k*BLOCK_SIZE + i] != INFINITY) {           
                        for(int j=0;j<num_edges;j++) {
                            index = A->col_ind[A->row_ptr[k*BLOCK_SIZE + i] + j];
                            {
                            C[k][index] = min(C[k][index], T[k*BLOCK_SIZE + i]+A->val[A->row_ptr[k*BLOCK_SIZE + i]+j]);                 
                            }
                        }
                    }       
                }   
            }
        }
    }    

    #pragma omp taskwait

    for(int i=0; i<N; i++) {
        {
            double minimum = INFINITY;
            for(int j=0;j<size;j++) {
                if(C[j][i] < minimum) {
                    minimum = C[j][i];
                }
            }
            tReq[i] = minimum;
        }
    }

    return tReq;
}

По существу, есть ли недостатки в использовании параллельного интерфейса по сравнению с подходом, основанным на задачах?

Answer 1

Вы правы в том, что у вас есть два основных варианта: защитить обновление данных или использовать специфичные для потока копии. Тем не менее, вы можете сделать намного лучше для каждого варианта:

При работе с защищенными обновлениями вы должны защищать только то, что и когда это абсолютно необходимо. Вы можете использовать начальную атомную проверку для предотвращения критических областей в большинстве случаев, аналогично двойной проверке шаблона блокировки.

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq) {
    initialize_T(tReq);
    #pragma omp parallel for schedule(static, BLOCK_SIZE)
    for(int i=0;i<N;i++) {
        int num_edges = A->row_ptr[i+1] - A->row_ptr[i];
        if (num_edges) {
            if(T[i] != INFINITY && tB[i] != INFINITY) {
                for(int j=0;j<num_edges;j++) {
                    // !WARNING! You MUST declare index within the parallel region
                    // or explicitly declare it private to avoid data races!
                    int index = A->col_ind[A->row_ptr[i] + j];
                    double tmp = T[i] + A->val[A->row_ptr[i]+j];
                    double old;
                    #pragma omp atomic
                    old = tReq[index];
                    if (tmp < old) {
                        #pragma omp critical
                        {
                            tmp = min(tReq[index], tmp);
                            // Another atomic ensures that the earlier read
                            // outside of critical works correctly
                            #pragma omp atomic
                            tReq[index] = tmp;
                        }
                    }
                }
            }
        }
    }
    return tReq;
}

Примечание. К сожалению, ни один OpenMP / C не поддерживает прямой атомарный минимум.

Альтернативой является сокращение, которое поддерживается даже самим стандартом. Таким образом, нет необходимости заново изобретать разделение работы и т. Д. Вы можете просто сделать следующее:

double *matrixVectorHadamard(CSR *A, double *T, double *tB, double *tReq) {
    initialize_T(tReq);
    #pragma omp parallel for schedule(static, BLOCK_SIZE) reduction(min:tReq)
    for(int i=0;i<N;i++) {
        int num_edges = A->row_ptr[i+1] - A->row_ptr[i];
        if (num_edges) {
            if(T[i] != INFINITY && tB[i] != INFINITY) {
                for(int j=0;j<num_edges;j++) {
                    // !WARNING! You MUST declare index within the parallel region
                    // or explicitly declare it private to avoid data races!
                    int index = A->col_ind[A->row_ptr[i] + j];
                    tReq[index] = min(tReq[index], T[i]+A->val[A->row_ptr[i]+j]);
                }
            }
        }
    }
    return tReq;
}

OpenMP волшебным образом создаст локальные для потока копии tReq и объединит (уменьшит) их в конце.

Какая версия лучше для вас, зависит от размера целевого массива и скорости записи. Если вы часто пишете, сокращение будет полезным, потому что оно не замедляется из-за плохого кэширования critical / atomic. Если у вас огромный целевой массив или не так много итераций обновления, первое решение становится более интересным из-за относительных накладных расходов на создание и уменьшение массива.