Skip to content

Instantly share code, notes, and snippets.

@SimonDanisch SimonDanisch/benchmark.jl Secret
Last active Jul 26, 2017

Embed
What would you like to do?
#-------------------------------------
# Function for Jacobi sweep with loops
#-------------------------------------
function regular_time_step(u::Matrix)
n, m = size(u)
err = 0.0
for i in 2:n-1
for j in 2:m-1
tmp = u[i,j]
u[i,j] = ((u[i-1, j] + u[i+1, j] +
u[i, j-1] + u[i, j+1])*4.0 +
u[i-1,j-1] + u[i-1,j+1] +
u[i+1,j-1] + u[i+1,j+1])/20.0
diff = u[i,j] - tmp
err = err + diff*diff
end
end
return sqrt(err)
end
function optimized_time_step{T}(u::Matrix{T})
n, m = size(u)
err = T(0.0) # not an optimization, but makes it work for all element types
@inbounds for i = 2:n-1
for j = 2:m-1
tmp = u[j, i];
u[j, i] = ( T(4) *(u[j-1, i] + u[j+1, i] +
u[j, i-1] + u[j, i+1]) +
u[j-1, i-1] + u[j+1, i+1] +
u[j+1, i-1] + u[j-1, i+1]) / T(20)
diff = u[j, i] - tmp
err += diff*diff;
end
end
return sqrt(err)
end
function optimized_time_step_simd{T}(u::Matrix{T})
n, m = size(u)
err = T(0.0) # not an optimization, but makes it work for all element types
@inbounds for i = 2:n-1
@simd for j = 2:m-1
tmp = u[j, i];
u[j, i] = ( T(4) *(u[j-1, i] + u[j+1, i] +
u[j, i-1] + u[j, i+1]) +
u[j-1, i-1] + u[j+1, i+1] +
u[j+1, i-1] + u[j-1, i+1]) / T(20)
diff = u[j, i] - tmp
err += diff*diff;
end
end
return sqrt(err)
end
#-------------------------------------
# Function for vectorized Jacobi sweep
#-------------------------------------
function optimized_vectorized_time_step(u::Matrix)
(n, m) = size(u)
u_old = copy(u)
@views u[2:n-1, 2:m-1] = ((u[1:n-2, 2:m-1] + u[3:n, 2:m-1] +
u[2:n-1, 1:m-2] + u[2:n-1, 3:m])*4.0 +
u[1:n-2, 1:m-2] + u[1:n-2, 3:m] +
u[3:n, 1:m-2] + u[3:n, 3:m])/20.0
v = vec(u - u_old)
return sqrt(dot(v, v))
end
# C version
const clib_name5 = tempname()
source = """
#include <math.h>
double c_time_step(int N, double *u[])
{
double tmp, diff;
int i, j;
double err;
err = 0.0;
for (i = 1; i < N-1; i++) {
for (j = 1; j < N-1; j++) {
tmp = u[j][i];
u[j][i] = ( 4.0 *(u[j-1][i] + u[j+1][i] +
u[j][i-1] + u[j][i+1]) +
u[j-1][i-1] + u[j+1][i+1] +
u[j+1][i-1] + u[j-1][i+1])/20.;
diff = u[j][i] - tmp;
err += diff*diff;
}
}
return sqrt(err);
}
"""
run(`gcc --version`)
# gcc (Ubuntu 5.4.1-2ubuntu1~16.04) 5.4.1 20160904
open(`gcc -fPIC -Ofast -O3 -xc -shared -o $(clib_name5).so -`, "w") do pipe
print(pipe, source)
end
function c_time_step(X::Vector{Vector{Float64}})
ccall(("c_time_step", clib_name5), Float64, (Cint, Ptr{Ptr{Float64}}), length(X), X)
end
function solver(u, time_step)
iter = 0
err = 2.0
while iter < 100000 && err > 1e-6
err = time_step(u)
iter += 1
end
u, err, iter
end
function init_array(numPoints = 200)
pi_c = 4.0 * atan(1.0)
v = zeros(Float64, numPoints, numPoints)
x = [(i-1)*pi_c/(numPoints-1) for i in 1:numPoints]
v[1,:] = [sin(a) for a in x]
v[numPoints,:] = [sin(a)*exp(-pi_c) for a in x]
v
end
init_c() = (v = init_array(); v_2 = [v[i, :] for i=1:size(v, 1)])
funcs = (regular_time_step, optimized_time_step, optimized_vectorized_time_step, optimized_time_step_simd) #
benchmarks = map(funcs) do time_step
v = init_array()
solver(v, time_step) # compile
v = init_array()
# should be enough iterations already to be okay with just a simple elapsed
result = @elapsed begin
u, err, it = solver(v, time_step)
end
Symbol(time_step) => (result, it)
end
v = init_c()
result = @elapsed begin
u, err, it = solver(v, c_time_step)
end
benchmarks = (benchmarks..., :c_time_step => (result, it))
for (name, (t, sweeps)) in benchmarks
println(name, ":")
println(" Number of sweeps: ", sweeps)
println(" Time Jacobi Iteration: $(t) Seconds")
end
#= with julia -O3 --math-mode=fast, which seems a fairer comparison to gcc -Ofast -O3
regular_time_step:
Number of sweeps: 30712
Time Jacobi Iteration: 10.331530467 Seconds
optimized_time_step:
Number of sweeps: 30727
Time Jacobi Iteration: 10.150175164 Seconds
optimized_vectorized_time_step:
Number of sweeps: 56877
Time Jacobi Iteration: 20.465307448 Seconds
optimized_time_step_simd:
Number of sweeps: 38836
Time Jacobi Iteration: 4.759529088 Seconds
c_time_step:
Number of sweeps: 30727
Time Jacobi Iteration: 15.875432439 Seconds
=#
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.