Can't mul! to a slice

New to Julia
1

After watching Brian Jackson’s talk on JuliaCon 2020 Adventures in Avoiding Allocations, I wanted to try LinearAlgebra.mul! out. But I’m guessing that I don’t get it as test1 and test2 don’t yield the same result. The problem seem to arise when I want to make the “out” variable of mul! be a slice. Can I use mul! to write to an array slice?

julia> import LinearAlgebra

julia> function test1(A, B, n)
           Y = Matrix{Float64}(undef, 2, n)
           LinearAlgebra.mul!(Y[1, :], A, B)
           return Y
       end
test1 (generic function with 1 method)

julia> function test2(A, B, n)
           Y = Matrix{Float64}(undef, 2, n)
           Y[1, :] .= A * B
           return Y
       end
test2 (generic function with 1 method)

julia> n1 = 2
2

julia> n2 = 3
3

julia> A = rand(n2, n1)
3×2 Matrix{Float64}:
 0.425788   0.748112
 0.817301   0.299914
 0.0576634  0.713061

julia> B = ones(n1)
2-element Vector{Float64}:
 1.0
 1.0

julia> Y1 = test1(A, B, n2)
2×3 Matrix{Float64}:
 8.93663e-316  8.93663e-316  8.93663e-316
 8.93663e-316  8.93663e-316  8.93663e-316

julia> Y2 = test2(A, B, n2)
2×3 Matrix{Float64}:
 1.1739  1.11721  0.770724
 0.0     0.0      0.0

julia> [Y1[1, :] Y2[1, :]]
3×2 Matrix{Float64}:
 8.93663e-316  1.1739
 8.93663e-316  1.11721
 8.93663e-316  0.770724
2

While I don’t understand why, the problem seems to go away if I use view to rewrite test1 as

function test1(A, B, n)
    Y = Matrix{Float64}(undef, 2, n)
    LinearAlgebra.mul!(view(Y, 1, :), A, B)
    return Y
end
3

Without the view you allocate a new array that is a copy of the slice of Y, and send this to mul to store the result in. This means Y is never changed in that case. When creating the view you actually use the memory of Y corresponding to that slice to send to mul to be written to.

4 Likes