# Help improving the performance of my implementation of a lower diagonal storage

**URL:** https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867
**Category:** Performance
**Created:** [October 3, 2025, 9:12pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867 "2025-10-03T21:12:57Z")
**Posts on this page:** 16
**Page:** 1

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 3, 2025, 9:12pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/1 "2025-10-03T21:12:57Z")

</div>

Hi!

I am trying to build a lower diagonal storage to reduce the memory footprint of some applications. In my case, when computing spherical harmonics, the coefficients are usual placed in a lower triangular matrix. This is useful for many applications like computing the geomagnetic field and gravity acceleration.

My current implementation is:

```julia
abstract type AbstractDataAlignment end
struct RowMajor <: AbstractDataAlignment end
struct ColumnMajor <: AbstractDataAlignment end

struct LowerTriangularStorage{Ta, Tt} <: AbstractMatrix{Tt}
    @static if VERSION >= v"1.11-"
        data::Memory{Tt}
    else
        data::Vector{Tt}
    end

    n::Int

    function LowerTriangularStorage{Ta, Tt}(n::Int) where {Ta<:AbstractDataAlignment, Tt}
        n < 1 && throw(ArgumentError("Matrix size must be positive"))
        len = (n * (n + 1)) ÷ 2

        data = @static if VERSION >= v"1.11-"
            Memory{Tt}(undef, len)
        else
            Vector{Tt}(undef, len)
        end

        new{Ta, Tt}(data, n)
    end

    function LowerTriangularStorage{Ta}(n::Int) where Ta <: AbstractDataAlignment
        return LowerTriangularStorage{Ta, Float64}(n)
    end

    function LowerTriangularStorage{T}(n::Int) where T
        return LowerTriangularStorage{ColumnMajor, T}(n)
    end

    function LowerTriangularStorage(n::Int)
        return LowerTriangularStorage{ColumnMajor, Float64}(n)
    end
end

@inline Base.@propagate_inbounds function Base.getindex(
    L::LowerTriangularStorage{Ta, Tt},
    i::Int,
    j::Int
) where {Ta <: AbstractDataAlignment, Tt}
    Base.@boundscheck (i < 1 || j < 1 || i > L.n || j > L.n || j > i) &&
        Base.throw_boundserror(L, (i, j))

    return (@inbounds L.data[_axes_to_index(L, i, j)])
end

@inline Base.@propagate_inbounds function Base.setindex!(
    L::LowerTriangularStorage{Ta, Tt},
    v,
    i::Int,
    j::Int
) where {Ta <: AbstractDataAlignment, Tt}
    Base.@boundscheck (i < 1 || j < 1 || i > L.n || j > L.n || j > i) &&
        Base.throw_boundserror(L, (i, j))

    @inbounds L.data[_axes_to_index(L, i, j)] = v
    return L
end

Base.size(L::LowerTriangularStorage) = (L.n, L.n)

function Base.zeros(::Type{T}, n::Int) where T <: LowerTriangularStorage
    L = T(n)
    L.data .= zero(eltype(L.data))
    return L
end

function Base.replace_in_print_matrix(::LowerTriangularStorage, i::Int, j::Int, s::AbstractString)
    i >= j && return s
    s == "#undef" && return " ×" * repeat(" ", length(s) - 3)
    return repeat(" ", length(s))
end

@inline function _axes_to_index(::LowerTriangularStorage{RowMajor}, i::Integer, j::Integer)
    return (i * (i - 1)) ÷ 2 + j
end

@inline function _axes_to_index(L::LowerTriangularStorage{ColumnMajor}, i::Integer, j::Integer)
    return ((j - 1) * (2L.n - j) + 2i) ÷ 2
end

```

For **very large** matrices, we have a huge gain (50%) as we can see by computing the Legendre coefficients:

