# SymTrididiagonal matrices performance and views allocations

**URL:** <https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618>\
**Category:** Performance\
**Created:** [May 26, 2019, 8:30am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618 "2019-05-26T08:30:48Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [May 26, 2019, 8:30am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/1 "2019-05-26T08:30:48Z")

</div>

Hi,

I work on a toy 2D CFD solver and I wonder about the **idiomatic Julian way** to compute efficiently the product for all `j` \in `[1:ny]`

```julia
Pxy[:,j]=Lx*Sxy[:j]

```

where `Lx` is a a SymTridiagonal matrix of rank `nx` and `Pxy` and `Sxy` are two 2D arrays of size `(nx,ny)`.

I wrote a loop based function `myprod_basic` which seems to perform significantly faster than `myprod`. In addition `myprod` allocates memory:

```julia
 
 3.858 ms (0 allocations: 0 bytes) #myprod_basic
 8.958 ms (200000 allocations: 9.16 MiB) #myprod

```

Note that this snippet may eventually be part of some material for a Julia hands-on session.  
Thank you for your help.

Here is the MWE:

```julia
using BenchmarkTools
using LinearAlgebra
function myprod_basic(Pxy,Sxy,Lx)
    nx,ny=size(Pxy)

    d=Lx.dv
    e=Lx.ev

    for it=1:1000
        # Threads.@threads for j=1:ny
        for j=1:ny
            @inbounds Pxy[1,j]=d[1]*Sxy[1,j]+e[1]*Sxy[2,j]
            for i=2:nx-1
                @inbounds Pxy[i,j]=d[i]*Sxy[i,j]+e[i-1]*Sxy[i-1,j]+e[i]*Sxy[i+1,j]
            end
            @inbounds Pxy[nx,j]=d[nx]*Sxy[nx,j]+e[nx-1]*Sxy[nx-1,j]
        end
    end
end

function myprod(Pxy,Sxy,Lx)
    nx,ny=size(Pxy)
    for it=1:1000
        @inbounds for j=1:ny
            @views mul!(Pxy[:,j],Lx,Sxy[:,j])
        end
    end
end

function tprod(nx,ny)
    dv=rand(nx)
    ev=rand(nx)
    Lx=SymTridiagonal(dv,ev)

    Pxy=zeros(nx,ny)
    Pxyref=zeros(nx,ny)
    Sxy=rand(nx,ny)

    pj=Pxy[:,1]
    sj=Sxy[:,1]

    @btime mul!($pj,$Lx,$sj)
    @btime myprod_basic($Pxyref,$Sxy,$Lx)
    @btime myprod($Pxy,$Sxy,$Lx)

    @assert Pxyref≈Pxy

end

nx=100
ny=100
tprod(nx,ny)

```

---

<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:** [May 26, 2019, 11:23am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/2 "2019-05-26T11:23:02Z")

</div>

> [@LaurentPlagne](#):
>
> I work on a toy 2D CFD solver and I wonder about the **idiomatic Julian way** to compute efficiently the product for all `j` ∈\in `[1:ny]` `Pxy[:,j]=Lx*Sxy[:,j]`

Multiplying a matrix by all of the columns is equivalent to the matrix–matrix multiplication `Lx*Sxy`. So, just do

```julia
mul!(Pxy, Lx, Sxy)

```

to use the pre-allocated output matrix `Pxy`. No views or loops required.

> [@LaurentPlagne](#):
>
> In addition `myprod` allocates memory:

The memory allocation here seems excessive, and it looks like may arise because the `mul!` method for `SymTridiagonal` in the standard library was unnecessarily restricted to act on only a few matrix types. Should be fixed here:

> <https://github.com/JuliaLang/julia/pull/32147>
>
> It looks like several of the methods for tridiagonal methods were unnecessarily …restricted to act on \`StridedVecOrMat\` rather than \`AbstractVecOrMat\`.
> 
> I noticed this because someone reported that \`mul!\` with a view was still allocating lots of memory (https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations)

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [May 26, 2019, 2:54pm UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/3 "2019-05-26T14:54:07Z")

</div>

Thank you very much !  
I did not realized that it can be seen as a Matrix-Matrix product 😊

The timings improve and the allocations disappear but it is still slower than the loop based version

```julia
3.879 ms (0 allocations: 0 bytes) #myprod_basic
8.899 ms (200000 allocations: 9.16 MiB) #myprod
7.663 ms (0 allocations: 0 bytes) #myprod SJ with mul!(Pxy, Lx, Sxy)

```

Looking at the native code, it looks that the latter fails to vectorize.

Thank you for the mul! fix !

Follows the updated MWE.

> **Summary**
>
> ```julia
> using BenchmarkTools
> using LinearAlgebra
> function myprod_basic(Pxy,Sxy,Lx)
> nx,ny=size(Pxy)
> 
> d=Lx.dv
> e=Lx.ev
> 
> for it=1:1000
> # Threads.@threads for j=1:ny
> for j=1:ny
> @inbounds Pxy[1,j]=d[1]*Sxy[1,j]+e[1]*Sxy[2,j]
> for i=2:nx-1
> @inbounds Pxy[i,j]=d[i]*Sxy[i,j]+e[i-1]*Sxy[i-1,j]+e[i]*Sxy[i+1,j]
> end
> @inbounds Pxy[nx,j]=d[nx]*Sxy[nx,j]+e[nx-1]*Sxy[nx-1,j]
> end
> end
> end
> 
> function myprod(Pxy,Sxy,Lx)
> nx,ny=size(Pxy)
> for it=1:1000
> @inbounds for j=1:ny
> @views mul!(Pxy[:,j],Lx,Sxy[:,j])
> end
> end
> end
> 
> function myprod2(Pxy,Sxy,Lx)
> for it=1:1000
> mul!(Pxy,Lx,Sxy)
> end
> end
> 
> function tprod(nx,ny)
> dv=rand(nx)
> ev=rand(nx)
> Lx=SymTridiagonal(dv,ev)
> 
> Pxy=zeros(nx,ny)
> Pxyref=zeros(nx,ny)
> Sxy=rand(nx,ny)
> 
> pj=Pxy[:,1]
> sj=Sxy[:,1]
> 
> @btime mul!($pj,$Lx,$sj)
> @btime myprod_basic($Pxyref,$Sxy,$Lx)
> @btime myprod($Pxy,$Sxy,$Lx)
> 
> @assert Pxyref≈Pxy
> 
> @btime myprod2($Pxy,$Sxy,$Lx)
> @assert Pxyref≈Pxy
> 
> end
> 
> nx=100
> ny=100
> tprod(nx,ny)
> 
> ```

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [May 26, 2019, 3:26pm UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/4 "2019-05-26T15:26:12Z")

</div>

That’s interesting. It seems the julia version (in tridiag.jl) tries to minimize the memory accesses, which interferes with the vectorization. The code is pretty old (@andreasnoack in 2014), so things might have been different then. Looks like you can make a PR with your version and improve things for everybody! Possibly the same applies to the multiplications in bidiag.jl also.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [May 26, 2019, 4:06pm UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/5 "2019-05-26T16:06:33Z")

</div>

Well my biggest concern was about the allocations.  
Before a `SymTridiagonal mul!` PR I should think more about it…

- First: my loop based version lacks some corner case tests like :

```julia
 d=Lx.dv
 e=Lx.ev         
 nx>0 && @inbounds Pxy[1,j]=d[1]*Sxy[1,j]+e[1]*Sxy[2,j]
 for i=2:nx-1
      @inbounds Pxy[i,j]=d[i]*Sxy[i,j]+e[i-1]*Sxy[i-1,j]+e[i]*Sxy[i+1,j]
 end
 nx>1 && @inbounds Pxy[nx,j]=d[nx]*Sxy[nx,j]+e[nx-1]*Sxy[nx-1,j]

```

The 2 extra tests do not seem to cause a noticeable overhead for the 100x100 case.

- Second: I have discussed with Andreas about the opportunity to add a more general matrix type like  
`SymNBanded{T,N}` where `N` would stand for the matrix halfbandwidth (N=1 for Tridiagonal).  
In this case, the product and the solvers (LDLT factorization) can be written once for tri/penta/hepta…diagonal matrices (for higher order discretization schemes).

A few months ago, @ffevotte and I have written such a (proto-)type (François took good care of the unrolling macros) but we still postponed a PR proposal because we (especially I) are still climbing the Julia learning curve 😉

Julia is very easy to enter and makes you want to write scientific software like no other. Nevertheless, it takes time to know enough about Julia to propose library with a quality that can match with existing ones…

---

<div class="post-metadata">

**Author:** ![dkarrasch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dkarrasch/32/7410_2.png) [@dkarrasch](https://discourse.julialang.org/u/dkarrasch)\
**Post date:** [May 26, 2019, 11:08pm UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/6 "2019-05-26T23:08:41Z")

</div>

I’m not sure the `StridedVecOrMat` restriction applies here. When I try

```julia
v = @views Pxy[:,1]
v isa LinearAlgebra.StridedVecOrMat # true
@which mul!(v, Lx, v) # returns the method in tridiag.jl, line 166

```

So there is no fallback acting. Why would a relaxation to `AbstractVecOrMat` avoid the allocations?

---

<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:** [May 27, 2019, 12:15am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/7 "2019-05-27T00:15:09Z")

</div>

> [@dkarrasch](#):
>
> So there is no fallback acting.

Right, sorry: that was a red herring, since `StridedVecOrMat` already includes subarrays.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [May 27, 2019, 5:47am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/8 "2019-05-27T05:47:51Z")

</div>

> [@LaurentPlagne](#):
>
> Second: I have discussed with Andreas about the opportunity to add a more general matrix type like  
> `SymNBanded{T,N}`

You might find BandedMatrices.jl, which supports `Symmetric{T,<:BandedMatrix}` .

---

<div class="post-metadata">

**Author:** ![dkarrasch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dkarrasch/32/7410_2.png) [@dkarrasch](https://discourse.julialang.org/u/dkarrasch)\
**Post date:** [May 27, 2019, 5:52am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/9 "2019-05-27T05:52:54Z")

</div>

I think the allocations are due to generating views. If you divide the allocated memory by the number of allocations, then it’s just `48 bytes` per allocation.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [May 27, 2019, 6:08am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/10 "2019-05-27T06:08:13Z")

</div>

I think that the fixed bandwidth `N` defined in the type is important for performance.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [May 27, 2019, 6:16am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/11 "2019-05-27T06:16:47Z")

</div>

Perhaps, especially with OpenBLAS whose Banded routines are dreadfully slow. But a new Banded type would make a good PR, and you get all the default banded implementations for free.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [May 27, 2019, 6:23am UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/12 "2019-05-27T06:23:44Z")

</div>

Yes, extend BandedMatrices.jl seems logical.

---

<div class="post-metadata">

**Author:** ![dkarrasch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dkarrasch/32/7410_2.png) [@dkarrasch](https://discourse.julialang.org/u/dkarrasch)\
**Post date:** [November 8, 2019, 4:54pm UTC](https://discourse.julialang.org/t/symtrididiagonal-matrices-performance-and-views-allocations/24618/13 "2019-11-08T16:54:10Z")

</div>

I don’t know what has changed recently, but I’m not seeing any allocations in any of the three methods in Julia v1.3, in contrast to v1.2 and below. Anyway, the number of allocations seems to indicate that this is indeed due to the generation of views exclusively. There are 200,000 allocations, which equals 2 (for the 2 views in each `mul!` call) times 100 (for going through the 100 columns) times 1000 (for the number of iterations in `myprod`).