如何解决使用OpenMP线程和std::( experimentalal::) simd计算Mandelbrot集
我希望使用不同类型的HPC范例来实现一个简单的Mandelbrot集绘图仪,以显示它们的优缺点以及实现的难易程度。考虑一下GPGPU(CUDA / OpenACC / OpenMP4.5),线程/ OpenMP和MPI。并使用这些示例为刚接触HPC的程序员提供帮助,并了解可能的情况。代码的清晰性比从硬件中获得绝对最高的性能更为重要,这是第二步;)
由于该问题很难并行化,并且现代CPU可以使用矢量指令获得大量性能,因此我也想将OpenMP和SIMD结合起来。不幸的是,简单地添加#pragma omp simd
并不能产生令人满意的结果,并且使用内在函数不是非常用户友好或未来的证明。或pretty。
幸运的是,正在按照C ++标准进行工作,以便更容易地通用实现矢量指令,如TS:"Extensions for parallelism,version 2"中所述,尤其是关于数据并行类型的第9节。可以在here的基础上找到here的WIP实现。
假设我有以下课程(已更改为使其更简单一些)
#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()
中更改double的数据类型,但是不确定。 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,32
jne .L6
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。