```julia
julia> using SatelliteToolboxLegendre, BenchmarkTools

julia> P = LowerTriangularStorage{RowMajor}(2200);

julia> dP = LowerTriangularStorage{RowMajor}(2200);

julia> Pm = zeros(2200, 2200);

julia> dPm = zeros(2200, 2200);

julia> legendre!(Val(:full), P, 0.4)

julia> legendre!(Val(:full), Pm, 0.4)

julia> @benchmark dlegendre!(Val(:full), $dPm, 0.4, $Pm)
BenchmarkTools.Trial: 574 samples with 1 evaluation per sample.
 Range (min … max): 8.503 ms … 10.539 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 8.602 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 8.719 ms ± 290.209 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

   ▃▇█▆▃▄▄▂▂▁                                                  
  ▄██████████▆▆█▆▇▇▇▄▅▅▇██▆▄▁▄▄▅▅▄▆▁▁▅▅▅▇▁▁▄▆▁▁▆▅▄▅▆▄▄▁▄▁▆▁▄▄ ▇
  8.5 ms Histogram: log(frequency) by time 9.85 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark dlegendre!(Val(:full), $dP, 0.4, $P)
BenchmarkTools.Trial: 1172 samples with 1 evaluation per sample.
 Range (min … max): 4.156 ms … 5.287 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 4.253 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 4.267 ms ± 83.894 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

         ▃▆▃▆▅█▃                                              
  ▃▄▃▃▄▅█████████▆▅▅▄▄▄▃▄▃▃▃▂▂▂▂▂▂▂▁▁▁▁▂▁▂▁▂▂▁▁▁▁▁▂▂▁▂▂▂▂▁▂▂ ▃
  4.16 ms Histogram: frequency by time 4.65 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

However, for small matrices, the new implementation leads to a performance degradation:

```julia
julia> using SatelliteToolboxLegendre, BenchmarkTools

julia> P = LowerTriangularStorage{RowMajor}(20);

julia> dP = LowerTriangularStorage{RowMajor}(20);

julia> Pm = zeros(20, 20);

julia> dPm = zeros(20, 20);

julia> @benchmark dlegendre!(Val(:full), $dPm, 0.4, $Pm)
BenchmarkTools.Trial: 10000 samples with 285 evaluations per sample.
 Range (min … max): 282.600 ns … 514.326 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 283.042 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 284.933 ns ± 7.142 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

  █▆▂ ▃▁ ▁
  ████▇▆▆▅▅▅▄▄▃▄▄▃▃▃▃▁▃▇██▆▆█▅▆█▆▆▅▆▃▅▆▅▅▅▅▆▆▇█▆▅▆▇▅▅▇▇▅▅▅▅▄▅▅▆ █
  283 ns Histogram: log(frequency) by time 316 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark dlegendre!(Val(:full), $dP, 0.4, $P)
BenchmarkTools.Trial: 10000 samples with 214 evaluations per sample.
 Range (min … max): 349.103 ns … 748.051 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 351.827 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 356.279 ns ± 13.982 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▅▅█▅▁ ▃▂▂▂ ▁
  █████▆▄▅▅▅▆▅▆▄▄▄███████▆▇▆▆▇▆▆▆▆▆▇▇██▇█▇▇█▇▄▅▆▆▅▆▄▆▅▅▃▅▄▄▄▄▆▅ █
  349 ns Histogram: log(frequency) by time 408 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

I am wondering if there is anything I can do to avoid this performance degradation on small matrices so that I can use this new storage type by default.

---

<div class="post-metadata">

### Author: ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [October 3, 2025, 10:06pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/2 "2025-10-03T22:06:45Z")

</div>

Make the size a type parameter and dispatch to whatever is faster?

---

<div class="post-metadata">

### Author: ![DanimirD](https://avatars.discourse-cdn.com/v4/letter/d/d2c977/32.png) [@DanimirD](https://discourse.julialang.org/u/DanimirD)
#### Post date: [October 3, 2025, 10:18pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/3 "2025-10-03T22:18:09Z")

</div>

> [@Ronis\_BR](#):
>
> ` L.data .= zero(eltype(L.data))`

You could check @code\_warntype for this line. I guess there is a few things that could go wrong

---

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 3, 2025, 10:25pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/4 "2025-10-03T22:25:00Z")

</div>

