vel*_*s14 1 c++ parallel-processing optimization performance cuda
我想让 C++ 函数更快。我正在向您询问可能的方法。
我最多可以使用 32 个 OMP 线程。
我可以使用 NVidia GPU。
该函数的 MWE 为：
#include <iostream>
#include <complex>
#include <cmath>
typedef std::numeric_limits<double> dbl;
#define _USE_MATH_DEFINES
#include <omp.h>
const std::complex<double> I(0.0, 1.0); // imaginary unit, I*I = -1
std::complex<double> zero_imag (0.0, 0.0);
const int N_rs = 1500;
const int l_max = 70;
const int lmax = 70;
const int N_thetas = l_max + 1;
const int N_phis = 2 * l_max + 2;
const int N_ps = 600;
const int nphi = 2 * l_max + 2;
const double sqrt_of_2_over_pi = sqrt( 2.0 / M_PI );
void rtop(std::complex<double> * Psi_outer_spec,
std::complex<double> * Psi_outer_spec_plm,
double * BJ,
double * wrk,
std::complex<double> * wrk2,
double * ris_without_ends,
double * r_primes_without_ends,
double * weights_Lobatto_without_ends
)
{
int l, kk, kkk, m;
long int idx, idxx, idxxx;
// #pragma omp parallel for firstprivate (wrk2) private(l, kkk, idx, m, kk, idxx, idxxx) schedule(static)
// #pragma omp target teams distribute parallel for firstprivate(wrk2) private(l, kkk, idx, m, kk, idxx, idxxx)
for (int i = 0; i <= (N_ps - 1); i++) { // THIS IS THE BOTTLENECK !!!
std::complex<double> sum1 = std::complex<double> (0.0, 0.0); // each thread creates a sum1 on its own
for (l = 0; l <= lmax; l++) {
for (kkk = 0; kkk <= (N_rs-1); kkk++) {
idx = i * (N_rs*(l_max+1)) + kkk * (l_max+1) + l;
wrk2[kkk] = pow(-I, l) * BJ[idx] * wrk[kkk];
}
for (m = 0; m <= (nphi-1); m++) {
sum1 = zero_imag;
for (kk = 0; kk <= (N_rs-1); kk++) {
idxx = kk * (N_thetas*N_phis) + l * N_phis + m;
sum1 += Psi_outer_spec[idxx] * wrk2[kk];
}
idxxx = i * (N_thetas*N_phis) + l * N_phis + m;
Psi_outer_spec_plm[idxxx] = sum1 * sqrt_of_2_over_pi;
}
// END for m loop
}
// END for l loop
}
// END for i loop
}
int main() {
double * wrk = new double [N_rs];
std::complex<double> * wrk2 = new std::complex<double> [N_rs];
double * ris_without_ends = new double [N_rs];
double * r_primes_without_ends = new double [N_rs];
double * weights_Lobatto_without_ends = new double [N_rs];
double * BJ = new double [N_ps * N_rs * (l_max+1)];
std::complex<double> * Psi_outer_spec = new std::complex<double> [N_rs * N_thetas * N_phis];
std::complex<double> * Psi_outer_spec_plm = new std::complex<double> [N_ps * N_thetas * N_phis];
rtop(Psi_outer_spec, Psi_outer_spec_plm, BJ, wrk, wrk2, ris_without_ends, r_primes_without_ends, weights_Lobatto_without_ends);
return 0;
}
关联的 CMakeLists.txt 是：
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
set(CMAKE_VERBOSE_MAKEFILE ON)
set(CMAKE_C_COMPILER "gcc")
set(CMAKE_CXX_COMPILER "g++")
project(trial)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -pedantic -Wall")
find_package(OpenMP)
add_executable(trial trial.cpp)
if(OpenMP_CXX_FOUND)
target_link_libraries(trial PUBLIC OpenMP::OpenMP_CXX)
endif()
set_property(TARGET trial PROPERTY CXX_STANDARD 17)
编译为：$ cmake .. 那么$ cmake --build . --config Release.
我的输出是：
-- The C compiler identification is GNU 11.3.0
-- The CXX compiler identification is GNU 11.3.0
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: /apps20/sw/eb/software/GCCcore/11.3.0/bin/gcc - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /apps20/sw/eb/software/GCCcore/11.3.0/bin/g++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Found OpenMP_C: -fopenmp (found version "4.5")
-- Found OpenMP_CXX: -fopenmp (found version "4.5")
-- Found OpenMP: TRUE (found version "4.5")
-- Configuring done
-- Generating done
-- Build files have been written to: /work4/clf/ouatu/trial_for_SO/build
然后进行构建：
[ 50%] Building CXX object CMakeFiles/trial.dir/trial.cpp.o
[100%] Linking CXX executable trial
[100%] Built target trial
我尝试过的：
通过 OpenMP 并行，我确实获得了加速。
我的 OpenMP GPU 卸载失败（看来我的编译器标志无法实现卸载）。CMakeLists.txt（这些标志在该 MWE 的显示中是隐藏的）
我愿意接受任何其他建议。
例如，rtop作为 CUDA 内核会受益吗？做到这一点很难吗？
谢谢你！
我建议使用经过一些优化和调整的 OpenMP 版本。快速回顾一些变化以及需要注意的事项：
整个事情wrk2[kkk] = pow(-I, l) * ... 是双重冗余的。首先，pow(-I, l)它是一种优雅但昂贵的表达 4 种不同值的方式。其次，它仅用作点积中的因子。您可以将整个结果折叠到最后的乘法中sum1 * sqrt_of_2_over_pi。这也允许wrk2实值化，这也将最内层循环从复数-复数点积转变为复数-实数点积。
像这样的多维索引计算idx = i * (N_rs*(l_max+1)) + kkk * (l_max+1) + l应该遵循Horner方法来完成，以避免冗余乘法。更多挑剔，但也更清晰。例如这里idx = (i * N_rs + kkk) * (l_max+1) + l。当我们这样做时，请小心您的索引变量。他们都是int。特别是 3 维数组的大小可以快速增长到多个 GiB，此时您将遇到整数溢出。std::ptrdiff_t如果您担心这可能会成为问题，请切换到。
BJ和的迭代顺序Psi_outer_spec_plm并不理想。如果可能，BJ应该交换两个内部维度以获得更好的数据局部性，这也将允许循环初始化的矢量化wrk2。Psi_outer_spec更糟糕的是，因为您在最内层循环中沿着外部维度进行迭代。然而，我认为选择这个顺序是为了让它与 with 相同Psi_outer_spec_plm，并且因为它是好的。无论如何，这种更高的步幅会阻止矢量化。
我不明白为什么您在使用计数器和索引变量的范围之外声明它们。即使现代 C 标准也允许在 for 循环内声明它们，更不用说 C++ 了。对于并行化，您需要限制共享或意外共享变量的数量。
说到共享数据，据我所知，线程可能重叠的唯一共享内存是数组wrk2。可以简单地为每个线程分配它，这将我们带入最终的实现。
# pragma omp parallel
{
auto wrk2 = std::make_unique<double[]>(N_rs);
# pragma omp for collapse(2) nowait
for (int i = 0; i <= (N_ps - 1); i++) {
for (int l = 0; l <= lmax; l++) {
for (int kkk = 0; kkk <= (N_rs-1); kkk++) {
int idx = (i * N_rs + kkk) * (lmax + 1) + l;
wrk2[kkk] = BJ[idx] * wrk[kkk];
}
constexpr std::complex<double> I(0., 1.);
std::complex<double> factor(-sqrt_of_2_over_pi);
if(l & 1)
factor *= I;
if(l & 2)
factor = -factor;
for (int m = 0; m <= (N_phis-1); m++) {
std::complex<double> sum1;
for (int kk = 0; kk <= (N_rs-1); kk++) {
int idx = (kk * N_thetas + l) * N_phis + m;
sum1 += Psi_outer_spec[idx] * wrk2[kk];
}
int idx = (i * N_thetas + l) * N_phis + m;
Psi_outer_spec_plm[idx] = sum1 * factor;
}
}
}
}
请注意如何将通常pragma omp parallel for分为omp parallel和 单独的omp for以允许分配临时内存。这collapse(2)意味着两个外循环都是并行的。
其他需要考虑的事项：
m循环更改为矩阵向量乘积，这可能会通过 BLAS 库解决一些向量化/内存访问问题-march=native或者您想要的任何基线架构，那么这里应该是值得的。-mavx2 -mfma可能是一个很好的折衷方案，可以处理所有相对较新的 CPU，而无需过多地专门化二进制文件回到将循环卸载m到矩阵向量乘积的想法，我们必须重新解释Psi_outer_spec我们用作矩阵的切片。我选择列主矩阵，因为我想在这一步中使用 Eigen3。
N_phi（循环计数器m）N_rs（循环计数器kk）N_phi * N_thetal * N_phis假设这是正确的，我们可以将数组映射到特征向量和矩阵，并让它处理转置访问。这会将wrk2初始化下面的所有内容都变成此代码
using MatrixMap = Eigen::Map<const Eigen::MatrixXcd,
Eigen::Unaligned, Eigen::OuterStride<>>;
MatrixMap Psi_slice(
Psi_outer_spec + l * N_phis /*top left corner*/,
N_phis /*rows*/, N_rs /*cols*/,
Eigen::OuterStride<>(N_phis * N_thetas));
const auto wrk2_mapped = Eigen::VectorXd::Map(wrk2.get(), N_rs);
auto Psi_plm_mapped = Eigen::VectorXcd::Map(
Psi_outer_spec_plm + (i * N_thetas + l) * N_phis, N_phis);
Psi_plm_mapped.noalias() = Psi_slice * wrk2_mapped * factor;
现在，这一步显然提出了一个问题，即我们是否可以通过一些预处理或后处理将整个事物转变为矩阵-矩阵乘积，这可能会处理整个并行化和潜在的 GPU 卸载。这就是为什么我要求提供数学描述，而不是通过代码进行这种徒劳的追逐
确实可以将其重写为矩阵-矩阵乘积。Psi_outer_spec诀窍是独立于 的观察i。因此，如果我们切换两个外部循环，我们就可以在一次操作中计算出其中一个的所有l值i。
在这样做的同时，我又回到了wrk2复杂化并包含了这个因素。从技术上讲，这需要更多的计算时间和内存，但对于矩阵-矩阵产品，您可能希望直接分派到 BLAS 后端，例如通过OpenBLAS、通过Eigen 的后端甚至 GPU 加速（例如CuBLAS）。为此，您需要复数乘法。
Eigen::MatrixXcd wrk2mat(N_rs, N_ps);
for (int l = 0; l <= lmax; l++) {
std::complex<double> factor(-sqrt_of_2_over_pi);
if(l & 1)
factor *= I;
if(l & 2)
factor = -factor;
# pragma omp parallel for
for (int i = 0; i <= N_ps - 1; i++) {
for (int k = 0; k <= N_rs - 1; ++k) {
int idx = (i * N_rs + k) * (lmax + 1) + l;
wrk2mat(k, i) = BJ[idx] * wrk[k] * factor;
}
}
using ConstMatrixMap = Eigen::Map<const Eigen::MatrixXcd,
Eigen::Unaligned, Eigen::OuterStride<>>;
ConstMatrixMap Psi_slice(
Psi_outer_spec + l * N_phis /*top left corner*/,
N_phis /*rows*/, N_rs /*cols*/,
Eigen::OuterStride<>(N_phis * N_thetas));
using MatrixMap = Eigen::Map<Eigen::MatrixXcd,
Eigen::Unaligned, Eigen::OuterStride<>>;
MatrixMap Psi_plm_mapped(
Psi_outer_spec_plm + l * N_phis,
N_phis, N_ps,
Eigen::OuterStride<>((lmax + 1) * N_phis));
Psi_plm_mapped.noalias() = Psi_slice * wrk2mat;
}
只要矩阵足够大，矩阵-矩阵乘积就应该在内部并行化。如果情况并非总是如此，您可以将整个事情包装到运行时可选的并行块中。大致是这样的：
bool small_matrices = ...;
#pragma omp parallel if(small_matrices)
{
Eigen::MatrixXcd wrk2mat(N_rs, N_ps);
# pragma omp for nowait
for (int l = 0; l <= lmax; l++) {
...
}
}
由于 OpenMP 通常会停用嵌套并行化，因此这将自动停用所有内部parallel部分并按顺序运行它们。