Memory allocation free array slices ?

[Nitin Arora]

Hello,

I realize that the current best way to do array slices is to either use ArrayViews or SubArrays. 

But will it ever be possible (in future) to have slice of contiguous arrays without memory allocation ? How do other compiled languages like Fortran or C deal with this issue ?

For example, given: 

M1 = rand(7,3)
M = view(M1,1:6,1:3) # some memory allocation occurs here
a= rand(3)
c1 = zeros(6)
c =view(c1) # some memory allocation occurs here

would it be possible to do this step without array allocation in future ?

c = M*a # 2 memory allocations

I don't have enough of a background to know the limits of the array indexing/slicing algorithms. I am designing an optimization software in Julia and I am happy using ArrayViews for now but just wanted to know what can be expected in the future.

thanks,

Nitin

[Tim Holy]

I imagine this will be possible some day, but I can't predict when.

Currently, if your slices are contiguous (step size 1 in each range), there is
in fact another option: encode the slice with a CartesianRange. See
http://julialang.org/blog/2016/02/iteration/

Example:

julia> function foo(A, R)
s = zero(eltype(A))
for I in R
s += A[I]
end
s
end
foo (generic function with 1 method)

julia> function bar(A)
s = zero(eltype(A))
for j = 1:size(A,2)
s += foo(A, CartesianRange(CartesianIndex((1,j)),
CartesianIndex((size(A,1),j))))
end
s
end
bar (generic function with 1 method)

julia> A = rand(4,10^4);

julia> bar(A)
20026.547698847495

julia> @time 1
0.000008 seconds (148 allocations: 10.151 KB)
1

julia> @time bar
0.000001 seconds (3 allocations: 144 bytes)
bar (generic function with 1 method)


Look ma, no allocation!

--Tim

[Kevin Squire]

julia> @time bar
  0.000001 seconds (3 allocations: 144 bytes)
bar (generic function with 1 method)

I think this needs to be @time bar(A).  I get

julia> @time bar(A)

  0.000269 seconds (5 allocations: 176 bytes)

20010.937886591404

[Tim Holy]

On Friday, February 05, 2016 08:17:21 AM Kevin Squire wrote:
> I think this needs to be @time bar(A).

Yeah, sorry for the typo.

> I get
>
> julia> @time bar(A)
> 0.000269 seconds (5 allocations: 176 bytes)
> 20010.937886591404

That's just REPL allocation. Since A has 10000 columns, if this strategy were
allocating you'd expect "10000+n allocations," where n comes from the REPL.

--Tim

[Nitin Arora]

Thanks Tim, this is very useful. I will probably use CartesianRange now. 

Is Julia 0.5 Arraypocalypse planning to address this issue as well ?

thanks,

Nitin

[Tim Holy]

On Friday, February 05, 2016 05:24:04 PM Nitin Arora wrote:
> Thanks Tim, this is very useful. I will probably use CartesianRange now.
>
> Is Julia 0.5 Arraypocalypse planning to address this issue as well ?

I don't think there's a way to solve this by changing our implementation of
views; I think it's more of a compiler issue, and the fact that

julia> isbits(CartesianRange((3,5)))
true

julia> immutable ArrayWrapper{M}
data::M
end

julia> z = rand(5,5);

julia> isbits(ArrayWrapper(z))
false


Since I don't work on the compiler, I can't speak for what's going to happen
there, but I'd be surprised if this changes in 0.5.

Best,
--Tim

[Nitin Arora]

I see, thanks for the information. I think, if possible, this feature will help the language a lot. 

I have recommended Julia to many of my colleagues (they all love it over Fortran and Matlab) and most of them seem to run into this issue. I think, a modern language with soo much elegance and potential, like Julia, should nail such major issues. I am sure by V-1.0 we will have the best scientific programming language ever :-)

thanks,

Nitin

[Cedric St-Jean]

M1 = rand(7,3)
M = view(M1,1:6,1:3) # some memory allocation occurs here
a= rand(3)
c1 = zeros(6)
c =view(c1) # some memory allocation occurs here

I don't think that either the M or c line allocates anything (on the heap). ArrayViews are immutables, so they are stack allocated (as far as I can understand) and virtually cost-free.

[Tim Holy]

I filed an issue, https://github.com/JuliaLang/julia/issues/14955, which you
can check to learn about progress on this problem.

Best,
--Tim

[Stefan Karpinski]

Being able to stack-allocate objects that refer to the heap is an important case that we need to address, but doing so is non-trivial and hasn't been done yet.

[Nitin Arora]

Thanks and Tim and Stefan, if there is anyway I can contribute to help, do let me know.

thanks,

Nitin

[mschauer]

A further possibility is to work with vector of FixedSizeArrays (a julia package). You can cast a dxn array to an n-array of immutable d-FixedSizeVectors (they have the same memory layout) and then access the columns directly using stack allocation (just the same as you would work with a complex vector instead of a 2xn array)

[Daan Huybrechs]

Wow, this trick to add CartesianRange's (or other isbits-immutable-iterators I guess) without allocation was very good to learn about, thanks a lot!

Pending more efficient views, this is very helpful to me as well.