Laplacian operator in Julia and C. Version 2.

laplace.c

#include "laplace.h"

void laplace3D (double *del_field, double *field, long n1, long n2, long n3) {

	long i1, i2, i3, idx1, idx2, idx3, idx4, idx5, idx6, idx7;

	for (i3 = 1; i3 < n3-1; i3++) {
		for (i2 = 1; i2 < n2-1; i2++) {
			for (i1 = 1; i1 < n1-1; i1++) {
				idx1 = i3 * n2 * n1 + i2 * n1 + i1;
				idx2 = i3 * n2 * n1 + i2 * n1 + (i1 + 1);
				idx3 = i3 * n2 * n1 + i2 * n1 + (i1 - 1);
				idx4 = i3 * n2 * n1 + (i2 + 1) * n1 + i1;
				idx5 = i3 * n2 * n1 + (i2 - 1) * n1 + i1;
				idx6 = (i3 + 1) * n2 * n1 + i2 * n1 + i1;
				idx7 = (i3 - 1) * n2 * n1 + i2 * n1 + i1;
				
				del_field[idx1] = 
					field[idx2] + field[idx3] +
					field[idx4] + field[idx5] +
					field[idx6] + field[idx7] -
					6.0 * field[idx1];
			}
		}

	}

}

laplace.h

#ifndef __LAPLACE__
#define __LAPLACE__

void laplace3D (double *del_field, double *field, long n1, long n2, long n3);

#endif

laplace.jl

function laplace_idiomatic (del_x::Array{Float64, 3}, x::Array{Float64, 3})

        n1, n2, n3 = size (x)
        for i3 = 2:n3-1, i2 = 2:n2-1, i1 = 2:n1-1

                del_x[i1,i2,i3] =
                        x[i1+1,i2,i3] + x[i1-1,i2,i3] +
                        x[i1,i2+1,i3] + x[i1,i2-1,i3] +
                        x[i1,i2,i3+1] + x[i1,i2,i3-1] -
                        6.0 * x[i1,i2,i3]
        end

end

function laplace_jameson (del_x::Array{Float64, 3}, x::Array{Float64, 3})
                
        n1, n2, n3 = size (x)
        for i3 = 2:n3-1, i2 = 2:n2-1, i1 = 2:n1-1
                a =
                        x[i1+1,i2,i3] + x[i1-1,i2,i3] +
                        x[i1,i2+1,i3] + x[i1,i2-1,i3]
                b =
                        x[i1,i2,i3+1] + x[i1,i2,i3-1] -
                        6.0 * x[i1,i2,i3]
                del_x[i1,i2,i3] = a + b
        end

end

function laplace_ariel (del_x::Array{Float64, 3}, x::Array{Float64, 3})

        n1, n2, n3 = size (x)
        n12 = n1 * n2
        for i3 = 2:n3-1, i2 = 2:n2-1, i1 = 2:n1-1
                index = i1 + n1 * ((i2 - 1) + n2 * (i3 - 1))
                tmp1 = -6.0 * x[index] + x[index - 1] + x[index + 1]
                tmp2 = x[index - n1] + x[index + n1] + x[index - n12]
                del_x[index] = x[index + n12] + tmp1 + tmp2
        end

end

function laplace_gunnar_no_bounds (del_x::Array{Float64, 3}, x::Array{Float64, 3})

    n1, n2, n3 = size (x)
    n12 = n1 * n2
    xp = pointer(x)
    del_xp = pointer(del_x)
    for i3 = 2:n3-1, i2 = 2:n2-1, i1 = 2:n1-1
        index = i1 + n1 * ((i2 - 1) + n2 * (i3 - 1))
        t = -6.0 * unsafe_load(xp, index)
        t += unsafe_load(xp, index - 1)
        t += unsafe_load(xp, index + 1)
        t += unsafe_load(xp, index - n1)
        t += unsafe_load(xp, index + n1)
        t += unsafe_load(xp, index - n12)
        t += unsafe_load(xp, index + n12)
        unsafe_store!(del_xp, t, index)
    end

end

function laplace_c (del_x::Array{Float64, 3}, x::Array{Float64, 3})

        n1, n2, n3 = size (x)

        ccall ( (:laplace3D, "liblaplace"),
                Void,
                (Ptr{Float64}, Ptr{Float64}, Int64, Int64, Int64),
                del_x, x, n1, n2, n3 )

end

laplace_bench.jl

require ("laplace.jl")

function main ()
        niters = 10
        n = 300
        field = init_field (n)
        del_field = zeros (size (field))

        # warmup Julia
        laplace_idiomatic (del_field, field)
        laplace_jameson (del_field, field)
        laplace_ariel (del_field, field)
        laplace_gunnar_no_bounds (del_field, field)

        # timing idiomatic Julia version
        tic ()
        for n = 1:niters
                laplace_idiomatic (del_field, field)
        end
        write (STDOUT, "Julia (idomatic): ")
        toc ()

        # timing Jameson's Julia version
        tic ()
        for n = 1:niters
                laplace_jameson (del_field, field)
        end
        write (STDOUT, "Julia (Jameson): ")
        toc ()
        
        # timing Ariel's Julia version
        tic ()
        for n = 1:niters
                laplace_ariel (del_field, field)
        end
        write (STDOUT, "Julia (Ariel): ")
        toc ()
        
        # timing Gunnar's (no bounds checking) Julia version
        tic ()
        for n = 1:niters
                laplace_gunnar_no_bounds (del_field, field)
        end
        write (STDOUT, "Julia (Gunnar, no bounds): ")
        toc ()

        # timing C (serial) version
        tic ()
        for n = 1:niters
                laplace_c (del_field, field)
        end
        write (STDOUT, "C (serial): ")
        toc ()

end

function init_field (n::Int64)

        field = zeros (n, n, n)
        for i1 = 1:n, i2 = 1:n, i3 = 1:n
                if i1 == 1 || i2 == 1 || i3 == 1 || i1 == n || i2 == n || i3 == n
                        field[i1,i2,i3] = 1.0
                end
        end
        field

end

main ()
        

Makefile

all : liblaplace.so

liblaplace.so : laplace.h laplace.c
	gcc -O3 -fPIC -c laplace.c
	gcc -O3 -fPIC -shared -o liblaplace.so laplace.o
	rm -f laplace.o

clean:
	rm -f liblaplace.so