朱莉娅将性能与R进行比较的例子似乎特别令人费解. https://github.com/JuliaLang/julia/blob/master/test/perf/perf.R
你可以从下面的两种算法中获得最快的性能(最好解释你改变了什么以使它更像R)？
## mandel
mandel = function(z) {
c = z
maxiter = 80
for (n in 1:maxiter) {
if (Mod(z) > 2) return(n-1)
z = z^2+c
}
return(maxiter)
}
mandelperf = function() {
re = seq(-2,0.5,.1)
im = seq(-1,1,.1)
M = matrix(0.0,nrow=length(re),ncol=length(im))
count = 1
for (r in re) {
for (i in im) {
M[count] = mandel(complex(real=r,imag=i))
count = count + 1
}
}
return(M)
}
assert(sum(mandelperf()) == 14791)
## quicksort ##
qsort_kernel = function(a, …对于我在R&Rcpp + OpenMP&Shiny中制作的实时交互式Mandelbrot查看器我正在寻找一种高性能的方式来显示1920x1080矩阵作为光栅图像,以期能够实现ca. 5-10 fps(计算Mandelbrot图像本身现在在中等变焦时达到约20-30 fps,当然滚动应该快速).使用image()with选项useRaster=TRUE,plot.raster或者甚至grid.raster()还没有完全削减它,所以我正在寻找更高性能的选项,理想情况下使用OpenGL加速.
我注意到有qt 包装器qtutils和qtpaint
http://finzi.psych.upenn.edu/R/library/qtutils/html/sceneDevice.html
你可以在其中设置参数opengl=TRUE和
http://finzi.psych.upenn.edu/ R/library/qtpaint/html/qplotView.html
再次使用参数opengl=TRUE和http://finzi.psych.upenn.edu/R/library/qtpaint/html/painting.html.
我还注意到应该能够使用rdyncall软件包调用SDL和GL/OpenGL函数(从https://cran.r-project.org/src/contrib/Archive/rdyncall/安装,从https://安装SDL)www.libsdl.org/download-1.2.php)`,演示可在http://hg.dyncall.org/pub/dyncall/bindings/file/87fd9f34eaa0/R/rdyncall/demo/00Index获得,例如http:// hg.dyncall.org/pub/dyncall/bindings/file/87fd9f34eaa0/R/rdyncall/demo/randomfield.R).
我是否正确使用这些软件包,应该能够使用opengl加速显示2D图像栅格？如果是这样,有没有人想过如何做到这一点(我问,因为我不是qt或SDL/OpenGL的专家)？
非OpenGL选项的某些时间对我的应用来说太慢了:
# some example data & desired colour mapping of [0-1] ranged data matrix
library(RColorBrewer)
ncol=1080
cols=colorRampPalette(RColorBrewer::brewer.pal(11, "RdYlBu"))(ncol)
colfun=colorRamp(RColorBrewer::brewer.pal(11, "RdYlBu"))
col = rgb(colfun(seq(0,1, length.out …在R中给出这样的数据帧:
+---+---+
| X | Y |
+---+---+
| 1 | 2 |
| 2 | 4 |
| 4 | 5 |
+---+---+
如果对此数据帧执行矢量化操作,如下所示:
data$Z <- data$X * data$Y
这会利用处理器的单指令多数据(SIMD)功能来优化性能吗？这似乎是一个完美的案例,但我找不到任何证实我的预感的东西.
当我第一次使用Haswell处理器时,我尝试使用FMA来确定Mandelbrot集.主要算法是这样的:
intn = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn x2 = square(x), y2 = square(y); //square(x) = x*x
floatn r2 = x2 + y2;
booln mask = r2<cut; //booln is in the float domain non integer domain
if(!horizontal_or(mask)) break; //_mm256_testz_pd(mask)
n -= mask
floatn t = x*y; mul2(t); //mul2(t): t*=2
x = x2 - y2 + cx;
y = t + cy;
}
这确定n像素是否在Mandelbrot集中.因此对于双浮点,它运行超过4个像素(floatn = __m256d,intn = __m256i).这需要4个SIMD浮点乘法和4个SIMD浮点加法.
然后我修改了这个就像这样使用FMA
intn n = 0; …我正在寻找在AVX元件(单精度浮点)上运行的指数函数的有效(快速)近似.即 - __m256 _mm256_exp_ps( __m256 x )没有SVML.
相对精度应该类似于~1e-6,或~20个尾数位(1 ^ 2 ^ 20).
如果用英特尔内在函数用C风格编写,我会很高兴.
代码应该是可移植的(Windows,macOS,Linux,MSVC,ICC,GCC等).
这类似于使用SSE的指数函数的最快实现,但是这个问题寻求非常快速且精度低(当前的答案提供了大约1e-3的精度).
此外,这个问题是寻找AVX/AVX2(和FMA).但请注意,这两个问题的答案很容易在SSE4 __m128或AVX2 之间移植__m256,因此未来读者应根据所需的精度/性能权衡进行选择.
我如何编写一个可移植的GNU C内置向量版本,它不依赖于x86 set1内在函数？
typedef uint16_t v8su __attribute__((vector_size(16)));
v8su set1_u16_x86(uint16_t scalar) {
return (v8su)_mm_set1_epi16(scalar); // cast needed for gcc
}
当然必须有一个更好的方式
v8su set1_u16(uint16_t s) {
return (v8su){s,s,s,s, s,s,s,s};
}
我不想写一个用于广播单个字节的AVX2版本!
对于你想要分配给变量而不是仅仅使用二元运算符的操作数(这与gcc一起工作,见下文),即使只对gcc-only或clang-only这一部分的答案也会很有趣.
如果我想使用广播标量作为二元运算符的一个操作数,则可以使用gcc(如手册中所述),但不能使用clang:
v8su vecdiv10(v8su v) { return v / 10; } // doesn't compile with clang
有了clang,如果我只针对x86并且只使用本机向量语法来让编译器为我生成模数乘法逆常数和指令,我可以写:
v8su vecdiv_set1(v8su v) {
return v / (v8su)_mm_set1_epi16(10); // gcc needs the cast
}
但是如果我加宽向量(to _mm256_set1_epi16),我必须改变内在,而不是通过改变到vector_size(32)一个地方(对于不需要改组的纯垂直SIMD )将整个代码转换为AVX2 .它也会破坏本机向量的部分目的,因为它不会为ARM或任何非x86目标编译.
丑陋的演员阵容是必需的,因为与clang不同,gcc不v8us {aka __vector(8) short …
我看到如下代码:
#include "stdio.h"
#define VECTOR_SIZE 4
typedef float v4sf __attribute__ ((vector_size(sizeof(float)*VECTOR_SIZE)));
// vector of four single floats
typedef union f4vector
{
v4sf v;
float f[VECTOR_SIZE];
} f4vector;
void print_vector (f4vector *v)
{
printf("%f,%f,%f,%f\n", v->f[0], v->f[1], v->f[2], v->f[3]);
}
int main()
{
union f4vector a, b, c;
a.v = (v4sf){1.2, 2.3, 3.4, 4.5};
b.v = (v4sf){5., 6., 7., 8.};
c.v = a.v + b.v;
print_vector(&a);
print_vector(&b);
print_vector(&c);
}
这段代码构建良好,并且预期使用gcc(它是内置的SSE/MMX扩展和向量数据类型.这段代码使用4个单浮点进行SIMD向量加法.
我想详细了解每个关键字/函数调用此typedef行的含义和含义:
typedef float v4sf __attribute__ ((vector_size(sizeof(float)*VECTOR_SIZE)));
什么是vector_size()函数返回;
__attribute__关键字是什么？
这是浮点数据类型是为vfsf类型定义的类型？
我理解其余部分.
谢谢, …
有人可以解释一下本文中描述的扰动是如何加速渲染Mandelbrot集的吗？
我知道如何使用传统方法渲染Mandelbrot集合,其中对每个像素执行许多迭代,但我不太明白该文章中描述的是什么.
我像这样计算参考轨道:
std::complex<double> Xo(some_x, some_y);
std::complex<double> Xn(0,0);
for (int n = 0; n < maxIterations; ++n) {
orbit.push_back(Xn);
Xn = Xn * Xn + Xo;
}
那是对的吗？那么我如何使用参考轨道计算所有其他像素？
有人可以帮助我改进这些代码并给我一些提示.我试图自己创建一个OpenMP版本的Mandelbrot.我是OpenMP初学者,在这里我没有加快速度,这可能是因为#pragma omp critical我现在想不出更好的主意.
int main()
{
// picture resolution
int iX,iY;
const int ImageWidth = 1000;
const int ImageHeight = 1000;
double Cx,Cy;
const double CxMin=-2.5;
const double CxMax=1.5;
const double CyMin=-2.0;
const double CyMax=2.0;
double PixelWidth=(CxMax-CxMin)/ImageWidth; /* scaled x coordinate of pixel (must be scaled to lie somewhere in the Mandelbrot
X scale (-2.5, 1.5) */
double PixelHeight=(CyMax-CyMin)/ImageHeight;/* scaled y coordinate of pixel (must be scaled to lie somewhere in the Mandelbrot
Y scale (-2.0, 2.0) */ …我有一个函数计算mandelbrot集我试图使用它并行化openMP.
我#pragma omp parallel for private在每个之前for
static void calculer (Image * im, int nb_iter, double x_min, double x_max, double y_min, double y_max) {
/* Discretisation de l'ensemble */
double pasx = (x_max - x_min) / im -> nb_col;
double pasy = (y_max - y_min) / im -> nb_lig;
double cy = y_min;
double new_zx;
unsigned int l,c;
// Calcul
#pragma omp parallel for private ( pasx, pasy, im,nb_iter,x_min,x_max,y_min, y_max)
for (l = 0; …