Why is my Julia implementation for computing Euclidean distances in 3D faster than my C implementation

Question

I am comparing the time it takes Julia to compute the Euclidean distances between two sets of points in 3D space against an equivalent implementation in C. I was very surprised to observe that (for this particular case and my particular implementations) Julia is 22% faster than C. When I also included @fastmath in the Julia version, it would be even 83% faster than C.

This leads to my question: why? Either Julia is more amazing than I originally thought or I am doing something very inefficient in C. I am betting my money on the latter.

Some particulars about the implementation:

• In Julia I use 2D arrays of Float64.
• In C I use dynamically allocated 1D arrays of double.
• In C I use the sqrt function from math.h.
• The computations are very fast, therefore I compute them a 1000 times to avoid comparing on the micro/millisecond level.

Some particulars about the compilation:

• Compiler: gcc 5.4.0
• Optimisation flags: -O3 -ffast-math

Timings:

• Julia (without @fastmath): 90 s
• Julia (with @fastmath): 20 s
• C: 116 s
• I use the bash command time for the timings
  • $ time ./particleDistance.jl (with shebang in file)
  • $ time ./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;
}

Makefile

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

Answer 1

Maybe it should be a comment, but the point is that Julia is indeed pretty optimized. In the Julia web page you can see that it can beat C in some cases (mandel).

I see that you are using -ffast-math in your compilation. But, maybe you could do some optimizations in your code (although nowadays compilers are pretty smart and this might not solve the issue).

1. Instead of using int for your indexes, try to use unsigned int, this allows you to maybe try the following thing;
2. Instead of multiply by 3, if you use an unsigned you can do a shift and add. This can save some computation time;
3. In accessing the elements like x[J], maybe try using pointers directly and access the elements in a sequential manner like x+=3 (?);
4. Instead of int n and int m, try to set them as macros. If they are known in advance, you can take advantage of that.
5. Does the malloc make difference in this case? If n and m are known, fixed size arrays would reduce the time spent for the OS allocate memory.

There might be a few other things, but Julia is pretty optimized with real time compilation, so everything that is constant and is known in advance is used in favor of it. I have tried Julia with no regrets.

Answer 2

Your index calculation in C is rather slow

Try something like the following (I did not compiled it, it may have still errors, just too visualize the idea):

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;
    }
}
}

Alternatively you could try, it might be a little bit faster (one increment instead of 3)

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

Y[x] is the same as *(Y+x).

Good luck