> [@lmiq](#):
>
> Make the size a type parameter and dispatch to whatever is faster?

It can be a solution, but I want to understand what I am doing wrong. I really expect that the performance will increase in any case due to the caching.

> [@DanimirD](#):
>
> You could check @code\_warntype for this line. I guess there is a few things that could go wrong

I checked and it is type stable.

---

<div class="post-metadata">

### Author: ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [October 3, 2025, 10:28pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/5 "2025-10-03T22:28:12Z")

</div>

> [@lmiq](#):
>
> Make the size a type parameter and dispatch to whatever is faster?

I’d probably go with a runtime branch instead for easier type stability and reduced compilation, but figuring out the cutoff isn’t straightforward. It could easily vary by platform. Considering the slowdown is a modest 0.806x at sub-microsecond scales, I question if it’s even worth the effort, especially since there is doubt about the implementation.

> [@DanimirD](#):
>
> You could check @code\_warntype for this line

JET and Cthulhu are recursive, former automatically and the latter manually.

---

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 3, 2025, 10:29pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/6 "2025-10-03T22:29:52Z")

</div>

> [@Benny](#):
>
> I’d probably go with a runtime branch instead for easier type stability and reduced compilation, but figuring out the cutoff isn’t straightforward. It could easily vary by platform. Considering the slowdown is a modest 0.806x at sub-microsecond scales, I question if it’s even worth the effort, especially since there is doubt about the implementation.

Usually this code is called inside a simulator many times per integration step. Hence, this minor degradation can lead to a substantial slowdown.

---

<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: [October 3, 2025, 11:02pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/7 "2025-10-03T23:02:47Z")

</div>

See also the discussion in this thread: [Triangular Arrays in julia? - #27 by mkborregaard](https://discourse.julialang.org/t/triangular-arrays-in-julia/930/27)

My feeling is that the factor of 2 storage benefit is not going to be worth it for the performance degradation of giving up BLAS, _if_ you want to do linear algebra (e.g. triangular solves) with the matrix. YMMV depending on what you want to do with the array.

---

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 3, 2025, 11:09pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/8 "2025-10-03T23:09:24Z")

</div>

Actually I will use it only as storage. Thanks for the thread!

---

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 4, 2025, 2:19pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/9 "2025-10-04T14:19:41Z")

</div>

Now I saw something **very** strange and I really want to know why 🙂

The original code for the derivative of the Legendre coefficients was:

```julia
            if m == 0
                ...
            elseif m == 1
                ...
            elseif n != m
                a_nm = +√(T(n + m) * T(n - m + 1)) / 2
                b_nm = -√(T(n + m + 1) * T(n - m)) / 2

                dP_nm = a_nm * P[i₀ + n, j₀ + m - 1] + b_nm * P[i₀ + n, j₀ + m + 1]

            else
                a_nm = +√(T(n + m) * T(n - m + 1)) / 2

                dP_nm = a_nm * P[i₀ + n, j₀ + m - 1]
            end

```

Now, if I change to:

```julia
            if m == 0
                ...
            elseif m == 1
                ...
            else
                a_nm = +√(T(n + m) * T(n - m + 1)) / 2
                dP_nm = a_nm * P[i₀ + n, j₀ + m - 1]

                if n != m
                    b_nm = -√(T(n + m + 1) * T(n - m)) / 2
                    dP_nm += b_nm * P[i₀ + n, j₀ + m + 1]
                end
            end

```

The performance of both implementations are equal for small matrices. Is this suppose to happen?

---

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 4, 2025, 9:01pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/10 "2025-10-04T21:01:53Z")

</div>

After **a lot** of tests, this simple change in the `if` allowed the new triangular storage to be faster in all scenarios. Its performance computing very big degrees is even better. I just hope that a compiler expert can tell me why this happened 🙂 I will learn a very useful experience for other scenarios.

---

<div class="post-metadata">

### Author: ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [October 4, 2025, 9:10pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/11 "2025-10-04T21:10:55Z")

</div>

I don’t understand the question. Probably you want to provide a reproducer. Also maybe open a new topic.

---

<div class="post-metadata">

### Author: ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [October 4, 2025, 9:48pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/12 "2025-10-04T21:48:13Z")

</div>

Can’t say much without something to tinker with, but it’s plausible at first glance that changed some code motion optimization or compacted the routine to fit in the instruction cache a bit better.

---

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 4, 2025, 9:58pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/13 "2025-10-04T21:58:49Z")

</div>

The problem is that the only way to actually reproduce the problem is copying that definition of the `LowerTriangularStorage` and execute it together with SatelliteToolboxLegendre.jl as I showed in the first post.

---

<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: [October 5, 2025, 1:17am UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/14 "2025-10-05T01:17:32Z")

</div>

I’m skeptical that this is really the _only_ way to see what is going on. But it can take a fair amount of effort to reduce a problem to a minimal working example that illustrates the issue.

---

<div class="post-metadata">

### Author: ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)
#### Post date: [October 5, 2025, 1:32am UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/15 "2025-10-05T01:32:38Z")

</div>

> [@stevengj](#):
>
> I’m skeptical that this is really the _only_ way to see what is going on.

I completely agree! I meant that I was unable to reproduce the problem in any other way 😃 I will keep trying.

---

<div class="post-metadata">

### Author: ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [October 6, 2025, 3:29pm UTC](https://discourse.julialang.org/t/help-improving-the-performance-of-my-implementation-of-a-lower-diagonal-storage/132867/16 "2025-10-06T15:29:55Z")

</div>

Another interesting packed format is Recangular Full Packed (see RectangularFullPacked.jl). It splits a triangular matrix into a trapezoid and triangle, transposes one, and then stacks those into a rectangle. This means it uses the same space as triangular packed formats but has a more BLAS-friendly layout (it looks like it takes advantage of this for solves/inverses but not yet multiplication? could be a PR). Be sure to test the functionality you utilize as the package still seems a touch raw in a few spots.