A_mul_B! overwrites elements in an array

General Usage
1

I am using DifferentialEquations and I have following preallocated arrays of the type Array{Array{Float64,2},1}:

lay1mul = fill(zeros(n,n),g+1);
views = fill(zeros(n,n),g+1);

Then, in the derivative function I am filling those arrays:

for i=1:1:g+1
  views[i] = view(x,:,:,i)
  A_mul_B!(lay1mul[i], views[i], views[i])
end

(x is generated with ‘rand(n,n,g+1)’ )

I expected lay1mul to contain g+1 possibly different arrays of type Array{Float64,2}.
But all g+1 elements are the same.
It means that running ‘A_mul_B!(lay1mul[1], views[1], views[1])’ fills all elements of lay1mul to the same value.

P.S. I am still using Julia 0.6.4.

2

See https://github.com/JuliaLang/julia/issues/29168, views = fill(zeros(n,n),g+1) creates a vector with the same matrix in each position, e.g.

julia> A = fill(zeros(2,2), 2);

julia> A[1] === A[2]
true

and thus when you change A[1] A[2] (which is the same object) will also change:

julia> A[1][1] = 42; A
2-element Array{Array{Float64,2},1}:
 [42.0 0.0; 0.0 0.0]
 [42.0 0.0; 0.0 0.0]
3 Likes
3

To get a vector of g+1 separate arrays you can use

lay1mul = [zeros(n,n) for _=1:g+1]
1 Like
4

Ow, you are right, that initialization was a problem.

But I have one more question. Looking more at the variables I noticed that

julia> lay1mul[1]===lay1mul[2]
true

But views (after for loop) contains different elements and

julia> views[1]===views[2]
false

What is the reason for this behavior? Is it because after using view() elements of views point at different parts of memory and they are not linked anymore?

@yha Thanks! I’ve been looking for 15 minutes for one-liner to do it.

5

This line is maybe not doing what you are thinking it is doing.
Unless I am mistaken, this is actually doing a deepcopy:

views[i] = view(x,:,:,i)

It is creating a view, but then a new array is created in views, to which the data from the view is deepcopied:

julia> n = 500
500

julia> g = 3
3

julia> x = rand(n,n,g+1);

julia> views = fill(zeros(n,n),g+1);

julia> i = 1
1

julia> @time views[i] = view(x,:,:,i);
  0.000760 seconds (10 allocations: 1.908 MiB)

julia> @time views[i] = view(x,:,:,i);
  0.000387 seconds (10 allocations: 1.908 MiB)

julia> @time v = view(x,:,:,i);
  0.000016 seconds (8 allocations: 304 bytes)

julia> @time v = view(x,:,:,i);
  0.000014 seconds (8 allocations: 304 bytes)

julia> @time views[i] = v;
  0.000451 seconds (6 allocations: 1.908 MiB)

julia> @time views[i] = v;
  0.000561 seconds (6 allocations: 1.908 MiB)

Even if you actually created the views, by doing something like

views = Array{Any}(undef, g+1)

You would still end up with distinct views which might be pointing to the same underlying memory due to the issue listed above.

1 Like