Как я могу использовать openmp и AVX2 одновременно с идеальным ответом? - PullRequest
Новая
       18

Как я могу использовать openmp и AVX2 одновременно с идеальным ответом?

0 голосов
/ 01 июля 2018

Я написал программу продукта Matrix-Vector с использованием OpenMP и AVX2.

Однако я получил неправильный ответ из-за OpenMP. Истинный ответ - все значение массива c станет 100.

Мой ответ был смесью 98, 99 и 100.

Фактический код указан ниже.

Я скомпилировал Clang с -fopenmp, -mavx, -mfma. 

#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include "omp.h"
#include "x86intrin.h"

void mv(double *a,double *b,double *c, int m, int n, int l)
{
    int k;
#pragma omp parallel
    {
        __m256d va,vb,vc;
        int i;
#pragma omp for private(i, va, vb, vc) schedule(static)
        for (k = 0; k < l; k++) {
            vb = _mm256_broadcast_sd(&b[k]);
            for (i = 0; i < m; i+=4) {
                va = _mm256_loadu_pd(&a[m*k+i]);
                vc = _mm256_loadu_pd(&c[i]);

                vc = _mm256_fmadd_pd(vc, va, vb);

                _mm256_storeu_pd( &c[i], vc );
            }
        }
    }
}
int main(int argc, char* argv[]) {

    // set variables
    int m;
    double* a;
    double* b;
    double* c;
    int i;

    m=100;
    // main program

    // set vector or matrix
    a=(double *)malloc(sizeof(double) * m*m);
    b=(double *)malloc(sizeof(double) * m*1);
    c=(double *)malloc(sizeof(double) * m*1);
    //preset
    for (i=0;i<m;i++) {
        a[i]=1;
        b[i]=1;
        c[i]=0.0;
    }
    for (i=m;i<m*m;i++) {
        a[i]=1;
    }

    mv(a, b, c, m, 1, m);

    for (i=0;i<m;i++) {
        printf("%e\n", c[i]);
    }
    free(a);
    free(b);
    free(c);
    return 0;
}

Я знаю, что критический раздел поможет. Однако критическая секция была медленной.

Итак, как мне решить проблему?

Ответы [ 2 ]

0 голосов
/ 02 июля 2018

Основная операция, которую вы хотите это 

c[i] = a[i,k]*b[k]

Если вы используете хранилище с мажорными строками , оно становится 

c[i] = a[i*l + k]*b[k]

Если вы используете хранилище ордеров с основным столбцом, оно становится 

c[i] = a[k*m + i]*b[k]

Для порядка строк вы можете распараллелить так: 

#pragma omp parallel for
for(int i=0; i<m; i++) {
  for(int k=0; k<l; k++) {
    c[i] += a[i*l+k]*b[k];
  }
}

Для порядка столбцов вы можете распараллелить так: 

#pragma omp parallel
for(int k=0; k<l; k++) {
  #pragma omp for
  for(int i=0; i<m; i++) {
    c[i] += a[k*m+i]*b[k];
  }
}

Матрично-векторные операции - это операции уровня 2, которые являются операциями с ограниченной пропускной способностью памяти. Операции уровня 1 и уровня 2 не масштабируются, например, с количеством ядер. Это только операции уровня 3 (например, умножение плотных матриц), которые масштабируются https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms#Level_3.

0 голосов
/ 01 июля 2018

Проблема не в ваших встроенных функциях AVX, давайте посмотрим на код без встроенных функций на минуту: 

void mv(double *a,double *b,double *c, int m, int n, int l)
{
    #pragma omp parallel for schedule(static)
    for (int k = 0; k < l; k++) {
        double xb = b[k];
        for (int i = 0; i < m; i++) {
            double xa = a[m*k+i];
            double xc = c[i];
            xc = xc + xa * xb;
            c[i] = xc;
        }
    }
}

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

Состояние гонки в вашем коде - c[i] - который пытаются обновить несколько потоков. Теперь, даже если бы вы могли защитить это, скажем, с помощью атомарного обновления, производительность была бы ужасной: не только из-за защиты, но и потому, что данные c[i] должны постоянно перемещаться между кэшами разных ядер.

Одна вещь, которую вы можете сделать с этим, это использовать сокращение массива на c. Это делает личную копию c для каждого потока, и в конце они объединяются: 

void mv(double *a,double *b,double *c, int m, int n, int l)
{
    #pragma omp parallel for schedule(static) reduction(+:c[:m])
    for (int k = 0; k < l; k++) {
        for (int i = 0; i < m; i++) {
            c[i] += a[m*k+i] * b[k];
        }
    }
}

Это должно быть достаточно эффективным, если два m -вектора вписываются в ваш кеш, но вы все равно можете получить много накладных расходов из-за накладных расходов на управление потоками. В конечном итоге вы будете ограничены пропускной способностью памяти, потому что при умножении векторной матрицы у вас есть только одно вычисление на элемент, считанный из a.

В любом случае, вы, конечно, можете поменять местами циклы i и k и сохранить сокращение, но тогда ваш шаблон доступа к памяти на a будет неэффективным (пошаговым) - так что вам следует block цикл, чтобы избежать этого.

Теперь, если вы посмотрите на вывод любого современного компилятора , он сгенерирует SIMD-код сам по себе. Конечно, вы можете применить свои собственные SIMD-функции, если хотите. Но убедитесь, что вы правильно обрабатываете крайние случаи, если m не делится на 4 (вы не делали это в исходной версии).

В конце концов, если вы действительно хотите повысить производительность - используйте функции из библиотеки BLAS (например, MKL). Если вы хотите поэкспериментировать с оптимизацией, у вас будет множество возможностей углубиться в детали.

...