作为OpenMP
& Rcpp
性能测试我想检查使用最直接和简单的方法在 R 中计算 Mandelbrot 集的速度有多快 Rcpp
+ OpenMP
执行。目前我所做的是:
#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) collapse(2)
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 线程速度约为。比矢量化普通 R Mandelbrot 版本快 350 倍 here并且比(非多线程)Python/Numba 版本快约 7.3 倍 here (类似于 PyCUDA 或 PyOpenCL 速度),对此非常满意... Rasterizing/display now seems the bottleneck in R....
最佳答案
请勿不要将OpenMP与Rcpp的*Vector
一起使用或*Matrix
遮盖物体 SEXP
单线程的函数/内存分配。 OpenMP 是 multi-threaded approach 。
这就是代码崩溃的原因。
解决此限制的一种方法是使用非R数据结构来存储结果。以下之一就足够了:arma::mat
或Eigen::MatrixXd
或std::vector<T>
...由于我喜欢 Armadillo ,我将更改 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)
对于:
关于multithreading - 使用 Rcpp 和 OpenMP 在 R 中实现多线程和 SIMD 矢量化 Mandelbrot,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/48069990/