# \[ANN\] RecursiveFactorization.jl

**URL:** <https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213>\
**Category:** Package Announcements\
**Tags:** package, linearalgebra, recursion\
**Created:** [May 10, 2020, 6:01am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213 "2020-05-10T06:01:43Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [May 10, 2020, 6:01am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/1 "2020-05-10T06:01:43Z")

</div>

RecursiveFactorization allows you to perform any factorization you’d like so long as it is `lu`.

If you like LU, DifferentialEquations.jl, or recursion, you’ll love RecursiveFactorization.jl.

DifferentialEquations.jl uses it as the default method for solving matrices with dimension less than 500 x 500. This package will greatly boost the speed of implicit stiff differential equation solvers, as most of the runtime of implicit solvers is spent on the matrix factorization (the only O(n^3) operation in stiff solvers).

For the following plots, the reference is `LinearAlgebra.generic_lufact!` (a triple loop LU factorization with partial pivoting), `OpenBLAS` is invoked by `LinearAlgebra.lu!` , and `RF` denotes `RecursiveFactorization`. “RF fully iterative” means the base case of the recursion, not to be confused with an [iterative solver](https://github.com/JuliaMath/IterativeSolvers.jl). “RF fully recursive” refers to the algorithm that uses the recursive method for all sizes greater than 10, but the base case calls the optimized triple loop LU for blocks smaller than 16x16. The `threshold` is the cutoff point of using the recursive algorithm. (Sorry for the confusing labels!)

**OpenBLAS 0.3.9**  
**CPU: Intel(R) Xeon(R) Silver 4114 CPU @ 2.20GHz**

 ![](https://global.discourse-cdn.com/julialang/original/3X/3/a/3a4f52d9ee27b57fc5a4ad87c8534fa8e7fb4516.png)

**OpenBLAS 0.3.9**  
**CPU: Intel(R) Core™ i9-10980XE CPU @ 3.00GHz**

 ![](https://global.discourse-cdn.com/julialang/original/3X/c/b/cb9adb7f5a8eea2af6bc61f9902f48292d182590.png)

The first two examples have optimally tuned meta parameter.

**OpenBLAS 0.3.9**  
**CPU: Intel(R) Xeon(R) CPU E5-2680 v3 @ 2.50GHz**

 ![](https://global.discourse-cdn.com/julialang/original/3X/2/5/2571470bf4c9ec3a0e182bf0415053a131f589c6.png)

This example has a badly tuned meta parameter. The default threshold value is 192 here, but one can see the optimal threshold value should be 80. You can override it via

```julia
RecursiveFactorization.RECURSION_THRESHOLD[] = 80

```

These plots were generated via:

```julia
using Pkg
Pkg.add(["RecursiveFactorization", "VegaLite", "DataFrames", "BenchmarkTools"])
using RecursiveFactorization
include(joinpath(pkgdir(RecursiveFactorization), "perf", "lu.jl"))

```

It will save the file in `joinpath(homedir(), "Pictures", "lu_float64.png")`. You want the threshold to be where the fully recursive version passes the fully iterative version.

If you run this, please post the results and your CPU info here, as this will help us develop a better default algorithm!

### Acknowledgment

This package uses @Elrod’s `LoopVectorization.jl`.

---

<div class="post-metadata">

**Author:** ![saschatimme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/saschatimme/32/10313_2.png) [@saschatimme](https://discourse.julialang.org/u/saschatimme)\
**Post date:** [May 10, 2020, 6:49am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/2 "2020-05-10T06:49:00Z")

</div>

Very nice work 🙂 Does the package also work for `Matrix{ComplexF64}`?

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [May 10, 2020, 7:10am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/3 "2020-05-10T07:10:57Z")

</div>

Yes, see [RecursiveFactorization.jl/runtests.jl at a7901459875c5d1ad39e969ce9e16ad764de9557 · JuliaLinearAlgebra/RecursiveFactorization.jl · GitHub](https://github.com/YingboMa/RecursiveFactorization.jl/blob/a7901459875c5d1ad39e969ce9e16ad764de9557/test/runtests.jl#L19) . It is tested for `Float64, Float32, ComplexF64, ComplexF32, Real`. Though I haven’t benchmarked `Matrix{ComplexF64}`, I imagine it will be slower.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [May 10, 2020, 7:24am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/4 "2020-05-10T07:24:02Z")

</div>

FWIW, I ran this script on my laptop:

```julia
:openblas64
Julia Version 1.5.0-DEV.854
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i5-7360U CPU @ 2.30GHz
  JULIA_NUM_THREADS = 4

```

 ![lu_float64_laptop](https://global.discourse-cdn.com/julialang/original/3X/f/b/fb13563831f0677ac73f313af1059cc2b546d929.png)

and desktop:

```julia
:mkl
  CPU: Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz
  JULIA_NUM_THREADS = 12

```

 ![lu_float64_desktop](https://global.discourse-cdn.com/julialang/original/3X/c/e/ce7e468bfb86708f80012e74746200a2a372ca3d.png)

Edit: same desktop on Julia 1.4.1, without MKL:

 ![lu_float64_desktop_1.4](https://global.discourse-cdn.com/julialang/original/3X/7/1/7160726d15bbe7ef51e1002cd2a2f258e0e421b7.png)

I have run into this before, but I’m surprised just how large the gap from openblas is here:

```julia
julia> rr = rand(50,50);

julia> @btime RecursiveFactorization.lu($rr); # on Julia 1.5
  8.313 μs (4 allocations: 20.16 KiB)

julia> @btime lu($rr); # MKL, 1.5
  12.471 μs (4 allocations: 20.16 KiB)

julia> @btime lu($rr); # openblas, 1.4.1
  1.460 ms (4 allocations: 20.16 KiB)

julia> 1460/8.3
175.9036144578313

```

On my laptop (with 1.4 or 1.5) it’s around `20 μs` which is more reasonable.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 10, 2020, 7:46am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/5 "2020-05-10T07:46:34Z")

</div>

On an old desktop

```julia
julia> versioninfo()
Julia Version 1.4.1
Commit 381693d3df (2020-04-14 17:20 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i5-2500K CPU @ 3.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, sandybridge)
Environment:
  JULIA_EDITOR = atom -a
  JULIA_NUM_THREADS = 4
  JULIA_CUDA_VERBOSE = true

```

 ![lu_float64](https://global.discourse-cdn.com/julialang/original/3X/6/6/66f2ca00248d3d05042769d41dc96d04ee9726bb.png)

---

<div class="post-metadata">

**Author:** ![JonasIsensee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonasisensee/32/4704_2.png) [@JonasIsensee](https://discourse.julialang.org/u/JonasIsensee)\
**Post date:** [May 10, 2020, 8:57am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/6 "2020-05-10T08:57:06Z")

</div>

![lu_float64](https://global.discourse-cdn.com/julialang/original/3X/3/0/30796dad76ceb1c87aba44277d943aa5406f88ee.png)

```julia
julia> versioninfo()
Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 2

```

---

<div class="post-metadata">

**Author:** ![zekeriya.sari](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zekeriya.sari/32/13695_2.png) [@zekeriya.sari](https://discourse.julialang.org/u/zekeriya.sari)\
**Post date:** [May 10, 2020, 12:21pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/7 "2020-05-10T12:21:13Z")

</div>

```julia
julia> versioninfo()
Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)

```

 ![lu_float64](https://global.discourse-cdn.com/julialang/original/3X/0/6/06ac689df29e0a05dc108904466d6c2370c761a6.png)

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [May 10, 2020, 4:01pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/8 "2020-05-10T16:01:05Z")

</div>

```julia
julia> versioninfo()
Julia Version 1.4.0
Commit b8e9a9ecc6 (2020-03-21 16:36 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Ryzen Threadripper 2950X 16-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, znver1)

```

 ![lu_float64](https://global.discourse-cdn.com/julialang/original/3X/b/f/bf2ed01c70cd4d5111c7f60c688d9304b24697a1.png)

---

<div class="post-metadata">

**Author:** ![hellemo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hellemo/32/324_2.png) [@hellemo](https://discourse.julialang.org/u/hellemo)\
**Post date:** [May 10, 2020, 5:05pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/9 "2020-05-10T17:05:44Z")

</div>

```julia
julia> versioninfo()
Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8665U CPU @ 1.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 4

```

 ![lu_float64_1.4.1](https://global.discourse-cdn.com/julialang/original/3X/2/f/2fc09baf2db40b46420f55c4b7a98397c53a0beb.png)

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [May 10, 2020, 9:01pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/10 "2020-05-10T21:01:17Z")

</div>

Thank you for running the benchmark!

Note that `RecursiveFactorization` calls [`ldiv!`](https://github.com/YingboMa/RecursiveFactorization.jl/blob/02bd68fe619284d95211cef7db06360c1c5c3f81/src/lu.jl#L48) which in turn invokes `*trtrs` provided by your BLAS. Hence, RecursiveFactorization’s performance is all highly dependent on how well OpenBLAS performs. Unfortunately, from the benchmarks above, it seems that the `threshold` values is highly dependent on the performance of `*trtrs`, as only the recursive algorithm calls `ldiv!`. This is makes `threshold` almost impossible to tune for all computers. The next step might be to implement fast `*trtrs/*trsm` in Julia.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [May 10, 2020, 11:05pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/11 "2020-05-10T23:05:09Z")

</div>

I think it’s also worth pointing out that the primary goal was to perform better under “default settings” for sizes \< 100, which it looks like it is doing fairly consistently, often by large margins.

“Default settings” included not setting `BLAS.set_num_threads(1)`. OpenBLAS is making very poor use of threads here: making it use only a single thread substantially boosted its performance on the computers I tested!  
This is also probably a factor in why OpenBLAS did worse on the higher core-count systems, but did fairly well on laptops and old desktops.

MKL, on the other hand, makes very good use of threads. My tests were consistent with what mcabbott found:

> [@mcabbott](#):
>
> ![lu_float64_desktop](https://global.discourse-cdn.com/julialang/original/3X/c/e/ce7e468bfb86708f80012e74746200a2a372ca3d.png)
> 
> Edit: same desktop on Julia 1.4.1, without MKL:
> 
> ![lu_float64_desktop_1.4](https://global.discourse-cdn.com/julialang/original/3X/7/1/7160726d15bbe7ef51e1002cd2a2f258e0e421b7.png)
> 
> I have run into this before, but I’m surprised just how large the gap from openblas is here:

Also, because RecursiveFactorization is using `*trtrs`, you can see that it actually does better when paired with MKL itself (around 55 vs 45 GFLOPS in the two plots above).

Setting 1 thread will probably shrink the gap at the upper end of the size range in all cases: MKL will no longer be quite so far ahead, and (on systems where it was tripping over its own threads) OpenBLAS won’t be quite so far behind.

Aside from making performance more predictable, so that setting the recursion threshold will be easier, I do think Julia `*trtrs` will further help the package’s primary objective of increasing smaller size (\<100) performance over OpenBLAS and even being competitive with MKL over that range regardless of threading and without the DiffEq ecosystem having to mess with the user’s settings.

---

<div class="post-metadata">

**Author:** ![perseid28](https://avatars.discourse-cdn.com/v4/letter/p/3ab097/32.png) [@perseid28](https://discourse.julialang.org/u/perseid28)\
**Post date:** [June 27, 2020, 9:15am UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/12 "2020-06-27T09:15:00Z")

</div>

Here are my results on a MacBook Pro Ice Lake.

```julia
Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, goldmont)

```

 ![lu_float64](https://global.discourse-cdn.com/julialang/original/3X/b/f/bfcca5b208a74054b0c21adf587c3a6bf193fd8f.png)

And I repeated it with 1.6:

```julia
Julia Version 1.6.0-DEV.306
Commit 59b8dde7c1 (2020-06-26 09:21 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, icelake-client)

```

 ![lu_float64](https://global.discourse-cdn.com/julialang/original/3X/b/c/bce69208ca08702eec702141d7c2d8b5a7c7cf38.png)

---

<div class="post-metadata">

**Author:** ![jebej](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jebej/32/1784_2.png) [@jebej](https://discourse.julialang.org/u/jebej)\
**Post date:** [July 8, 2020, 3:17pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/13 "2020-07-08T15:17:45Z")

</div>

```julia
Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8086K CPU @ 4.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = atom.cmd -a
  JULIA_NUM_THREADS = 6
  JULIA_PKG_SERVER = pkg.julialang.org

```

 ![lu_float64](https://global.discourse-cdn.com/julialang/original/3X/f/a/fa295fbe8cc1d421b31f3164b0d49ff986de1492.png)

---

<div class="post-metadata">

**Author:** ![karei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karei/32/214809_2.png) [@karei](https://discourse.julialang.org/u/karei)\
**Post date:** [February 12, 2026, 1:04pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/14 "2026-02-12T13:04:13Z")

</div>

I just reran the benchmark recently and found that RecursiveFactorization.jl still clearly outperforms the default OpenBLAS path for small–medium matrices on my machine.

On current package versions, the benchmark script needs `VectorizationBase` as an extra dependency. I installed deps via:  
`Pkg.add(["RecursiveFactorization","VegaLite","DataFrames","BenchmarkTools","VectorizationBase"])`

```julia-auto
julia> versioninfo()
Julia Version 1.12.5
Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × Intel(R) Core(TM) i7-10875H CPU @ 2.30GHz
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, skylake)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 16 virtual cores)

julia> LinearAlgebra.BLAS.vendor()
:lbt

julia> BLAS.get_num_threads()
1

```

 ![lu_float64_1.12.5_skylake_1cores_OpenBLAS](https://global.discourse-cdn.com/julialang/original/3X/4/c/4cf8ea0f35a4f900753a8861d4f9e3ec17af0eb9.png)

---

<div class="post-metadata">

**Author:** ![linwaytin](https://avatars.discourse-cdn.com/v4/letter/l/898d66/32.png) [@linwaytin](https://discourse.julialang.org/u/linwaytin)\
**Post date:** [February 13, 2026, 5:40pm UTC](https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213/15 "2026-02-13T17:40:28Z")

</div>

I just know the existence of this package. Its performance is amazing for matrices with size \< 300. More than 50% than the `lu` in LinearAlgebra, which is incredible!