c++ - 使用 OpenMP 线程和 std::(experimental::)simd 计算 Mandelbrot 集
问题描述
我希望使用不同类型的 HPC 范例来实现一个简单的 Mandelbrot 集绘图仪,展示它们的优缺点以及实现的难易程度。想想 GPGPU (CUDA/OpenACC/OpenMP4.5)、线程/OpenMP 和 MPI。并使用这些示例为 HPC 新手提供帮助并了解可能性。代码的清晰性比从硬件中获得绝对顶级性能更重要,这是第二步;)
因为并行化这个问题很简单,而且现代 CPU 可以使用向量指令获得大量性能,所以我还想结合 OpenMP 和 SIMD。不幸的是,简单地添加 a#pragma omp simd
并不能产生令人满意的结果,并且使用内在函数不是非常用户友好或未来证明。还是蛮漂亮的。
幸运的是,正在对 C++ 标准进行工作,以便更容易通用地实现向量指令,如 TS: "Extensions for parallelism, version 2"中所述,特别是第 9 节关于数据并行类型。WIP 实现可以在这里找到,它基于 VC,可以在这里找到。
假设我有以下课程(已更改以使其更简单)
#include <stddef.h>
using Range = std::pair<double, double>;
using Resolution = std::pair<std::size_t, std::size_t>;
class Mandelbrot
{
double* d_iters;
Range d_xrange;
Range d_yrange;
Resolution d_res;
std::size_t d_maxIter;
public:
Mandelbrot(Range xrange, Range yrange, Resolution res, std::size_t maxIter);
~Mandelbrot();
void writeImage(std::string const& fileName);
void computeMandelbrot();
private:
void calculateColors();
};
以及以下computeMandelbrot()
使用 OpenMP的实现
void Mandelbrot::computeMandelbrot()
{
double dx = (d_xrange.second - d_xrange.first) / d_res.first;
double dy = (d_yrange.second - d_yrange.first) / d_res.second;
#pragma omp parallel for schedule(dynamic)
for (std::size_t row = 0; row != d_res.second; ++row)
{
double c_imag = d_yrange.first + row * dy;
for (std::size_t col = 0; col != d_res.first; ++col)
{
double real = 0.0;
double imag = 0.0;
double realSquared = 0.0;
double imagSquared = 0.0;
double c_real = d_xrange.first + col * dx;
std::size_t iter = 0;
while (iter < d_maxIter && realSquared + imagSquared < 4.0)
{
realSquared = real * real;
imagSquared = imag * imag;
imag = 2 * real * imag + c_imag;
real = realSquared - imagSquared + c_real;
++iter;
}
d_iters[row * d_res.first + col] = iter;
}
}
}
我们可以假设 x 和 y 方向的分辨率都是 2/4/8/.. 的倍数,具体取决于我们使用的 SIMD 指令。
不幸的是,在std::experimental::simd
. 据我所知,也没有任何重要的例子。
在 Vc git 存储库中,有一个 Mandelbrot 集合计算器的实现,但它非常复杂,并且由于缺少注释而很难理解。
很明显,我应该更改函数中双精度数的数据类型computeMandelbrot()
,但我不确定是什么。TS 提到了某些类型 T 的两个主要新数据类型,
native_simd = std::experimental::simd<T, std::experimental::simd_abi::native>;
和
fixed_size_simd = std::experimental::simd<T, std::experimental::simd_abi::fixed_size<N>>;
使用native_simd
是最有意义的,因为我在编译时不知道我的界限。但是我不清楚这些类型代表什么,是native_simd<double>
单个双精度数还是执行向量指令的双精度数集合?那么这个系列中有多少双打?
如果有人可以指出使用这些概念的示例,或者给我一些关于如何使用 std::experimental::simd 实现矢量指令的指示,我将非常感激。
解决方案
这是一个非常基本的实现,它有效(据我所知)。测试向量的哪些元素的绝对值大于 2 的方法非常繁琐且效率低下。必须有更好的方法来做到这一点,但我还没有找到它。
我在 AMD Ryzen 5 3600 上获得了大约 72% 的性能提升,并为 g++ 提供了选项-march=znver2
,这低于预期。
template <class T>
void mandelbrot(T xstart, T xend,
T ystart, T yend)
{
namespace stdx = std::experimental;
constexpr auto simdSize = stdx::native_simd<T>().size();
constexpr unsigned size = 4096;
constexpr unsigned maxIter = 250;
assert(size % simdSize == 0);
unsigned* res = new unsigned[size * size];
T dx = (xend - xstart) / size;
T dy = (yend - ystart) / size;
for (std::size_t row = 0; row != size; ++row)
{
T c_imag = ystart + row * dy;
for (std::size_t col = 0; col != size; col += simdSize)
{
stdx::native_simd<T> real{0};
stdx::native_simd<T> imag{0};
stdx::native_simd<T> realSquared{0};
stdx::native_simd<T> imagSquared{0};
stdx::fixed_size_simd<unsigned, simdSize> iters{0};
stdx::native_simd<T> c_real;
for (int idx = 0; idx != simdSize; ++idx)
{
c_real[idx] = xstart + (col + idx) * dx;
}
for (unsigned iter = 0; iter != maxIter; ++iter)
{
realSquared = real * real;
imagSquared = imag * imag;
auto isInside = realSquared + imagSquared > stdx::native_simd<T>{4};
for (int idx = 0; idx != simdSize; ++idx)
{
// if not bigger than 4, increase iters
if (!isInside[idx])
{
iters[idx] += 1;
}
else
{
// prevent that they become inf/nan
real[idx] = static_cast<T>(4);
imag[idx] = static_cast<T>(4);
}
}
if (stdx::all_of(isInside) )
{
break;
}
imag = static_cast<T>(2.0) * real * imag + c_imag;
real = realSquared - imagSquared + c_real;
}
iters.copy_to(res + row * size + col, stdx::element_aligned);
}
}
delete[] res;
}
整个测试代码(从 开始auto test = (...)
)编译为
.L9:
vmulps ymm1, ymm1, ymm1
vmulps ymm13, ymm2, ymm2
xor eax, eax
vaddps ymm2, ymm13, ymm1
vcmpltps ymm2, ymm5, ymm2
vmovaps YMMWORD PTR [rsp+160], ymm2
jmp .L6
.L3:
vmovss DWORD PTR [rsp+32+rax], xmm0
vmovss DWORD PTR [rsp+64+rax], xmm0
add rax, 4
cmp rax, 32
je .L22
.L6:
vucomiss xmm3, DWORD PTR [rsp+160+rax]
jp .L3
jne .L3
inc DWORD PTR [rsp+96+rax]
add rax, 4
cmp rax, 32
jne .L6
推荐阅读
- python - 使用 Python 获取 SVN 存储库的根目录
- ansible - 如何遍历列表的字典并将变量传递给 Ansible 角色?
- windows - 在 Buildbot 中设置和获取属性
- xamarin - Xamarin 表单 - 动态更改网格行可见性
- javascript - VueJS 首先导航到根路由
- google-chrome - Chrome 扩展:无法监听通过服务工作者进行的网络调用
- javascript - VueJs - Retrieve objects from array into new object array
- android - 防止 Jacoco 在默认 gradle 任务中运行(在每个本地构建中都会发生)
- ruby - Rspec mongoid - 测试嵌入式文档回调(after_save)
- c# - 将存储在项目文件夹中的图像添加到 Microsoft Interop Word 文档