# Code for last post. Christoph
type SiteLinear
stencil::Array{Int, 1}
coeffs::Array{Float64, 3}
end
function evalSiteFunction_without_type!(coeffs::Array{Float64, 3},
stencil::Array{Int, 1},
J::Array{Int, 1},
Y::Array{Float64,2},
Frc::Array{Float64,2})
# extract dimension information (d1, d2 \in \{2, 3\})
d1, d2, nneig = size(coeffs)
# allocate some variables
kX = 0::Int; jX = 0::Int; n = 0::Int
# loop over sites
for jX in J # length(J) = 11268
# initialise force to 0
Frc[:,jX] = 0.0
# loop over neighbouring sites
for n in 1:nneig # nneig = 37
# get the X-index of the current neighbour (highly simplified)
kX = jX + stencil[n]
# evaluate the term (devectorized and unrolled for performance)
# f[:] += slice(sp.coeffs, :, :, n) * Y[:, kX]
for i = 1:d1, j = 1:d2
Frc[i,jX] += coeffs[i,j,n]*Y[j,kX]
end
end
end
end
function evalSiteFunction!(sp::SiteLinear,
J::Array{Int, 1},
Y::Array{Float64,2},
Frc::Array{Float64,2})
# extract dimension information (d1, d2 \in \{2, 3\})
d1, d2, nneig = size(coeffs)
# allocate some variables
kX = 0::Int; jX = 0::Int; n = 0::Int
# loop over sites
for jX in J # length(J) = 11268
# initialise force to 0
Frc[:,jX] = 0.0
# loop over neighbouring sites
for n in 1:nneig # nneig = 37
# get the X-index of the current neighbour (highly simplified)
kX = jX + sp.stencil[n]
# evaluate the term (devectorized and unrolled for performance)
# f[:] += slice(sp.coeffs, :, :, n) * Y[:, kX]
for i = 1:d1, j = 1:d2
Frc[i,jX] += sp.coeffs[i,j,n]*Y[j,kX]
end
end
end
end
# problem size
nneigs = 18
nsites = 11268
dim = 2
# data to be passed to function (simplified)
coeffs = rand(dim,dim,2*nneigs+1)
stencil = int(linspace(-nneigs, nneigs, 2*nneigs+1))
Y = rand(2, nsites + 2*nneigs)
Frc = zeros(2, nsites + 2*nneigs)
J = int(linspace(nneigs+1, nsites+nneigs, nsites))
evalSiteFunction_without_type!(coeffs, stencil, J, Y, Frc)
@time evalSiteFunction_without_type!(coeffs, stencil, J, Y, Frc)
# now try same with
sp = SiteLinear(stencil, coeffs)
evalSiteFunction!(sp, J, Y, Frc)
@time evalSiteFunction!(sp, J, Y, Frc)
