# Ldiv for Cholesky is slower than two substitutions

**URL:** <https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092>\
**Category:** Performance\
**Created:** [June 21, 2025, 3:37pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092 "2025-06-21T15:37:26Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [June 21, 2025, 3:37pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/1 "2025-06-21T15:37:26Z")

</div>

I’ve faced that `ldiv` for `Cholesky` is unexpectedly slower than two substitutions.

Here is a test problem with a random 1000-by-1000 `Matrix{Float64}`.

```julia-repl
julia> using BenchmarkTools, LinearAlgebra

julia> n = 1000; A = rand(n, n); A = A' * A; b = rand(n);

```

Results for `ldiv` on factorization object.

```julia-repl
julia> Achol = cholesky(A);

julia> @benchmark $Achol \ $b
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 315.875 μs … 531.959 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 317.833 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 319.034 μs ± 4.421 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

   ▄██▇▇▇
  ▃██████▇▃▁▂▄▆▅▇█▆▄▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  316 μs Histogram: frequency by time 336 μs <

 Memory estimate: 8.00 KiB, allocs estimate: 1.

```

Results for two substitutions (factors of `Achol` are `Upper`/`LowerTriangular`).

```julia-repl
julia> norm(Achol.U \ (Achol.L \ b) - Achol \ b, Inf)
1.7053025658242404e-12

julia> @benchmark U \ (L \ $b) setup=(L=$Achol.L; U=$Achol.U)
BenchmarkTools.Trial: 4178 samples with 1 evaluation.
 Range (min … max): 189.958 μs … 341.667 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 204.875 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 205.505 μs ± 7.942 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

             ▂▅▇▆██▆▇▇▅▄▆▅▄▄▄▄▆▄▂▆▃▁▅▅▁▃▅▂▃▂▂▄▁▁
  ▂▁▂▃▂▃▃▄▆▆████████████████████████████████████▆▆▄▄▃▃▂▃▁▂▁▂▁▂▁ ▅
  190 μs Histogram: frequency by time 224 μs <

 Memory estimate: 16.00 KiB, allocs estimate: 2.

```

Performance of the above example is close to `ldiv` for `LU`.

```julia-repl
julia> Alu = lu(A);

julia> @benchmark $Alu \ $b
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 176.708 μs … 451.208 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 178.750 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 188.390 μs ± 32.983 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  █▅▂ ▃▃▃ ▁
  ████▇█▇▆▄▃▃▁▁▃▁▁███▇▃▃▅▁▄▃▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▆▆▃▅▆▇ █
  177 μs Histogram: log(frequency) by time 400 μs <

 Memory estimate: 8.00 KiB, allocs estimate: 1.

```

The examples above I run on my MacBook. Here is corresponding `versioninfo`.

```julia-repl
julia> versioninfo()
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 10 × Apple M2 Pro
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 6 virtual cores)

```

But, I’ve faced similar situation with the following setups (can’t provide `versioninfo` for them)

- Julia 1.10, Ubuntu 22.04, build for Generic Linux on x86 (so OpenBLAS should be under the hood)
- Julia 1.11 on Windows 10, 64-bit installer

I expect that `Achol \ b` should be faster than `Alu \ b`. Do I measuring performance wrong or `Achol \ b` is buggy?

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [June 21, 2025, 5:02pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/2 "2025-06-21T17:02:10Z")

</div>

Using the code below, I get that the former is much faster than the second for n=1000. The difference is that you’re using the same matrix every time.

```julia
function cholera( n )
           A = rand( n, n )
           C = cholesky( A'A )
           b = rand( n )
           C, b
end

f1( Achol, b ) = Achol \ b

f2( Achol, b ) = Achol.U \ (Achol.L \ b )

@benchmark f1(C,b) setup=( (C,b) = cholera(n); )

@benchmark f2(C,b) setup=( (C,b) = cholera(n); )

```

```julia
BenchmarkTools.Trial: 293 samples with 1 evaluation per sample.
 Range (min … max): 574.176 μs … 751.006 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 601.859 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 606.899 μs ± 27.591 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

      ▂▃ ▅▇█▃▃▂▂                                                
  ▃▁▃████▇███████▆▅▂▃▃▃▃▃▂▂▁▁▁▂▁▂▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▁▂▂▁▃ ▃
  574 μs Histogram: frequency by time 739 μs <

```

```julia
BenchmarkTools.Trial: 246 samples with 1 evaluation per sample.
 Range (min … max): 2.490 ms … 8.196 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 3.565 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 4.447 ms ± 1.719 ms ┊ GC (mean ± σ): 11.05% ± 16.52%

  ▅▂ ▂▆ █ ▁▂ ▂▄  
  ██▁▁▁▁▁▁██▆██▄▄▁▁▁▁▁▁▁▁▄▄▁▁▁▇▆███▄▆▇▆▄▁▁▁▁▇▇█▁▁▁▁▁▄▁▄▄▆██ ▆
  2.49 ms Histogram: log(frequency) by time 7.79 ms <

```

---

<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:** [June 21, 2025, 5:11pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/3 "2025-06-21T17:11:09Z")

</div>

> [@Joris\_Pinkse](#):
>
> `f2( Achol, b ) = Achol.U \ (Achol.L \ b )`

Note that this allocates a matrix copy (in Julia 1.11) to obtain the`Achol.L` matrix each time, and it also allocates additional intermediate vectors, which will slow things down.

The [Cholesky `ldiv !` method](https://github.com/JuliaLang/LinearAlgebra.jl/blob/b9d884339efcb08d4cbe4c7b6503557fc6ae5ba4/src/cholesky.jl#L739) callsl the [LAPACK `dpotrs`](https://netlib.org/lapack/explore-html-3.6.1/d1/d7a/group__double_p_ocomputational_ga167aa0166c4ce726385f65e4ab05e7c1.html), whereas the [`UpperTriangular` and `LowerTriangular`](https://github.com/JuliaLang/LinearAlgebra.jl/blob/b9d884339efcb08d4cbe4c7b6503557fc6ae5ba4/src/triangular.jl#L1280-L1286) methods for `ldiv!` call the [LAPACK `dtrtrs`](https://netlib.org/lapack/explore-html/d4/dc1/group__trtrs_gab0b6a7438a7eb98fe2ab28e6c4d84b21.html) function.

(You might also want to turn off BLAS multi-threading, as that tends to make benchmark results more confusing.)

---

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [June 22, 2025, 6:38am UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/5 "2025-06-22T06:38:11Z")

</div>

Oh, it seems that LAPACK routines are optimized for matrix-matrix system `A X = B`, instead of matrix-vector `A x = b`. Performance for the matrix-matrix case looks correct.

```julia-repl
julia> using BenchmarkTools, LinearAlgebra

julia> Threads.nthreads()
1

julia> n = 1000; A = rand(n, n); A = A' * A; B = rand(n, n);

julia> Achol = cholesky(A); Alu = lu(A);

julia> @btime $Achol \ $B;
  8.391 ms (2 allocations: 7.63 MiB)

julia> @btime $Alu \ $B;
  8.499 ms (2 allocations: 7.63 MiB)

julia> @btime U \ (L \ $B) setup=(L=$Achol.L; U=$Achol.U);
  8.774 ms (4 allocations: 15.26 MiB)

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [June 22, 2025, 4:31pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/6 "2025-06-22T16:31:24Z")

</div>

> [@stepanzh](#):
>
> `julia> @benchmark U \ (L \ $b) setup=(L=$Achol.L; U=$Achol.U)`

To elaborate on what @stevengj said: You’re pulling the most costly part of the computation out into the setup here, so you’re only measuring a small part of the actual cost of the operation. In a fair benchmark, `Achol \ b` wins handily. Here’s what I find:

```julia-repl
julia> @benchmark $Achol \ $b
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max): 248.917 μs … 388.625 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 261.416 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 263.493 μs ± 9.988 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

   ▁▂ ▃██▄▂▃▃▃▂▁
  ▃██▇▆▆████████████▆▇▇▇▇▆▆▆▇▆▇▇▇▇▇▇▇▇▇▆▅▅▅▆▅▅▄▄▄▄▃▃▂▂▂▃▂▂▂▂▂▂▂ ▄
  249 μs Histogram: frequency by time 289 μs <

 Memory estimate: 8.00 KiB, allocs estimate: 1.

julia> @benchmark U \ (L \ $b) setup=(U = $Achol.U; L = $Achol.L) # Unfair
BenchmarkTools.Trial: 5858 samples with 1 evaluation per sample.
 Range (min … max): 132.750 μs … 227.459 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 145.292 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 145.427 μs ± 5.723 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

             ▁▁▄▄▄▄▄▄▃▂▄▄▄▄▆▅█▆▅▄▆▄▄▅▅▂▃
  ▁▁▁▁▂▃▄▅▄▇▇███████████████████████████▇▇▅▆▄▃▃▃▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁ ▅
  133 μs Histogram: frequency by time 162 μs <

 Memory estimate: 16.00 KiB, allocs estimate: 2.

julia> myldiv(chol, x) = chol.U \ (chol.L \ x)
myldiv (generic function with 1 method)

julia> @benchmark myldiv($Achol, $b) # Fair, but slower
BenchmarkTools.Trial: 5796 samples with 1 evaluation per sample.
 Range (min … max): 659.125 μs … 2.103 ms ┊ GC (min … max): 0.00% … 57.96%
 Time (median): 765.250 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 862.006 μs ± 209.153 μs ┊ GC (mean ± σ): 9.99% ± 14.68%

    ▁▂▂▆██▆▄▄▄▄▃▂▂▂▁▁ ▁▁▂▂▃▃▃▂▂▂▁▁▁ ▂
  ▆▆█████████████████▇▆▇▅▁▄▄▃▃▁▁▁▁▁▁▁▁▄▆████████████████▇█▆▆▆▅▄ █
  659 μs Histogram: log(frequency) by time 1.5 ms <

 Memory estimate: 7.65 MiB, allocs estimate: 4.

```

---

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [June 22, 2025, 5:25pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/7 "2025-06-22T17:25:23Z")

</div>

> [@danielwe](#):
>
> You’re pulling the most costly part of the computation out into the setup here

Let me clear the context. For actual solver I have “infinite” time budget to prepare, but limited time budget for real-time computations (offline/online computations). E.g. for me it’s fair to compute `cholesky(A)` or `Achol.L`/`Achol.U` and store it on a hard drive. I’m interested only in performance of real-time computations of `A \ b`. That’s why putting `Achol.L` and `Achol.U` into `setup` is valid and no-cost for me.

---

<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:** [June 22, 2025, 5:31pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/8 "2025-06-22T17:31:32Z")

</div>

> [@stevengj](#):
>
> The [Cholesky `ldiv !` method](https://github.com/JuliaLang/LinearAlgebra.jl/blob/b9d884339efcb08d4cbe4c7b6503557fc6ae5ba4/src/cholesky.jl#L739) callsl the [LAPACK `dpotrs`](https://netlib.org/lapack/explore-html-3.6.1/d1/d7a/group__double_p_ocomputational_ga167aa0166c4ce726385f65e4ab05e7c1.html), whereas the [`UpperTriangular` and `LowerTriangular`](https://github.com/JuliaLang/LinearAlgebra.jl/blob/b9d884339efcb08d4cbe4c7b6503557fc6ae5ba4/src/triangular.jl#L1280-L1286) methods for `ldiv!` call the [LAPACK `dtrtrs`](https://netlib.org/lapack/explore-html/d4/dc1/group__trtrs_gab0b6a7438a7eb98fe2ab28e6c4d84b21.html) function.

To be clear, the source code of the reference `dpotrs` routine simply calls `dtrsm` twice, whereas `dtrtrs` also calls `dtrsm`. So, in principle, there should be no advantage to doing the triangular solves yourself.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [June 22, 2025, 7:44pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/9 "2025-06-22T19:44:16Z")

</div>

> [@stepanzh](#):
>
> I have “infinite” time budget to prepare

Turns out it’s faster to avoid the extra preparation anyway. You can use an adjoint wrapper instead of materializing `L`, and the resulting computations are even faster:

```julia-repl
julia> function myldiv(chol, x)
           if chol.uplo == 'U'
               return ldiv!(chol.U, chol.U' \ x)
           elseif chol.uplo == 'L'
               return ldiv!(chol.L', chol.L \ x)
           else
               error("unreachable")
           end
       end
myldiv (generic function with 1 method)

julia> @benchmark myldiv($Achol, $b)
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max): 120.500 μs … 276.959 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 129.959 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 130.370 μs ± 8.761 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▂▁ █
  ██▆▄▅▅▅▄▄▄▃▃▃▃▃▃▃▄▄██▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂ ▃
  120 μs Histogram: frequency by time 166 μs <

 Memory estimate: 8.00 KiB, allocs estimate: 1.

```

Compare this to my timings above for the `U \ (L \ $b)` version on the same computer:

```julia-repl
 Range (min … max): 132.750 μs … 227.459 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 145.292 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 145.427 μs ± 5.723 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

             ▁▁▄▄▄▄▄▄▃▂▄▄▄▄▆▅█▆▅▄▆▄▄▅▅▂▃
  ▁▁▁▁▂▃▄▅▄▇▇███████████████████████████▇▇▅▆▄▃▃▃▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁ ▅
  133 μs Histogram: frequency by time 162 μs <

 Memory estimate: 16.00 KiB, allocs estimate: 2.

```

(Note that I used `ldiv!` to avoid an extra allocation in the new version, but that doesn’t significantly affect the timings, so the comparison is still valid. The O(N) allocation cost is dwarfed by the O(N^2) computation cost at this size.)

---

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [June 22, 2025, 8:07pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/10 "2025-06-22T20:07:12Z")

</div>

Thanks for this snippet! It works for me. In your environment, `Achol \ b` is slower than `myldiv(Achol, b)`?

---

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [June 22, 2025, 8:11pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/11 "2025-06-22T20:11:05Z")

</div>

Considering answer of @danielwe…

> [@Ldiv for Cholesky is slower than two substitutions](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/9):
>
> Turns out it’s faster to avoid the extra preparation anyway. You can use an adjoint wrapper instead of materializing L, and the resulting computations are even faster: julia\> function myldiv(chol, x) if chol.uplo == 'U' return ldiv!(chol.U, chol.U' \ x) elseif chol.uplo == 'L' return ldiv!(chol.L', chol.L \ x) else error("unreachable") end end myldiv (generic function with 1 method) julia\> @benchma…

Should `LinearAlgebra` have a method `ldiv!(::Cholesky, ::StridedVector)`?

---

<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:** [June 22, 2025, 8:14pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/12 "2025-06-22T20:14:08Z")

</div>

> [@stepanzh](#):
>
> Thanks for this snippet! It works for me. In your environment, `Achol \ b` is slower than `myldiv(Achol, b)`?

Yes, on mine (for the 1000x1000 test case), benchmarking with 1 thread:

```julia
julia> BLAS.set_num_threads(1)

julia> @btime $Achol \ $b;
  314.125 μs (3 allocations: 8.06 KiB)

julia> @btime myldiv($Achol, $b);
  147.167 μs (3 allocations: 8.06 KiB)

```

Seems like an issue in OpenBLAS.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [June 22, 2025, 8:23pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/13 "2025-06-22T20:23:53Z")

</div>

I think the difference here is between OpenBLAS’ `trsv` and `trsm`. As @stevengj pointed out, `Achol \ b` will dispatch to `LAPACK.potrs!` which calls `trsm`. While [reference LAPACK’ trtrs](https://www.netlib.org/lapack/explore-html/d4/dc1/group__trtrs_gab0b6a7438a7eb98fe2ab28e6c4d84b21.html#gab0b6a7438a7eb98fe2ab28e6c4d84b21) simply calls `trsm`, this is not the case for OpenBLAS’ trtrs [which branches on n == 1](https://github.com/OpenMathLib/OpenBLAS/blob/b4945057b7a72483538b20cb85d96e7d9491969e/lapack/trtrs/trtrs_single.c#L70-L74). So the main difference in performance for the two benchmarked solutions seems to be due to differences in OpenBLAS’ `trsv` and `trsm` which is confirmed below. I ran this on an Intel Mac. It would be interesting to see these numbers for MKL.

```julia
julia> @benchmark BLAS.trsv!('U', 'N', 'N', $(Achol).factors, copy($b))
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max): 77.639 μs … 311.568 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 83.043 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 87.832 μs ± 14.097 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▁▃▆▇█▇▅▃▂▄▄▄▂▁▃▃▂▁ ▁▁ ▁ ▁ ▂
  ███████████████████▇▆▇▅▆▇▇▇██▇██▇▇▅▅█▆▅▆▇▄▃▆██▅▆██▆▅▅▄▃▄▃▄▄▄ █
  77.6 μs Histogram: log(frequency) by time 144 μs <

 Memory estimate: 8.06 KiB, allocs estimate: 3.

julia> @benchmark BLAS.trsm!('L', 'U', 'N', 'N', 1.0, $(Achol).factors, copy($B))
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max): 245.243 μs … 855.663 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 267.437 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 276.065 μs ± 40.048 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▆▃█▆▅▃▇▆▆▄▁▁▃▃▂▁▁▂▃▂▂▃▂▁▂ ▁ ▂ ▂
  ████████████████████████████████▆▇▇▇▇█▆▇▆▆▇▄▅▆▅▅▇▅▅▆▅▄▆▃▃▄▃▅▄ █
  245 μs Histogram: log(frequency) by time 445 μs <

 Memory estimate: 8.08 KiB, allocs estimate: 3.

```

(I should mentioned that `B = reshape(b, length(b), 1)`).

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [June 22, 2025, 8:41pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/14 "2025-06-22T20:41:26Z")

</div>

> [@stepanzh](#):
>
> In your environment, `Achol \ b` is slower than `myldiv(Achol, b)`?

Yes. Here are my timings for `Achol \ b` (copied from my original post above). `myldiv` is about twice as fast.

```julia-repl
 Range (min … max): 248.917 μs … 388.625 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 261.416 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 263.493 μs ± 9.988 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

```

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [July 2, 2025, 9:09pm UTC](https://discourse.julialang.org/t/ldiv-for-cholesky-is-slower-than-two-substitutions/130092/15 "2025-07-02T21:09:15Z")

</div>

For me it’s 399 versus 112, so the latter is almost four times faster. (AMD, 1 BLAS thread).