The most julian way of doing this is to use the dimensionality (often called
N) as a parameter, either for the type or functions of the type:
type MyContainer{T,N}
data::Array{T,N}
end
There are a ton of examples of this kind of trick in base; reading through
those files is a great resource.
--Tim
On Saturday, January 03, 2015 02:19:34 PM Christoph Ortner wrote:
> 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
