# Avoiding unnecessary mallocs while building then solving a linear system inside a loop can only be achieved by routines that modify in-place?

**URL:** <https://discourse.julialang.org/t/avoiding-unnecessary-mallocs-while-building-then-solving-a-linear-system-inside-a-loop-can-only-be-achieved-by-routines-that-modify-in-place/29210>\
**Category:** Performance\
**Tags:** memory-allocation\
**Created:** [September 26, 2019, 9:22pm UTC](https://discourse.julialang.org/t/avoiding-unnecessary-mallocs-while-building-then-solving-a-linear-system-inside-a-loop-can-only-be-achieved-by-routines-that-modify-in-place/29210 "2019-09-26T21:22:53Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![pedrohnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pedrohnv/32/212503_2.png) [@pedrohnv](https://discourse.julialang.org/u/pedrohnv)\
**Post date:** [September 26, 2019, 9:22pm UTC](https://discourse.julialang.org/t/avoiding-unnecessary-mallocs-while-building-then-solving-a-linear-system-inside-a-loop-can-only-be-achieved-by-routines-that-modify-in-place/29210/1 "2019-09-26T21:22:53Z")

</div>

Hello,

I have some functions that build a linear system. I then want to build and solve the system in a loop, while modifying a parameter. My first, straightforward, approach resulted in many allocations being done unnecessarily. It was as follows:

```julia
function f(n::Integer)
    x = Array{Float64}(undef, n, n)
    @inbounds for i = 1:n
        @inbounds for k = 1:n
            x[k,i] = randn()
        end
    end
    return x
end

function g(n::Integer, m::Integer)
    y = Array{Float64}(undef, m)
    b = zeros(Float64, n)
    b[1] = 1.0
    @inbounds for i = 1:m
        a = f(n)
        x = a\b
        y[i] = x[1]
    end
    return y
end

```

The profiling:

```julia
@time g(2, 4) # 0.000021 seconds (26 allocations: 2.109 KiB)
@time g(2, 5) # 0.000025 seconds (31 allocations: 2.563 KiB)

```

The only way I could think of that would avoid those unnecessary allocations was to use functions that modify the inputs in-place, as follows:

```julia
using LinearAlgebra

function f!(x, n::Integer)
    for i = 1:n
        for k = 1:n
            x[k,i] = randn()
        end
    end
end

function g!(y, n::Integer, m::Integer)
    a = Array{Float64}(undef, n, n)
    b = Array{Float64}(undef, n)
    for i = 1:m
        b[1] = 1.0
        for k = 2:n
            b[k] = 0.0
        end
        f!(a, n)
        LAPACK.gesv!(a,b)
        y[i] = b[1]
    end
end

```

The profiling:

```julia
y = zeros(Float64, 5)
@time g!(y, 2, 4) # 0.000018 seconds (10 allocations: 752 bytes)
@time g!(y, 2, 5) # 0.000021 seconds (11 allocations: 848 bytes)

```

Is there a way to optimize the first approach with a `f` that does not need to modify the `a` matrix in-place? I am concerned with the users of my code: engineers with weak programming skills who just know a bit of MATLAB. I fear that the second approach will be too alien to them, weakening my efforts to proselytize Julia.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [September 26, 2019, 9:51pm UTC](https://discourse.julialang.org/t/avoiding-unnecessary-mallocs-while-building-then-solving-a-linear-system-inside-a-loop-can-only-be-achieved-by-routines-that-modify-in-place/29210/2 "2019-09-26T21:51:50Z")

</div>

One thing you could do would be to use the `@!` macro from [https://github.com/simonbyrne/InplaceOps.jl](https://github.com/simonbyrne/InplaceOps.jl) to let you write `a \ b` but have the actual code run in-place. Whether this is more or less clear to your users is hard to say: a macro can make the syntax look more like what they would expect, but it also introduces a new layer of indirection that they may eventually have to understand.

On the other hand, I think your current code is actually fine: there is no magic in it, and it accomplishes its task efficiently. Perhaps your effort would be better spent writing comments to help your readers understand. After all, Julia and Matlab are quite different languages, and good Julia code will usually end up looking quite different from Matlab code.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 27, 2019, 2:50am UTC](https://discourse.julialang.org/t/avoiding-unnecessary-mallocs-while-building-then-solving-a-linear-system-inside-a-loop-can-only-be-achieved-by-routines-that-modify-in-place/29210/3 "2019-09-27T02:50:19Z")

</div>

You can do in-place matrix operations with the `lu!` and `ldiv!` functions from the `LinearAlgebra` module, without dropping down to “raw” LAPACK calls. You can also use in-place “dot assignments” like `.=`. For example, your function could be rewritten:

```julia
function g!(y, n::Integer, m::Integer)
    a = Array{Float64}(undef, n, n)
    b = Array{Float64}(undef, n)
    for i = 1:m  
        b .= 0
        b[1] = 1
        a .= randn.()
        ldiv!(lu!(a), b)
        y[i] = b[1]
    end
end

```

However, if you are doing zillions of solves and other operations with tiny arrays, e.g. the 2x2 matrices in your code, it should generally be _much_ faster to use the `StaticArrays` package, which both avoids heap allocation and has optimized inlined implementations for such small sizes. `StaticArrays` take a bit of getting used to in order to get the full benefit (especially to make sure the dimensions are known when a critical function is compiled), however.

PS. Better to use `@btime` from the BenchmarkTools package, e.g. `@btime g!($y, 2, 5)`, for benchmarking, especially for benchmarking such inexpensive functions.