Улучшение производительности кэша и распараллеливание openMP ⇐ C++

[Anonymous]

Я работаю с большими трехмерными векторными полями, которые представляют некоторые физические величины в физическом пространстве или пространстве Фурье. Векторные поля преобразуются Фурье с использованием библиотеки FFTW3. Я пытаюсь улучшить общую производительность кода.
Компоненты векторного поля хранятся в
typedef double fftw_complex[2];

fftw_complex* cVec[3]
for (int dim = 0; dim < 3; dim++)
{
cVec[dim] = fftw_alloc_complex(dataSize);
}

Каждый компонент векторного поля хранится последовательно в памяти, которая «правильно» выровнена с помощью процедур выделения FFTW. На самом деле, я не проверял, есть ли падение производительности при нераспределении памяти с помощью fftw_malloc, собираюсь это сделать. Чтобы инициализировать это векторное поле, я бы позвонил
void fillUnAligned(fftw_complex* const c_inout[3], double value)
{
#ifdef USE_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < cDataSize; i++)
{
for (int dim = 0; dim < 3; dim++)
{
#ifdef USE_OMP
#pragma omp simd
#endif
for (int ccomp = 0; ccomp < 2; ccomp++)
{
c_inout[dim][ccomp] = value;
}
}
}
}

однако время выполнения остается неизменным независимо от количества потоков (до 24). Код скомпилирован с помощью
g++-11 -fopenmp -O3 -march=native

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

[*]выравнивании памяти и
[*]шаблонах доступа к памяти.

В большинстве случаев все три компонента вектора должны обрабатываться за одну итерацию над векторным полем. Память выровнена по 16 байт, а размер строки кэша составляет 64 байт. Если я не ошибаюсь, 4 элемента типа fftw_complex заполняют всю строку кэша и, таким образом, элементы векторного поля пересекают границы строки кэша. Из-за способа определения векторного поля (в основном для удобства применения БПФ) последовательные векторы не хранятся в памяти последовательно, в том смысле, что три компонента одного и того же вектора могут быть широко разделены. Чтобы решить эти проблемы, я определил простую структуру, представляющую трехмерный вектор в пространстве Фурье
#define CACHE_LINE 64

struct alignas(CACHE_LINE) ComplexVector
{
ComplexVector()
{
for (int dim = 0; dim < 3; dim++)
{
for (int ccomp = 0; ccomp < 2; ccomp++)
{
components[dim][ccomp] = 0.;
}
}
}

fftw_complex const& operator[](int index) const
{
return components[index];
}

fftw_complex& operator[](int index)
{
return components[index];
}

private:
fftw_complex components[3];
char padding[16];
};

Весь вектор в пространстве Фурье имеет размер 48 байт, поэтому, чтобы предотвратить пересечение векторами границ строки кэша, я добавил в конце отступ размером 16 байт. У меня нет опыта работы в CS, так что это может быть не самое лучшее решение. Является ли введение такого принудительного согласования хорошей/плохой/обычной практикой? Теперь векторное поле можно переопределить как
std::vector cVec(dataSize);

где каждый ComplexVector заполняет всю строку кэша, поэтому (я считаю) во время каждой итерации по векторному полю все компоненты могут обрабатываться более эффективно. Чтобы инициализировать векторное поле, я бы сейчас позвонил
void fillAligned(std::vector& c_inout, double value)
{
#ifdef USE_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < cDataSize; i++)
{
#ifdef USE_OMP
#pragma omp simd collapse(2)
#endif
for (int dim = 0; dim < 3; dim++)
{
for (int ccomp = 0; ccomp < 2; ccomp++)
{
c_inout[dim][ccomp] = value.;
}
}
}
}

К сожалению, второе решение (я считал его более эффективным) примерно на 40 % медленнее первого. OpenMP ничего не меняет. В чем причина этого? Как я могу контролировать производительность этого кода? Ниже приведен минимальный рабочий пример
#include
#include
#include
#include

using namespace std::chrono;

#define USE_OMP

#ifdef USE_OMP
#include
#define NUM_THREADS 24
#endif

#define CACHE_LINE 64
#define NUM_NODES 256

constexpr int cDataSize { NUM_NODES * NUM_NODES * (NUM_NODES / 2 + 1) };

typedef double fftw_complex[2];

struct alignas(CACHE_LINE) ComplexVector
{
ComplexVector()
{
for (int dim = 0; dim < 3; dim++)
{
for (int ccomp = 0; ccomp < 2; ccomp++)
{
components[dim][ccomp] = 0.;
}
}
}

fftw_complex const& operator[](int index) const
{
return components[index];
}

fftw_complex& operator[](int index)
{
return components[index];
}

private:
fftw_complex components[3];
char padding[16];
};

void fillUnAligned(fftw_complex* const c_inout[3], double value)
{
#ifdef USE_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < cDataSize; i++)
{
for (int dim = 0; dim < 3; dim++)
{
#ifdef USE_OMP
#pragma omp simd
#endif
for (int ccomp = 0; ccomp < 2; ccomp++)
{
c_inout[dim][ccomp] = value;
}
}
}
}

void fillAligned(std::vector& c_inout)
{
#ifdef USE_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < cDataSize; i++)
{
#ifdef USE_OMP
#pragma omp simd collapse(2)
#endif
for (int dim = 0; dim < 3; dim++)
{
for (int ccomp = 0; ccomp < 2; ccomp++)
{
c_inout[dim][ccomp] = 0.;
}
}
}
}

void copyAligned(std::vector& c_in,
std::vector& c_out)
{
#ifdef USE_OMP
#pragma omp parallel for
#endif
for (int i = 0; i < cDataSize; i++)
{
#ifdef USE_OMP
#pragma omp simd collapse(2)
#endif
for (int dim = 0; dim < 3; dim++)
{
for (int ccomp = 0; ccomp < 2; ccomp++)
{
c_out[dim][ccomp] = c_in[dim][ccomp];
}
}
}
}

void copyUnAligned(fftw_complex* const c_out[3],
const fftw_complex* const c_in[3])
{
const std::size_t size { 2 * cDataSize * sizeof(double) };

for (int dim = 0; dim < 3; dim++)
{
memcpy(c_out[dim], c_in[dim], size);
}
}

int main()
{
#ifdef USE_OMP
omp_set_num_threads(NUM_THREADS);
#endif

std::vector cVec1(cDataSize);
std::vector cVec2(cDataSize);
fftw_complex* c_Arr1[3];
fftw_complex* c_Arr2[3];

for (int dim = 0; dim < 3; dim++)
{
c_Arr1[dim] = new fftw_complex[cDataSize];
c_Arr2[dim] = new fftw_complex[cDataSize];
}

for (int dim = 0; dim < 3; dim++)
{
for (int i = 0; i < cDataSize; i++)
{
for (int ccomp = 0; ccomp < 2; ccomp++)
{
c_Arr1[dim][ccomp] = 0.;
c_Arr2[dim][ccomp] = 0.;
}
}
}

auto start_time = high_resolution_clock::now();

fillUnAligned(c_Arr1, 1.);
// copyUnAligned(c_Arr2, c_Arr1);

auto end_time = high_resolution_clock::now();

auto total_elapsed_time =
duration_cast(end_time - start_time).count() / 1e6;

std::cout

Подробнее здесь: https://stackoverflow.com/questions/781 ... lelization