# Memory Allocation when using mul! with sparse arrays and views

**URL:** <https://discourse.julialang.org/t/memory-allocation-when-using-mul-with-sparse-arrays-and-views/111276>\
**Category:** Performance\
**Tags:** blas, linearalgebra, memory-allocation, sparsearrays\
**Created:** [March 7, 2024, 4:45am UTC](https://discourse.julialang.org/t/memory-allocation-when-using-mul-with-sparse-arrays-and-views/111276 "2024-03-07T04:45:56Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![leespen1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leespen1/32/221928_2.png) [@leespen1](https://discourse.julialang.org/u/leespen1)\
**Post date:** [March 7, 2024, 4:45am UTC](https://discourse.julialang.org/t/memory-allocation-when-using-mul-with-sparse-arrays-and-views/111276/1 "2024-03-07T04:45:56Z")

</div>

I am writing code which performs a five argument `mul!` operation (from `LinearAlgebra`) inside a for loop where performance is critical.

When using a sparse array instead of a matrix, the operation allocates memory:

```julia
julia> using SparseArrays, LinearAlgebra

julia> A = rand(25,25); B = rand(25); C = rand(25);

julia> A_sparse = SparseMatrixCSC(A);

julia> @time mul!(C, A, B, 1, 1);
  0.137919 seconds (236.09 k allocations: 16.206 MiB, 32.65% gc time, 99.98% compilation time)

julia> @time mul!(C, A, B, 1, 1);
  0.000006 seconds

julia> @time mul!(C, A_sparse, B, 1, 1);
  0.514963 seconds (1.77 M allocations: 121.304 MiB, 5.29% gc time, 99.98% compilation time)

julia> @time mul!(C, A_sparse, B, 1, 1);
  0.000007 seconds (2 allocations: 80 bytes)

```

The allocation goes away when the `mul!` is performed in a loop inside a function:

```julia
julia> function f(A, B, C)
         for i in 1:100
           mul!(C, A, B, 1, 1)
         end
       end
f (generic function with 1 method)

julia> @time f(A,B,C);
  0.011613 seconds (3.78 k allocations: 253.445 KiB, 99.70% compilation time)

julia> @time f(A,B,C);
  0.000017 seconds

julia> @time f(A_sparse,B,C);
  0.004931 seconds (6.90 k allocations: 486.984 KiB, 98.89% compilation time)

julia> @time f(A_sparse,B,C);
  0.000039 seconds

julia> @time f(A_sparse,B,C);
  0.000039 seconds

```

However, there **is** allocation when the matrix is sparse and the vectors are views:

```julia
julia> Bview = @view B[:];

julia> Cview = @view C[:];

julia> @time f(A_sparse,Bview,Cview);
  0.320725 seconds (1.36 M allocations: 91.655 MiB, 24.13% gc time, 99.98% compilation time)

julia> @time f(A_sparse,Bview,Cview);
  0.000039 seconds (100 allocations: 4.688 KiB)

```

On the other hand, there **is not** allocation when the matrix is dense and the vectors are views:

```julia
julia> @time f(A,Bview,Cview);
  0.092996 seconds (179.02 k allocations: 12.220 MiB, 99.96% compilation time)

julia> @time f(A,Bview,Cview);
  0.000016 seconds

```

Could someone explain to me why the matrix being sparse and the vectors being views causes memory allocation, and how I could avoid this?

I know `mul!` has warning against some of the arguments aliasing the same memory as other arguments, but that is clearly not happening here.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 7, 2024, 7:33am UTC](https://discourse.julialang.org/t/memory-allocation-when-using-mul-with-sparse-arrays-and-views/111276/2 "2024-03-07T07:33:59Z")

</div>

I think this might be the same as

> [@Why does \`mul!(u, A, v)\` allocate when \`A\` is sparse and \`u, v\` are views?](https://discourse.julialang.org/t/why-does-mul-u-a-v-allocate-when-a-is-sparse-and-u-v-are-views/105995):
>
> I’m wondering if there is a way to avoid allocations in the following scenario: julia\> using BenchmarkTools, LinearAlgebra, SparseArrays julia\> A = sprand(2, 3, 0.5); julia\> u, v = rand(2), rand(3); julia\> @btime mul!($u, $A, $v); 21.139 ns (0 allocations: 0 bytes) julia\> u, v = view(rand(2, 3), :, 1), view(rand(2, 3), 1, :); julia\> @btime mul!($u, $A, $v); 33.731 ns (1 allocation: 48 bytes) If it helps, the size of the allocation does not increase with the size of the arrays

However, note that to benchmark functions properly, you need to use BenchmarkTools.jl and [interpolate variables](https://juliaci.github.io/BenchmarkTools.jl/stable/manual/#Interpolating-values-into-benchmark-expressions) with a dollar sign. Otherwise your allocation count will be biased by the presence of global variables.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 7, 2024, 7:34am UTC](https://discourse.julialang.org/t/memory-allocation-when-using-mul-with-sparse-arrays-and-views/111276/3 "2024-03-07T07:34:53Z")

</div>

Related GitHub issue:

> <https://github.com/JuliaSparse/SparseArrays.jl/issues/469>
>
> This is about in-place multiplication \`mul!(b, A, x)\` of a sparse matrix \`A\` by …a vector \`x\`.
> From 1.9.3 to 1.10.0-rc1, this operation
> 
> \- has gotten significantly slower
> \- has started allocating when \`x\` is a view
> 
> \*\*MWE\*\*
> 
> \`\`\`julia
> using BenchmarkTools, LinearAlgebra, SparseArrays
> 
> function testmul(n)
> A = sparse(Float64, I, n, n)
> b = Vector{Float64}(undef, n)
> @btime mul!($b, $A, x) setup=(x=ones($n))
> @btime mul!($b, $A, x) setup=(x=view(ones($n, 1), :, 1))
> return nothing
> end
> 
> testmul(1000)
> \`\`\`
> 
> \*\*Results\*\*
> 
> | | Julia 1.9 | Julia 1.10 |
> |---|---|---|
> | \`x\` vector | 2.135 μs (0 allocations) | 3.298 μs (0 allocations) |
> | \`x\` view | 2.502 μs (0 allocations) | 4.087 μs (1 allocation) |

Apparently it’s fixed on 1.11 but I haven’t checked

---

<div class="post-metadata">

**Author:** ![leespen1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leespen1/32/221928_2.png) [@leespen1](https://discourse.julialang.org/u/leespen1)\
**Post date:** [March 7, 2024, 5:57pm UTC](https://discourse.julialang.org/t/memory-allocation-when-using-mul-with-sparse-arrays-and-views/111276/4 "2024-03-07T17:57:21Z")

</div>

I think you are right that that is the issue. My call stack looks very similar to the one they show for 1.10. I’ll roll back to 1.9 and see if that fixes is (I would rather not move forward to a version in pre-alpha to fix a bug)

---

<div class="post-metadata">

**Author:** ![natlampen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/natlampen/32/36640_2.png) [@natlampen](https://discourse.julialang.org/u/natlampen)\
**Post date:** [July 8, 2024, 6:30pm UTC](https://discourse.julialang.org/t/memory-allocation-when-using-mul-with-sparse-arrays-and-views/111276/5 "2024-07-08T18:30:07Z")

</div>

One thing I discovered today which might help you. MKLSparse.jl removed the allocations for me. Also it provided a close to 2x speed up for my usecase.