You can do a lot better by dispatching on the dimension of A. Something like
this:
mygetindex{T}(A::AbstractArray{T,2}, L, Lcol) = A[L[1,Lcol], L[2,Lcol]]
mygetindex{T}(A::AbstractArray{T,3}, L, Lcol) = A[L[1,Lcol], L[2,Lcol],
L[3,Lcol]]
Put your timing loops in a function so you take advantage of inference. For
example:
function accesseach(A, L)
dummy = 0.0
for Lcol = 1:size(L, 2)
dummy = mygetindex(A, L, Lcol)
end
dummy
end
julia> L = Array(Int, 3, N^3);
julia> idx = 1; for k = 1:N, j = 1:N, i = 1:N
L[1,idx] = i
L[2,idx] = j
L[3,idx] = k
idx += 1
end
julia> accesseach(A, L)
0.4109242495731802
julia> @time accesseach(A, L)
elapsed time: 0.006461321 seconds (13824 bytes allocated)
0.4109242495731802
As you can see, quite a bit faster.
--Tim
On Monday, November 17, 2014 12:01:08 PM Christoph Ortner wrote:
> I would appreciate advise how to best implement the following:
>
> I have an N^d array A, where the dimension d depends on the
> application (d \in \{1, 2, 3\}). Somewhere else I have a list L which is a
> d x M array or integers corresponding to points/elements in this array,
> e.g. if
> L[:, 3] == [5, 2, 8]
> then A would be N^3 and
> M[L[1,3], L[2,3], L[3,3]]
> would read
> A[5, 2, 8]
>
> The difficulty is that I would like to read from A in a *fast*
> dimension-independent way. My only idea was to define a type and overload
> getindex, setindex! ; see below. If I tested this correctly, then I only
> loose a factor of two in terms of both speed and memory.
>
> I am slightly worried that it is a non-standard access to a standard array.
> + I'd love to get that factor 2 back. Hence, I'd like to know whether
> there are faster / more elegant alternatives, or just other ways of doing
> this?
>
>
>
> module damod
>
>
>
> export darray
>
> immutable type darray
>
> data
>
> dim
>
> end
>
>
>
> function darray(A)
>
> return darray(A, length(size(A)))
>
> end
>
>
>
> export getindex
>
> function getindex(A::darray, ii)
>
> if A.dim == 2
>
> return A.data[ii[1], ii[2]]
>
> elseif A.dim == 3
>
> return A.data[ii[1], ii[2], ii[3]]
>
> end
>
> end
>
> end
>
>
>
> using damod
>
> N = 100
>
> a = darray(rand(N, N, N))
>
> @time for n = 1:N, m=1:N, k=1:N; dummy = a[[n,m,k]] end
>
> @time for n = 1:N, m=1:N, k=1:N; dummy = a.data[n,m,k] end
>
> @time for n = 1:N, m=1:N, k=1:N; dummy = a[[n,m,k]] end
>
> @time for n = 1:N, m=1:N, k=1:N; dummy = a.data[n,m,k] end
>
> @time for n = 1:N, m=1:N, k=1:N; dummy = a[[n,m,k]] end
>
> @time for n = 1:N, m=1:N, k=1:N; dummy = a.data[n,m,k] end
>
>
>
>
>
>
>
> elapsed time: 0.470031905 seconds (128633244 bytes allocated, 11.90% gc
> time) elapsed time: 0.209358484 seconds (48565624 bytes allocated, 8.98% gc
> time) elapsed time: 0.42617227 seconds (128641396 bytes allocated, 16.12%
> gc time) elapsed time: 0.21529278 seconds (48565624 bytes allocated, 10.65%
> gc time) elapsed time: 0.457094657 seconds (128565624 bytes allocated,
> 14.37% gc time) elapsed time: 0.201969234 seconds (48565624 bytes
> allocated, 10.24% gc time)
