c - 为什么我在 3D 中计算欧几里得距离的 Julia 实现比我的 C 实现更快

Question

我正在比较 Julia 计算 3D 空间中两组点之间的欧几里德距离与 C 中的等效实现所花费的时间。我很惊讶地观察到(对于这个特殊情况和我的特殊实现) Julia 比 C 快 22%。当我还在 Julia 版本中包含 @fastmath 时，它甚至比 C 快 83%。

这引出了我的问题:为什么？要么 Julia 比我原先想象的更神奇，要么 我正在用 C 做一些非常低效的事情。我把钱押在后者身上。

关于实现的一些细节:

• 在 Julia 中，我使用 Float64 的二维数组。
• 在 C 语言中，我使用动态分配的 double 一维数组。
• 在 C 中，我使用 math.h 中的 sqrt 函数。
• 计算速度非常快，因此我将它们计算 1000 次以避免在微/毫秒级别进行比较。

关于编译的一些细节:

• 编译器:gcc 5.4.0
• 优化标志:-O3 -ffast-math

时间:

• Julia(没有 @fastmath):90 秒
• Julia(使用 @fastmath):20 秒
• C: 116 秒
• 我使用 bash 命令 time 来计时
  • $ time ./particleDistance.jl(文件中有 shebang)
  • $时间./particleDistance

particleDistance.jl

#!/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;
}

生成文件

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

Answer 1

也许它应该是一个评论，但重点是 Julia 确实经过了优化。在 Julia 网页中，您可以看到它在某些情况下可以击败 C (mandel)。

我看到您在编译中使用了 -ffast-math。但是，也许您可​​以对代码进行一些优化(尽管现在的编译器非常智能，这可能无法解决问题)。

1. 尝试使用 unsigned int 而不是对索引使用 int，这允许您尝试以下操作；
2. 如果您使用无符号数，则可以进行移位和加法，而不是乘以 3。这可以节省一些计算时间；
3. 在访问像x[J]这样的元素时，可以尝试直接使用指针，像x+=3(?)这样的顺序访问元素；
4. 尝试将它们设置为宏，而不是 int n 和 int m。如果他们事先已知，您可以利用这一点。
5. malloc 在这种情况下有什么不同吗？如果已知 n 和 m，固定大小的数组将减少操作系统分配内存所花费的时间。

可能还有一些其他的东西，但 Julia 对实时编译进行了相当优化，所以所有常量和预先已知的东西都会被用来支持它。我毫无遗憾地尝试过 Julia。