It seems that you are doing all the right things in the inner loops. Can
you post a self-sufficient code snippet that one can run?
-viral
On Saturday, January 3, 2015 6:02:45 PM UTC+5:30, Christoph Ortner wrote:
>
> Hi Tim,
> Many thanks for the suggestion, unfortunately it did not fix the issue;
> see below updated code + timing/allocation.
> All the best, Christoph
>
>
>
> function evalSiteFunction!(sp::SiteLinear, geom::tbgeom,
> Y::Array{Float64,2}, Frc::Array{Float64,2})
> # extract dimension information (d1, d2 \in \{2, 3\})
> d1, d2, nneig = size(sp.coeffs)
> # allocate some variables
> kX = 0::Int; jX = 0::Int
> # assert type of the index set over which we are looping
> J = geom["iMM"]::Array{Int, 1}
> for jX in J
> # initialise force to 0
> Frc[:,jX] = 0.0
> # loop over neighbouring sites
> for n in 1:nneig
> # get the X-index of the current neighbour
> kX = geom.I[ geom.X2I[1,jX] + sp.stencil[1,n], geom.X2I[2,jX]
> + sp.stencil[2,n] ]
> # evaluate the term (devectorized 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
>
> @time TB.evalSiteFunction!(frcMM, geom, Y, Frc)
>
> elapsed time: 0.549490861 seconds (166123912 bytes allocated, 21.85% gc time)
>
>
>
> 1 abstractarray.jl; checkbounds; line: 62
> 1 array.jl; getindex; line: 247
> 420 task.jl; anonymous; line: 340
> 420 .../IJulia/src/IJulia.jl; eventloop; line: 123
> 420 ...rc/execute_request.jl; execute_request_0x535c5df2; line: 140
> 420 loading.jl; include_string; line: 97
> 420 profile.jl; anonymous; line: 14
> 45 ...qmmm/TBmultiscale.jl; evalSiteFunction!; line: 198
> 1 .../lib/julia/sys.dylib; +; (unknown line)
> 1 array.jl; getindex; line: 247
> 1 ...qmmm/TBmultiscale.jl; evalSiteFunction!; line: 202
> 372 ...qmmm/TBmultiscale.jl; evalSiteFunction!; line: 203
> 3 array.jl; getindex; line: 247
> 8 array.jl; setindex!; line: 308
>
>
