Dear All,
Thank you so much for the various hints and tips. After re-reading the
performance Tips I've now found all the problems with the code: in the end
there were two type-instabilities, one that was easy to resolve (declare an
array to be Int instead of Integer) but the other is still a problem for me
so I'd love to get more feedback on this. (Initially I was looking for the
wrong thing as I hadn't realised that type-instability can cause unneeded
memory allocation.)
So my remaining issue is this. To get the code to run efficiently I had to
tell the compile that a certain array is 2-dimensional.
I = geom.I::Array{Int, 2}
In the type for geom, it is only declared as I::Array{Int}, because it
could in fact be 1-, 2- or 3-dimensional. At the moment, I see two ways to
fix this:
* write three functions, and a wrapper which checks for the correct
dimension. (Or possibly wrap my head around meta-programming and do it this
way. This would probably have the advantage that I could also unroll the
inner loops and get some additional factors)
* replace I with a one-dimensional array and resolve the multi-dimensional
indexing manually, probably some slow-down because of checking the dimension
BUT, is there a quicker/easier/more readable fix?
Many thanks,
Christoph
function evalSiteFunction!(sp::SiteLinear, geom::tbgeom,
Y::Array{Float64,2},
Frc::Array{Float64,2})
# extract dimension information (d1, d2 \in \{2, 3\})
I = geom.I::Array{Int, 2}
d1, d2, nneig = size(sp.coeffs)
# assert type of the index set over which we are looping
for jX in geom["iMM"]::Array{Int, 1}
# initialise force to 0
Frc[:,jX] = 0.0
# loop over neighbouring sites
@inbounds for n in 1:nneig
# get the X-index of the current neighbour
kX = getI(I, geom.X2I, jX, sp.stencil, n)
# evaluate the term (devectorized for performance)
for i = 1:d1
Frc[i, jX] = sp.coeffs[i, 1, n] * Y[1, kX]
for j = 2:d2
Frc[i,jX] += sp.coeffs[i, j, n] * Y[j, kX]
end
end
end
end
end
