爲什麼我的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）。

祝你好運