【发布时间】：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/… 的两种方法...
标签： multithreading openmp rcpp simd mandelbrot