Dear Tim,
This will definitely fix it - thank you. I did not realise it could be done
this way.
As usual - thank you for all your help.
Christoph
On Sunday, 4 January 2015 01:56:40 UTC, Tim Holy wrote:
>
> 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
>
>
