如何使用 Eigen 在 C++ 中多线程处理此代码片段
How can I multithread this code snippet in C++ with Eigen
我正在尝试实现以下代码片段的更快版本：
Eigen::VectorXd dTX = (( (XPSF.array() - x0).square() + (ZPSF.array() - z0).square() ).sqrt() + txShift)*fs/c + t0*fs;
Eigen::VectorXd Zsq = ZPSF.array().square();
Eigen::MatrixXd idxt(XPSF.size(),nc);
for (int i = 0; i < nc; i++) {
idxt.col(i) = ((XPSF.array() - xe(i)).square() + Zsq.array()).sqrt()*fs/c + dTX.array();
idxt.col(i) = (abs(XPSF.array()-xe(i)) <= ZPSF.array()*0.5/fnumber).select(idxt.col(i),-1);
}
我现在使用的示例数组大小是：
XPSF：591*192个系数的列向量（列向量中总共有113,472个值）
ZPSF：与 XPSF 大小相同
xe：192个系数的RowVector
idxt：113,472x192 大小的矩阵
当前使用 gcc 和 -msse2 和 -o3 优化运行时，循环第一行的平均时间约为 0.08 秒，循环第二行的平均时间约为 0.03 秒。我知道运行时依赖于平台，但我相信这仍然可以更快。商业软件执行我在这里尝试做的操作的时间减少了两个数量级。另外，我怀疑我的代码现在有点业余！
我尝试阅读 Eigen 文档以了解矢量化的工作原理、实现位置以及此代码中有多少可能被 Eigen“隐式”并行化，但我一直在努力跟踪细节。总的来说，我对 C++ 也有点陌生，但我看过有关 std::thread 的文档和其他资源，并尝试将其与此代码结合使用，但没有取得太大成功。
如有任何建议，我们将不胜感激。
更新：
更新 2
我会赞成 Soleil 的回答，因为如果我有它的声誉分数，它包含有用的信息。但是，我应该澄清一下，我想首先弄清楚没有 GPU 我可以做哪些优化。我确信（尽管没有 OpenMP）Eigen 固有的多线程和矢量化不会进一步加快它的速度（除非生成不必要的临时文件）。我怎么能使用 std::thread 之类的东西来明确地并行化呢？为此，我正在努力将 std::thread 和 Eigen 结合起来。
OpenMP
如果您的 CPU 有足够多的内核和线程，通常简单快速的第一步是通过添加 pragma 来调用 OpenMP：
#pragma omp parallel for
for (int i = 0; i < nc; i++)
并使用 /openmp (cl) 或 -fopenmp (gcc) 进行编译，或者仅使用 gcc -ftree-parallelize-loops 进行编译，以便自动展开循环。
这将执行 map reduce，并且 map 将出现在您的 CPU 可以处理的并行线程数上（7700HQ 为 8 个线程）。
通常您也可以 set a clause num_threads(n) 其中 n 是所需的线程数：
#pragma omp parallel num_threads(8)
我使用 8 的地方是因为 7700HQ 可以处理 8 个并发线程。
TBB
您也可以 unroll your loop with TBB:
#pragma unroll
for (int i = 0; i < nc; i++)
与特征集成的线程
有了Eigen你can add
OMP_NUM_THREADS=n ./my_program
omp_set_num_threads(n);
Eigen::setNbThreads(n);
带有特征的多线程的备注
然而，在 FAQ:
currently Eigen parallelizes only general matrix-matrix products (bench), so it doesn't by itself take much advantage of parallel hardware."
一般来说，OpenMP 的改进并不总是存在，因此请对发布版本进行基准测试。另一种方法是确保您使用的是矢量化指令。
同样，来自 FAQ/vectorization：
How can I enable vectorization?
You just need to tell your compiler to enable the corresponding
instruction set, and Eigen will then detect it. If it is enabled by
default, then you don't need to do anything. On GCC and clang you can
simply pass -march=native to let the compiler enables all instruction
set that are supported by your CPU.
On the x86 architecture, SSE is not enabled by default by most
compilers. You need to enable SSE2 (or newer) manually. For example,
with GCC, you would pass the -msse2 command-line option.
On the x86-64 architecture, SSE2 is generally enabled by default, but
you can enable AVX and FMA for better performance
On PowerPC, you have to use the following flags: -maltivec
-mabi=altivec, for AltiVec, or -mvsx for VSX-capable systems.
On 32-bit ARM NEON, the following: -mfpu=neon -mfloat-abi=softfp|hard,
depending if you are on a softfp/hardfp system. Most current
distributions are using a hard floating-point ABI, so go for the
latter, or just leave the default and just pass -mfpu=neon.
On 64-bit ARM, SIMD is enabled by default, you don't have to do
anything extra.
On S390X SIMD (ZVector), you have to use a recent gcc (version >5.2.1)
compiler, and add the following flags: -march=z13 -mzvector.
使用 cuda 进行多线程处理
鉴于您的数组大小，您想尝试卸载到 GPU 以达到微秒级；在那种情况下，您将拥有（通常）与数组中的元素数量一样多的线程。
对于一个简单的开始，如果你有一个 nvidia 卡，你想看看 cublas，它也允许你使用上一代的张量寄存器（融合乘加等），不像常规内核。
因为 eigen 是一个只有头文件的库，所以 you could use it in a cuda kernel.
是有道理的
您也可以使用常规内核“手动”（即，没有本征）实现所有内容。这是工程上的废话，但是在一个education/university项目中的常见做法，才能看得懂。
使用 OneAPI 和英特尔 GPU 的多线程
既然你有一个skylake架构，你也可以unroll your loop on your CPU's GPU with OneAPI:
// Unroll loop as specified by the unroll factor.
#pragma unroll unroll_factor
for (int i = 0; i < nc; i++)
(来自sample).
我正在尝试实现以下代码片段的更快版本：
Eigen::VectorXd dTX = (( (XPSF.array() - x0).square() + (ZPSF.array() - z0).square() ).sqrt() + txShift)*fs/c + t0*fs;
Eigen::VectorXd Zsq = ZPSF.array().square();
Eigen::MatrixXd idxt(XPSF.size(),nc);
for (int i = 0; i < nc; i++) {
idxt.col(i) = ((XPSF.array() - xe(i)).square() + Zsq.array()).sqrt()*fs/c + dTX.array();
idxt.col(i) = (abs(XPSF.array()-xe(i)) <= ZPSF.array()*0.5/fnumber).select(idxt.col(i),-1);
}
我现在使用的示例数组大小是：
XPSF：591*192个系数的列向量（列向量中总共有113,472个值）
ZPSF：与 XPSF 大小相同
xe：192个系数的RowVector
idxt：113,472x192 大小的矩阵
当前使用 gcc 和 -msse2 和 -o3 优化运行时，循环第一行的平均时间约为 0.08 秒，循环第二行的平均时间约为 0.03 秒。我知道运行时依赖于平台，但我相信这仍然可以更快。商业软件执行我在这里尝试做的操作的时间减少了两个数量级。另外，我怀疑我的代码现在有点业余！
我尝试阅读 Eigen 文档以了解矢量化的工作原理、实现位置以及此代码中有多少可能被 Eigen“隐式”并行化，但我一直在努力跟踪细节。总的来说，我对 C++ 也有点陌生，但我看过有关 std::thread 的文档和其他资源，并尝试将其与此代码结合使用，但没有取得太大成功。
如有任何建议，我们将不胜感激。
更新：
更新 2
我会赞成 Soleil 的回答，因为如果我有它的声誉分数，它包含有用的信息。但是，我应该澄清一下，我想首先弄清楚没有 GPU 我可以做哪些优化。我确信（尽管没有 OpenMP）Eigen 固有的多线程和矢量化不会进一步加快它的速度（除非生成不必要的临时文件）。我怎么能使用 std::thread 之类的东西来明确地并行化呢？为此，我正在努力将 std::thread 和 Eigen 结合起来。
OpenMP
如果您的 CPU 有足够多的内核和线程，通常简单快速的第一步是通过添加 pragma 来调用 OpenMP：
#pragma omp parallel for
for (int i = 0; i < nc; i++)
并使用 /openmp (cl) 或 -fopenmp (gcc) 进行编译，或者仅使用 gcc -ftree-parallelize-loops 进行编译，以便自动展开循环。
这将执行 map reduce，并且 map 将出现在您的 CPU 可以处理的并行线程数上（7700HQ 为 8 个线程）。
通常您也可以 set a clause num_threads(n) 其中 n 是所需的线程数：
#pragma omp parallel num_threads(8)
我使用 8 的地方是因为 7700HQ 可以处理 8 个并发线程。
TBB
您也可以 unroll your loop with TBB:
#pragma unroll
for (int i = 0; i < nc; i++)
与特征集成的线程
有了Eigen你can add
OMP_NUM_THREADS=n ./my_program
omp_set_num_threads(n);
Eigen::setNbThreads(n);
带有特征的多线程的备注
然而，在 FAQ:
currently Eigen parallelizes only general matrix-matrix products (bench), so it doesn't by itself take much advantage of parallel hardware."
一般来说，OpenMP 的改进并不总是存在，因此请对发布版本进行基准测试。另一种方法是确保您使用的是矢量化指令。
同样，来自 FAQ/vectorization：
How can I enable vectorization?
You just need to tell your compiler to enable the corresponding instruction set, and Eigen will then detect it. If it is enabled by default, then you don't need to do anything. On GCC and clang you can simply pass -march=native to let the compiler enables all instruction set that are supported by your CPU.
On the x86 architecture, SSE is not enabled by default by most compilers. You need to enable SSE2 (or newer) manually. For example, with GCC, you would pass the -msse2 command-line option.
On the x86-64 architecture, SSE2 is generally enabled by default, but you can enable AVX and FMA for better performance
On PowerPC, you have to use the following flags: -maltivec -mabi=altivec, for AltiVec, or -mvsx for VSX-capable systems.
On 32-bit ARM NEON, the following: -mfpu=neon -mfloat-abi=softfp|hard, depending if you are on a softfp/hardfp system. Most current distributions are using a hard floating-point ABI, so go for the latter, or just leave the default and just pass -mfpu=neon.
On 64-bit ARM, SIMD is enabled by default, you don't have to do anything extra.
On S390X SIMD (ZVector), you have to use a recent gcc (version >5.2.1) compiler, and add the following flags: -march=z13 -mzvector.
使用 cuda 进行多线程处理
鉴于您的数组大小，您想尝试卸载到 GPU 以达到微秒级；在那种情况下，您将拥有（通常）与数组中的元素数量一样多的线程。
对于一个简单的开始，如果你有一个 nvidia 卡，你想看看 cublas，它也允许你使用上一代的张量寄存器（融合乘加等），不像常规内核。
因为 eigen 是一个只有头文件的库，所以 you could use it in a cuda kernel.
是有道理的您也可以使用常规内核“手动”（即，没有本征）实现所有内容。这是工程上的废话，但是在一个education/university项目中的常见做法，才能看得懂。
使用 OneAPI 和英特尔 GPU 的多线程
既然你有一个skylake架构，你也可以unroll your loop on your CPU's GPU with OneAPI:
// Unroll loop as specified by the unroll factor.
#pragma unroll unroll_factor
for (int i = 0; i < nc; i++)
(来自sample).