Odd behavior with mul!() and @view

General Usage
1

I’m trying to use mul!() to perform low-memory multiplications on subsets of larger matrices using code. My actual code loops through many slices of smaller views of a large matrix (large set of ODEs chunked for easier understanding), but I’ve reduced the problem to the MWE below.

Is there a way to obtain the same results between the two methods? For reference, I’m using Julia 1.10.

n = 10
m = 8
out = zeros(n, m)
A = ones(n,n)
b = ones(n)
out2 = @view out[:,3:4]
#Method 1 (returns incorrect zero vector)
mul!(out2[:,1], A, b)
#returns out2[:, 1] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

#Method 2 (returns correct nonzero vector)
out2[:,1] = A*b
@show out2[:,1]
#returns out2[:, 1] = [10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0]
2

Method 1 lowers to getindex i.e. makes another slice. Method 2 lowers to setindex!

julia> Meta.@lower mul!(out2[:,1], A, b)
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = Base.getindex(out2, :, 1)
│   %2 = A
│   %3 = mul!(%1, %2, b)
└──      return %3
))))

julia> Meta.@lower out2[:,1] = A*b
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = A * b
│        Base.setindex!(out2, %1, :, 1)
└──      return %1
))))

what is the motivation for using mul! specifically rather than * ? also I think if you use another view this will work as intended

# equivalently, mul!(view(out2, :, 1), A, b)
julia> @views mul!(out2[:,1], A, b)
10-element view(::Matrix{Float64}, :, 3) with eltype Float64:
 10.0
 10.0
 10.0
 10.0
 10.0
 10.0
 10.0
 10.0
 10.0
 10.0
3

Your @views method 3 does work!

My motivation was to reduce memory usage since for this toy example, mul!() uses 112 bytes and A*b uses 160 bytes. As the matrices scale, and since this step is repeated many times my code runs, I thought it would be a more efficient way to do things. At m=n=10000, mul!() still only uses 112 bytes, but A*b uses 78.19KiB. However, they do seem to take about the same amount of time, and it’s not that much memory, so I’m not sure if I’m being silly. In any case, your solution works, and I’ll look more into the difference between setindex!() and getindex() so I understand them better.

Thanks!

4

Just remember that slicing always makes a copy. That’s why

mul!(out2[:,1], A, b)

doesn’t work. You first create a copy of a part of out2 and then pass that to mul! which happily overwrites it, but of course that does not change out2.

Note that on the left side of an assignment like in your version 2

out2[:,1] = A*b

we don’t create a copy to write to but write into out2. The left hand-side of an assignment is the only place where slicing does not create a copy.

1 Like