为什么我的Julia实现用于计算3D中的欧几里德距离比我的C实现更快

Question

我正在比较Julia计算3D空间中两组点之间的欧几里德距离与C中的等效实现之间的距离。惊讶地看到，（对于这种特殊情况下，我的具体实现）Julia是比C快22％。当我还包括在朱莉娅版本@fastmath，这将是比C

这导致更快的甚至83％对我的问题：为什么？无论是朱莉娅比我原本想象的更加神奇或我在C做的事情效率非常低。我打赌我的钱在后者上。

有关实施细节的一些：

• 在朱莉娅我用的Float64二维数组。
• 在C中，我使用动态分配的double的一维数组。
• 在C中，我使用math.h中的sqrt函数。
• 计算速度非常快，因此我计算它们1000倍以避免在微/毫秒级别进行比较。
  • 编译器：GCC 5.4.0
  • 优化标志：-O3 -ffast-math

  时序：

  • 朱莉娅（不@fastmath）

  关于编译的一些细节：90 s

• 朱（与@fastmath）：20秒
• C：116章第
• 我使用bash命令time用于定时
  • $ time ./particleDistance.jl（与认领在文件）
  • $ time ./particleDistance

particleDistance.j升

#!/usr/local/bin/julia 

function distance!(x::Array{Float64, 2}, y::Array{Float64, 2}, r::Array{Float64, 2}) 
    nx = size(x, 1) 
    ny = size(y, 1) 

    for k = 1:1000 

     for j = 1:ny 

      @fastmath for i = 1:nx 
       @inbounds dx = y[j, 1] - x[i, 1] 
       @inbounds dy = y[j, 2] - x[i, 2] 
       @inbounds dz = y[j, 3] - x[i, 3] 

       rSq = dx*dx + dy*dy + dz*dz 

       @inbounds r[i, j] = sqrt(rSq) 
      end 

     end 

    end 

end 

function main() 
    n = 4096 
    m = 4096 

    x = rand(n, 3) 
    y = rand(m, 3) 
    r = zeros(n, m) 

    distance!(x, y, r) 

    println("r[n, m] = $(r[n, m])") 
end 

main() 

particleDistance.c

#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 

void distance(int n, int m, double* x, double* y, double* r) 
{ 
    int i, j, I, J; 
    double dx, dy, dz, rSq; 

    for (int k = 0; k < 1000; k++) 
    { 
     for (j = 0; j < m; j++) 
     { 
      J = 3*j; 

      for (i = 0; i < n; i++) 
      { 
       I = 3*i; 

       dx = y[J] - x[I]; 
       dy = y[J+1] - x[I+1]; 
       dz = y[J+2] - x[I+2]; 

       rSq = dx*dx + dy*dy + dz*dz; 

       r[j*n+i] = sqrt(rSq); 
      } 
     } 
    } 
} 

int main() 
{ 
    int i; 
    int n = 4096; 
    int m = 4096; 

    double *x, *y, *r; 

    size_t xbytes = 3*n*sizeof(double); 
    size_t ybytes = 3*m*sizeof(double); 

    x = (double*) malloc(xbytes); 
    y = (double*) malloc(ybytes); 
    r = (double*) malloc(xbytes*ybytes/9); 

    for (i = 0; i < 3*n; i++) 
    { 
     x[i] = (double) rand()/RAND_MAX*2.0-1.0; 
    } 

    for (i = 0; i < 3*m; i++) 
    { 
     y[i] = (double) rand()/RAND_MAX*2.0-1.0; 
    } 

    distance(n, m, x, y, r); 

    printf("r[n*m-1] = %f\n", r[n*m-1]); 

    free(x); 
    free(y); 
    free(r); 

    return 0; 
} 

的Makefile

all: particleDistance.c 
    gcc -o particleDistance particleDistance.c -O3 -ffast-math -lm 

Answer 1

也许这应该是一个评论，但问题是，朱莉娅确实是相当优化。在Julia网页中，您可以看到它在某些情况下可以击败C（mandel）。

我发现你在汇编中使用了-ffast-math。但是，也许你可以在你的代码中做一些优化（虽然现在的编译器非常聪明，这可能无法解决问题）。

1. 而不是int使用您的索引，尝试使用unsigned int，这可以让你尝试以下事情;
2. 而不是乘以3，如果你使用unsigned，你可以做一个移位并添加。这可以节省一些计算时间;
3. 在访问像x [J]这样的元素时，可能会尝试直接使用指针并以顺序方式访问元素，如x + = 3（？）;
4. 而不是int n和int m，尝试将它们设置为宏。如果他们事先知道，你可以利用这一点。
5. malloc在这种情况下是否有所不同？如果已知n和m，则固定大小的数组将减少OS分配内存所花费的时间。

可能有一些其他的东西，但朱莉娅是相当具有实时编译优化，让一切是恒定的，并事先知道的是有利于它的使用。我已经试过朱莉娅，没有遗憾。

Answer 2

在C该指数的计算是相当缓慢的

尝试类似如下（我没编译它，它可能仍然错误，太直观的想法）：

void distance(int n, int m, double* x, double* y, double* r) 
{ 
int i, j; 
double dx, dy, dz, rSq; 
double* X, *Y, *R; 


for (int k = 0; k < 1000; k++) 
{ 
    R = r; 
    Y = y; 
    for (j = 0; j < m; j++) 
    { 
     X = x; 

     for (i = 0; i < n; i++) 
     { 
      dx = Y[0] - *X++; 
      dy = Y[1] - *X++; 
      dz = Y[2] - *X++; 

      rSq = dx*dx + dy*dy + dz*dz; 

      *R++ = sqrt(rSq); 
     } 
     Y += 3; 
    } 
} 
} 

另外，您可以尝试，这可能是一个有点快（一个增量而不是3）

  dx = Y[0] - X[0]; 
      dy = Y[1] - X[1]; 
      dz = Y[2] - X[2]; 
      X+=3; 

Y [X]相同*（Y + X）。

祝你好运