# Improve performance of function that produces and hcats SMatrices

**URL:** <https://discourse.julialang.org/t/improve-performance-of-function-that-produces-and-hcats-smatrices/9140>\
**Category:** Performance\
**Tags:** question, performance\
**Created:** [February 18, 2018, 10:47am UTC](https://discourse.julialang.org/t/improve-performance-of-function-that-produces-and-hcats-smatrices/9140 "2018-02-18T10:47:25Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [February 18, 2018, 10:47am UTC](https://discourse.julialang.org/t/improve-performance-of-function-that-produces-and-hcats-smatrices/9140/1 "2018-02-18T10:47:25Z")

</div>

I have two functions, which for our MWE could be the following:

```julia
using BenchmarkTools, StaticArrays
@inline function f(x, p, n)
    @inbounds x1, x2, x3 = x[1], x[2], x[3]
    SVector( 3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3),
    0.1*( (x2+0.35)*(1-2*x3) - 1 )*(1 - 1.9*x1),
    3.78*x3*(1-x3)+0.2*x2 )
end
@inline function jac(x, p, n)
    @SMatrix [3.8*(1 - 2x[1]) -0.05*(1-2x[3]) 0.1*(x[2] + 0.35);
    -0.19((x[2] + 0.35)*(1-2x[3]) - 1) 0.1*(1-2x[3])*(1-1.9x[1]) -0.2*(x[2] + 0.35)*(1-1.9x[1]);
    0.0 0.2 3.78(1-2x[3]) ]
end

```

I don’t think the actual functions matter, as long as the first returns `SVector` and the latter returns `SMatrix`.  
Let’s do timings:

```julia
p = nothing
s = rand(SVector{3}); Q = rand(SMatrix{3,3})
@btime $f($s, $p, 0)
@btime $jac($s, $p, 0)
  6.158 ns (0 allocations: 0 bytes)
  5.337 ns (0 allocations: 0 bytes)

```

Now, I create a new function, based on the above:

```julia
        ws_index = SVector{3, Int}(2:4...)
        tangentf = (u, p, t) -> begin
            du = f(u[:, 1], p, t)
            J = jac(u[:, 1], p, t)
            dW = J*u[:, ws_index]
            return hcat(du, dW)
        end

```

and benchmark to see:

```julia
S = hcat(s, Q)
@btime $(tangentf)($S, $p, 0)

  109.685 ns (6 allocations: 496 bytes)
3×4 StaticArrays.SArray{Tuple{3,4},Float64,2,12}:
 0.567631 -1.8819 -1.73773 -1.12878 
 0.102564 0.353318 0.263594 0.288192
 0.656533 -1.68341 -0.506125 -1.94484 

```

There are allocations done, and the timing is _much larger_ than the timing of applying the 2 individual functions and performing an `hcat` (which has miniscule timing).

So, how can I improve this and most importantly: **what am I doing wrong**?

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [February 18, 2018, 11:07am UTC](https://discourse.julialang.org/t/improve-performance-of-function-that-produces-and-hcats-smatrices/9140/2 "2018-02-18T11:07:23Z")

</div>

> [@Datseris](#):
>
> `ws_index = SVector{3, Int}(2:4...)`

Don’t use global variables when benchmarking.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [February 18, 2018, 11:10am UTC](https://discourse.julialang.org/t/improve-performance-of-function-that-produces-and-hcats-smatrices/9140/3 "2018-02-18T11:10:01Z")

</div>

> [@kristoffer.carlsson](#):
>
> Don’t use global variables when benchmarking.

Thanks!

Doing instead

```julia
tangentf = let ws_index = SVector{3, Int}(2:4...)
        tangentf = (u, p, t) -> begin
            du = f(u[:, 1], p, t)
            J = jac(u[:, 1], p, t)
            dW = J*u[:, ws_index]
            return hcat(du, dW)
        end
      end

```

```julia

S = hcat(s, Q)
@btime $(tangentf)($S, $p, 0)

  36.128 ns (0 allocations: 0 bytes)
3×4 StaticArrays.SArray{Tuple{3,4},Float64,2,12}:
  0.872362 1.0066 0.474349 0.599955 
 -0.0237808 0.143415 0.0510413 0.0682469
  0.742914 0.813412 1.08132 1.91705  

```

Gives super massive improvements. Still it is not close to the individual function evaluations… Why is that?

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [February 18, 2018, 11:23am UTC](https://discourse.julialang.org/t/improve-performance-of-function-that-produces-and-hcats-smatrices/9140/4 "2018-02-18T11:23:00Z")

</div>

Using a dedicated `struct` gave even more performance boost:

```julia
struct TangentOOP{F, JAC, k}
    f::F
    jacobian::JAC
    ws::SVector{k, Int}
end
@inline function (tan::TangentOOP)(u, p, t)
    du = tan.f(u[:, 1], p, t)
    J = tan.jacobian(u[:, 1], p, t)
    dW = J*u[:, tan.ws]
    return hcat(du, dW)
end

```

benchmarking the later with

```julia
tan = TangentOOP(f, jac, SVector{3, Int}(2:4...))

@btime ($tan)($S, $p, 0);

  22.169 ns (0 allocations: 0 bytes)

```