使用Rcpp和OpenMP在R中进行多线程和SIMD矢量化Mandelbrot

Question

作为OpenMP和Rcpp性能测试，我想检查我使用最直接和简单的Rcpp + OpenMP实现计算R中Mandelbrot集的速度有多快。目前我所做的是：

#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 #pragma中使用simd选项似乎没有做任何事情）

PS起初我的代码崩溃但我后来发现这是通过用ret[r,c] = n;替换ret(r,c) = n;来解决的。使用Armadillo类，如下面的答案所示，使事情稍微快一点，虽然时间几乎相同。同时在x和y周围翻转，以便在用image()绘制时以正确的方向出现。使用8个线程的速度是ca.比矢量化普通R Mandelbrot版本here快350倍，并且比（非多线程）Python / Numba版本here（类似于PyCUDA或PyOpenCL速度）快7.3倍，所以非常满意... Rasterizing/display now seems the bottleneck in R....

Answer 1

不要将OpenMP与Rcpp的*Vector或*Matrix对象一起使用，因为它们会掩盖单线程的SEXP函数/内存分配。 OpenMP是一个multi-threaded approach。

这就是代码崩溃的原因。

解决此限制的一种方法是使用非R数据结构来存储结果。以下之一就足够了：arma::mat或Eigen::MatrixXd或std::vector<T> ......我赞成犰狳，我会将res矩阵从arma::mat改为Rcpp::NumericMatrix。因此，以下将并行执行您的代码：

#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)

对于：

Answer 2

我继续使用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实际上使用efficient code指令生成相对ptest（但not for AVX512），但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。

---

代码产生或多或少相同的图像。