sarah-ek

i have, yes. i can't speak for openblas or mkl, but im familiar with eigen and nalgebra's implementations to some extent

nalgebra doesn't use blocking, so decompositions are handled one column (or row) at a time. this is great for small matrices, but scales poorly for larger ones

eigen uses blocking for most decompositions, other than the eigendecomposition, but they don't have a proper threading framework. the only operation that is properly multithreaded is matrix multiplication using openmp (and the unstable tensor module using a custom thread pool)

because we can do better

fair enough. i try to stay in touch with the eigen and nalgebra developers so i have a good idea on how to write code efficiently with them. for openblas and mkl i've been trying recently to call into the lapack api (benches doing that are still unpublished at the moment), that way im using a universal interface for that kinda stuff.

and of course i do check the cpu utilization to make sure that all threads are spinning for multithreaded benchmarks, and occasionally check the assembly of the hot loops to make sure that the libraries were built properly and are dispatching to the right code. (avx2, avx512, etc) so overall i try to take it seriously and i'll give credit where credit is due when it turns out another library is faster

looks interesting! thanks for sharing

as far as i know, very little is shared. ndarray-linalg is mostly a lapack wrapper. nalgebra and faer both implement the algorithms from scratch, with nalgebra focusing more on smaller matrices

one issue that may be affecting the result is that openmp's threadpool doesn't play well with rayon's. i've seen some perf degradation in the past (on both sides) when both are used in the same program

i plan to address that after refactoring the benches by executing each library individually

funny you mention that, the full pivoting. it's one of the few benchmarks where faer wins by a huge margin

n faer mkl openblas

1024 27.06 ms 186.33 ms 793.26 ms

1536 73.57 ms 605.71 ms 2.65 s

2048 280.74 ms 1.53 s 8.99 s

2560 867.15 ms 3.31 s 17.06 s

3072 1.87 s 6.13 s 55.21 s

3584 3.42 s 10.18 s 71.56 s

4096 6.11 s 15.70 s 168.88 s

the long compile times are mostly because im instantiating every dense decomposition in the library in one translation unit, for several data types (f32, f64, f128, c32, c64, c128)

not yet. a tensor api is on the long term todo list, but it's a big undertaking and i'd like to focus on matrix operations for the time being

i've contributed to eigen in the past and know enough about the internals of the codebase to know my way around safe `auto` usage

im in the process of refactoring the benchmark code at the moment, and plan to include mkl in the benches soon.

overall, the results show that faer is usually faster, or even with openblas, and slower than mkl on my desktop

i did check, yes

author here, eigen is compiled with -fopenmp, which enables parallelism by default