使用 Rcpp 和 OpenMP 在 R 中多线程和 SIMD 矢量化 Mandelbrot
Posted
技术标签:
【中文标题】使用 Rcpp 和 OpenMP 在 R 中多线程和 SIMD 矢量化 Mandelbrot【英文标题】：Multithreaded & SIMD vectorized Mandelbrot in R using Rcpp & OpenMP 【发布时间】：2018-06-12 16:35:57 【问题描述】：作为OpenMP 和Rcpp 性能测试，我想检查使用最直接最简单的Rcpp+OpenMP 实现在R 中计算曼德布罗集的速度有多快。目前我所做的是：
#include <Rcpp.h>
#include <omp.h>
// [[Rcpp::plugins(openmp)]]
using namespace Rcpp;
// [[Rcpp::export]]
Rcpp::NumericMatrix mandelRcpp(const double x_min, const double x_max, const double y_min, const double y_max,
const int res_x, const int res_y, const int nb_iter)
Rcpp::NumericMatrix ret(res_x, res_y);
double x_step = (x_max - x_min) / res_x;
double y_step = (y_max - y_min) / res_y;
int r,c;
#pragma omp parallel for default(shared) private(c) schedule(dynamic,1)
for (r = 0; r < res_y; r++)
for (c = 0; c < res_x; c++)
double zx = 0.0, zy = 0.0, new_zx;
double cx = x_min + c*x_step, cy = y_min + r*y_step;
int n = 0;
for (n=0; (zx*zx + zy*zy < 4.0 ) && ( n < nb_iter ); n++ )
new_zx = zx*zx - zy*zy + cx;
zy = 2.0*zx*zy + cy;
zx = new_zx;
ret(c,r) = n;
return ret;
然后在 R 中：
library(Rcpp)
sourceCpp("mandelRcpp.cpp")
xlims=c(-0.74877,-0.74872);
ylims=c(0.065053,0.065103);
x_res=y_res=1080L; nb_iter=10000L;
system.time(m <- mandelRcpp(xlims[[1]], xlims[[2]], ylims[[1]], ylims[[2]], x_res, y_res, nb_iter))
# 0.92s
rainbow=c(rgb(0.47,0.11,0.53),rgb(0.27,0.18,0.73),rgb(0.25,0.39,0.81),rgb(0.30,0.57,0.75),rgb(0.39,0.67,0.60),rgb(0.51,0.73,0.44),rgb(0.67,0.74,0.32),rgb(0.81,0.71,0.26),rgb(0.89,0.60,0.22),rgb(0.89,0.39,0.18),rgb(0.86,0.13,0.13))
cols=c(colorRampPalette(rainbow)(100),rev(colorRampPalette(rainbow)(100)),"black") # palette
par(mar=c(0, 0, 0, 0))
system.time(image(m^(1/7), col=cols, asp=diff(ylims)/diff(xlims), axes=F, useRaster=T))
# 0.5s
我不确定除了 OpenMP 多线程之外是否还有其他明显的速度改进可以利用，例如通过simd 矢量化？ （在 openmp 中使用 simd 选项 #pragma 似乎没有做任何事情）
PS 一开始我的代码崩溃了，但后来我发现这是通过将ret[r,c] = n; 替换为ret(r,c) = n; 来解决的
使用下面答案中建议的犰狳类会使事情变得稍微快一些，尽管时间几乎相同。还翻转了x 和y，因此当使用image() 绘制时，它以正确的方向出现。使用 8 个线程的速度是 ca。比矢量化纯 R Mandelbrot 版本 here 快 350 倍，也比（非多线程）Python/Numba 版本 here 快约 7.3 倍（类似于 PyCUDA 或 PyOpenCL 速度），对此非常满意...... Rasterizing/display now seems the bottleneck in R....
【问题讨论】：
一般来说，我通过避免在相同轮廓区域和 M 集上的迭代来提高速度（C 与汇编程序迭代）。远离 M-Set 边界，轮廓内包含大面积区域，我开发了一种曲线缝合方法来遵循轮廓边界，然后将其填充。迭代越深，增益越好。当一个芽被意外剪断时可能会受到惩罚，而且我看不出这种方法在使用线程时会如何工作。在进行两倍缩放时可以找到另一个节省，其中 1/4 的点是已知的。 是的，但另一方面，我计划转向连续着色，其中第一种优化不再那么简单了。重用我打算在缩放时计算的像素...在像这样的高缩放策略en.wikipedia.org/wiki/… 可以极大地提高性能。但我的主要问题更多地集中在我的 Rcpp 代码上，而不是更多地关注可以做的进一步算法优化，当然这些优化很多......而在 R 中，主要瓶颈似乎只是显示 我从来没有用颜色填充轮廓区域，只有迭代。着色算法是另一回事。 不是真的，因为一个人不再使用简单的转义时间算法，一个人没有得到连续的数字，而不是固定的迭代次数，正如en.wikipedia.org/wiki/…中所解释的@ 查看这里的 Python 代码示例：ibm.com/developerworks/community/blogs/jfp/entry/… 的两种方法... 【参考方案1】：我继续使用 GCC 和 Clang 的矢量扩展对 OP 的代码进行矢量化。在展示我是如何做到这一点之前，让我展示一下使用以下硬件的性能：
Skylake (SKL) at 3.1 GHz with 4 cores
Knights Landing (KNL) at 1.5 GHz with 68 cores
ARMv8 Cortex-A57 arch64 (Nvidia Jetson TX1) 4 cores at ? GHz
nb_iter = 1000000
GCC Clang
SKL_scalar 6m5,422s
SKL_SSE41 3m18,058s
SKL_AVX2 1m37,843s 1m39,943s
SKL_scalar_omp 0m52,237s
SKL_SSE41_omp 0m29,624s 0m31,356s
SKL_AVX2_omp 0m14,156s 0m16,783s
ARM_scalar 15m28.285s
ARM_vector 9m26.384s
ARM_scalar_omp 3m54.242s
ARM_vector_omp 2m21.780s
KNL_scalar 19m34.121s
KNL_SSE41 11m30.280s
KNL_AVX2 5m0.005s 6m39.568s
KNL_AVX512 2m40.934s 6m20.061s
KNL_scalar_omp 0m9.108s
KNL_SSE41_omp 0m6.666s 0m6.992s
KNL_AVX2_omp 0m2.973s 0m3.988s
KNL_AVX512_omp 0m1.761s 0m3.335s
KNL 与 SKL 的理论加速是
(68 cores/4 cores)*(1.5 GHz/3.1 Ghz)*
(8 doubles per lane/4 doubles per lane) = 16.45
我详细介绍了 GCC 和 Clang 的矢量扩展功能 here。为了向量化 OP 的代码，我们需要定义三个额外的向量操作。
1.广播
对于向量 v 和标量 s GCC 不能做 v = s 但 Clang 可以。但我找到了一个适用于 GCC 和 Clang here 的不错的解决方案。例如
vsi v = s - (vsi);
2。一个any() 函数like in OpenCL 或类似R。
我想出的最好的是通用函数
static bool any(vli const & x)
for(int i=0; i<VLI_SIZE; i++) if(x[i]) return true;
return false;
Clang 实际上使用 ptest 指令（但 not for AVX512）为此相对生成 efficient code，但 GCC 不会。
3。压缩
计算以 64 位双精度数形式完成，但结果以 32 位整数形式写出。因此，使用 64 位整数完成两次计算，然后将两次计算压缩为一个 32 位整数向量。我想出了一个通用的解决方案，Clang 做得很好
static vsi compress(vli const & lo, vli const & hi)
vsi lo2 = (vsi)lo, hi2 = (vsi)hi, z;
for(int i=0; i<VLI_SIZE; i++) z[i+0*VLI_SIZE] = lo2[2*i];
for(int i=0; i<VLI_SIZE; i++) z[i+1*VLI_SIZE] = hi2[2*i];
return z;
以下解决方案有效better for GCC but is no better for Clang。但由于这个函数并不重要，我只使用通用版本。
static vsi compress(vli const & low, vli const & high)
#if defined(__clang__)
return __builtin_shufflevector((vsi)low, (vsi)high, MASK);
#else
return __builtin_shuffle((vsi)low, (vsi)high, (vsi)MASK);
#endif
这些定义不依赖于任何特定于 x86 的内容，并且代码（定义如下）可以为 ARM 处理器以及 GCC 和 Clang 编译。
现在这里定义了这些就是代码
#include <string.h>
#include <inttypes.h>
#include <Rcpp.h>
using namespace Rcpp;
#ifdef _OPENMP
#include <omp.h>
#endif
// [[Rcpp::plugins(openmp)]]
// [[Rcpp::plugins(cpp14)]]
#if defined ( __AVX512F__ ) || defined ( __AVX512__ )
static const int SIMD_SIZE = 64;
#elif defined ( __AVX2__ )
static const int SIMD_SIZE = 32;
#else
static const int SIMD_SIZE = 16;
#endif
static const int VSI_SIZE = SIMD_SIZE/sizeof(int32_t);
static const int VLI_SIZE = SIMD_SIZE/sizeof(int64_t);
static const int VDF_SIZE = SIMD_SIZE/sizeof(double);
#if defined(__clang__)
typedef int32_t vsi __attribute__ ((ext_vector_type(VSI_SIZE)));
typedef int64_t vli __attribute__ ((ext_vector_type(VLI_SIZE)));
typedef double vdf __attribute__ ((ext_vector_type(VDF_SIZE)));
#else
typedef int32_t vsi __attribute__ ((vector_size (SIMD_SIZE)));
typedef int64_t vli __attribute__ ((vector_size (SIMD_SIZE)));
typedef double vdf __attribute__ ((vector_size (SIMD_SIZE)));
#endif
static bool any(vli const & x)
for(int i=0; i<VLI_SIZE; i++) if(x[i]) return true;
return false;
static vsi compress(vli const & lo, vli const & hi)
vsi lo2 = (vsi)lo, hi2 = (vsi)hi, z;
for(int i=0; i<VLI_SIZE; i++) z[i+0*VLI_SIZE] = lo2[2*i];
for(int i=0; i<VLI_SIZE; i++) z[i+1*VLI_SIZE] = hi2[2*i];
return z;
// [[Rcpp::export]]
IntegerVector frac(double x_min, double x_max, double y_min, double y_max, int res_x, int res_y, int nb_iter)
IntegerVector out(res_x*res_y);
vdf x_minv = x_min - (vdf), y_minv = y_min - (vdf);
vdf x_stepv = (x_max - x_min)/res_x - (vdf), y_stepv = (y_max - y_min)/res_y - (vdf);
double a[VDF_SIZE] __attribute__ ((aligned(SIMD_SIZE)));
for(int i=0; i<VDF_SIZE; i++) a[i] = 1.0*i;
vdf vi0 = *(vdf*)a;
#pragma omp parallel for schedule(dynamic) collapse(2)
for (int r = 0; r < res_y; r++)
for (int c = 0; c < res_x/(VSI_SIZE); c++)
vli nv[2] = 0 - (vli), 0 - (vli);
for(int j=0; j<2; j++)
vdf c2 = 1.0*VDF_SIZE*(2*c+j) + vi0;
vdf zx = 0.0 - (vdf), zy = 0.0 - (vdf), new_zx;
vdf cx = x_minv + c2*x_stepv, cy = y_minv + r*y_stepv;
vli t = -1 - (vli);
for (int n = 0; any(t = zx*zx + zy*zy < 4.0) && n < nb_iter; n++, nv[j] -= t)
new_zx = zx*zx - zy*zy + cx;
zy = 2.0*zx*zy + cy;
zx = new_zx;
vsi sp = compress(nv[0], nv[1]);
memcpy(&out[r*res_x + VSI_SIZE*c], (int*)&sp, SIMD_SIZE);
return out;
R代码和OP的代码几乎一样
library(Rcpp)
sourceCpp("frac.cpp", verbose=TRUE, rebuild=TRUE)
xlims=c(-0.74877,-0.74872);
ylims=c(0.065053,0.065103);
x_res=y_res=1080L; nb_iter=100000L;
t = system.time(m <- frac(xlims[[1]], xlims[[2]], ylims[[1]], ylims[[2]], x_res, y_res, nb_iter))
print(t)
m2 = matrix(m, ncol = x_res)
rainbow = c(
rgb(0.47, 0.11, 0.53),
rgb(0.27, 0.18, 0.73),
rgb(0.25, 0.39, 0.81),
rgb(0.30, 0.57, 0.75),
rgb(0.39, 0.67, 0.60),
rgb(0.51, 0.73, 0.44),
rgb(0.67, 0.74, 0.32),
rgb(0.81, 0.71, 0.26),
rgb(0.89, 0.60, 0.22),
rgb(0.89, 0.39, 0.18),
rgb(0.86, 0.13, 0.13)
)
cols = c(colorRampPalette(rainbow)(100),
rev(colorRampPalette(rainbow)(100)),"black") # palette
par(mar = c(0, 0, 0, 0))
image(m2^(1/7), col=cols, asp=diff(ylims)/diff(xlims), axes=F, useRaster=T)
要为 GCC 或 Clang 编译，请将文件 ~/.R/Makevars 更改为
CXXFLAGS= -Wall -std=c++14 -O3 -march=native -ffp-contract=fast -fopenmp
#uncomment the following two lines for clang
#CXX=clang-5.0
#LDFLAGS= -lomp
如果您无法让 OpenMP 为 Clang 工作，请参阅 this。
代码产生或多或少相同的图像。
【讨论】：
感谢数百万的详细教程 - 非常有启发性和优雅！我还刚刚发现 github.com/bisqwit/cpp_parallelization_examples/blob/master/… youtube.com/watch?v=Pc8DfEyAxzg youtube.com/watch?v=MfEkOcMILDo youtube.com/watch?v=pCoxpKTmykA&t=232s 在算法方面仍有一些优化，但我认为在代码方面没有那么优雅...... @TomWenseleers 我不想要被接受的投票。能不能还给原来的人。 @TomWenseleers，这些链接和视频太棒了！我们所做的工作显然有很多重叠之处。实际上，我几年前用内在函数编写了这个，并将其放入 SDL 中进行实时渲染。我还为 GPU 的 OpenCL 编写了它。此外，我已经为双双计算实现了它以增加进动，因为即使使用双倍放大，您也会很快耗尽分辨率。***.com/questions/30573443/… @TomWenseleers 我主要写了这个答案，因为我一直想测试向量扩展。除了all 函数和compress 函数，至少对于x86 的性能来说，内部函数可能不是必需的，但显式矢量化是必需的。我昨天第一次在 ARM 上测试我的代码。我没有看过 ARM 程序集，但至少我在 ARM 上得到了不错的加速。
谢谢，听起来真的很棒很有趣！对于颜色，这只是一个简单的伽马颜色变换，以稍微均衡颜色 - 使用的最佳伽马系数可能会有所不同。为了避免这种情况，我最后改用直方图均衡en.wikipedia.org/wiki/Histogram_equalization，因为它总是返回令人愉悦的颜色渐变...使用fractalforums.com/fractal-exteme/… 或ibm.com/developerworks/community/blogs/jfp/entry/… 中的平滑阴影也不错。【参考方案2】：
不要将 OpenMP 与 Rcpp 的 *Vector 或 *Matrix 对象一起使用，因为它们掩盖了 SEXP 函数/内存分配是单线程的。 OpenMP 是multi-threaded approach。
这就是代码崩溃的原因。
解决此限制的一种方法是使用非R 数据结构来存储结果。以下之一就足够了：arma::mat 或 Eigen::MatrixXd 或 std::vector<T>... 由于我喜欢犰狳，我将把 res 矩阵从 Rcpp::NumericMatrix 更改为 arma::mat。因此，以下代码将并行执行您的代码：
#include <RcppArmadillo.h> // Note the changed include and new attribute
// [[Rcpp::depends(RcppArmadillo)]]
// Avoid including header if openmp not on system
#ifdef _OPENMP
#include <omp.h>
#endif
// [[Rcpp::plugins(openmp)]]
// Note the changed return type
// [[Rcpp::export]]
arma::mat mandelRcpp(const double x_min, const double x_max,
const double y_min, const double y_max,
const int res_x, const int res_y, const int nb_iter)
arma::mat ret(res_x, res_y); // note change
double x_step = (x_max - x_min) / res_x;
double y_step = (y_max - y_min) / res_y;
unsigned r,c;
#pragma omp parallel for shared(res)
for (r = 0; r < res_y; r++)
for (c = 0; c < res_x; c++)
double zx = 0.0, zy = 0.0, new_zx;
double cx = x_min + c*x_step, cy = y_min + r*y_step;
unsigned n = 0;
for (; (zx*zx + zy*zy < 4.0 ) && ( n < nb_iter ); n++ )
new_zx = zx*zx - zy*zy + cx;
zy = 2.0*zx*zy + cy;
zx = new_zx;
if(n == nb_iter)
n = 0;
ret(r, c) = n;
return ret;
使用测试代码（注意 y 和 x 没有定义，因此我假设 y = ylims 和 x = xlims）我们有：
xlims = ylims = c(-2.0, 2.0)
x_res = y_res = 400L
nb_iter = 256L
system.time(m <-
mandelRcpp(xlims[[1]], xlims[[2]],
ylims[[1]], ylims[[2]],
x_res, y_res, nb_iter))
rainbow = c(
rgb(0.47, 0.11, 0.53),
rgb(0.27, 0.18, 0.73),
rgb(0.25, 0.39, 0.81),
rgb(0.30, 0.57, 0.75),
rgb(0.39, 0.67, 0.60),
rgb(0.51, 0.73, 0.44),
rgb(0.67, 0.74, 0.32),
rgb(0.81, 0.71, 0.26),
rgb(0.89, 0.60, 0.22),
rgb(0.89, 0.39, 0.18),
rgb(0.86, 0.13, 0.13)
)
cols = c(colorRampPalette(rainbow)(100),
rev(colorRampPalette(rainbow)(100)),
"black") # palette
par(mar = c(0, 0, 0, 0))
image(m,
col = cols,
asp = diff(range(ylims)) / diff(range(xlims)),
axes = F)
为：
【讨论】：
非常感谢！与此同时，我发现使用 ret(r,c) = n;而不是 ret[r,c] = n; （并添加 return ret; 我已经愚蠢地忘记了）确实会产生正确的结果 - 将测试哪个是最快的！您是否还知道我是否可以在 #pragma 中添加 simd 选项？如果我能从将更多变量定义为私有变量中受益？ 您可以在这些对象上使用新的 SIMD 结构。关于私有变量，嗯......这些变量是私有的。因此，您在每个线程的内存中为每个私有变量显式创建一个单独的副本。不确定是否会有所收获。 哈，是的，我现在明白了——谢谢！我玩了#pragma omp parallel for simd #pragma omp for simd #pragma omp simd 但这些似乎都没有帮助这里的性能...... @TomWenseleers 您需要手动矢量化。这种优化对于编译器来说太高级了。您必须保留在其他人之前完成的像素，并使用蒙版查找所有完成的时间，然后移动到下一个。 听起来很酷 - 如果它可以以一种相当优雅的方式完成，我很乐意看到它工作！以上是关于使用 Rcpp 和 OpenMP 在 R 中多线程和 SIMD 矢量化 Mandelbrot的主要内容，如果未能解决你的问题，请参考以下文章
在 Windows 上使用 C++11 和 Rcpp 构建 R 包