"Fast" solvers are slow for dense, complex matrix

[simonp0420]

According to the documentation, RFLUFactorization is supposed to be faster than backslash for matrices on the order of a few hundred. I did some quick benchmarking and I'm not finding this to be the case. Am I doing something wrong?

using LinearSolve, BenchmarkTools

n = 200
A = rand(ComplexF64, n, n)
b = rand(ComplexF64, n)

@btime $A\$b ; # 474 μs

prob = LinearProblem(copy(A), copy(b))
@btime solve($prob, LUFactorization()); # 477 μs. As expected, similar to `\`

prob = LinearProblem(copy(A), copy(b))
@btime solve($prob); # 20.4 ms.  Not sure what method is used here.

prob = LinearProblem(copy(A), copy(b))
@btime solve($prob, RFLUFactorization()); # 21.8 ms.  Forty-some times slower than BLAS LU

prob = LinearProblem(copy(A), copy(b))
method = LinearSolve.FastLUFactorization()
@btime solve($prob, $method); # 14.9 ms .  Better but still way slower than BLAS

using LinearAlgebra
BLAS.get_config()
#=
LinearAlgebra.BLAS.LBTConfig
Libraries: 
└ [ILP64] libopenblas64_.so
=#

I'm running Julia v1.7.3 and LinearSolve v1.22.1

Also, I found that FastLUFactorization is not exported so I had to invoke it as LinearSolve.FastLUFactorization().

[ChrisRackauckas]

How many Julia threads are open on your system? Show versioninfo()

[ChrisRackauckas]

@chriselrod

[rayegun]

LUFact and FastLUFact should have approximately the same time for the initial solve I think? They're both getrf right?

[chriselrod]

Those are complex matrices.
You could make a PR to add actual complex support to RecursiveFactorization.

Or wait for the LoopVectorization rewrite to be finished.
Development is active there, but it will still be a while.

[ChrisRackauckas]

Oh RecursiveFactorization.jl doesn't have a specialization for complex matrices? I didn't know that. We should specialize the default algorithm on that fact and document it then. What exactly is the issue, that it's only fast for Float32 and Float64?

[simonp0420]

8 threads on my system:

julia> versioninfo()
Julia Version 1.7.3
Commit 742b9abb4d (2022-05-06 12:58 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-9700 CPU @ 3.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 8

[ChrisRackauckas]

LUFact and FastLUFact should have approximately the same time for the initial solve I think? They're both getrf right?

Yes, they should. I mean, FastLUFact is just an alternative interface to BLAS that skips rebuilding the workspace each time, so it's not that it should be as fast as BLAS, it is a BLAS call 😅 . So that's a pretty confusing result here that needs to get explored a bit more.

I opened a PR that tweaks the default algorithms a bit more: #160 . Though I think we should figure out what's going on with FastLUFactorization.

[chriselrod]

Oh RecursiveFactorization.jl doesn't have a specialization for complex matrices? I didn't know that. We should specialize the default algorithm on that fact and document it then. What exactly is the issue, that it's only fast for Float32 and Float64?

LoopVectorization.jl currently only supports Union{Bool,Base.HWReal}.

If we have a special need for it somewhere, I could take the time to hack on support to LoopVectorization.jl, otherwise it's better to wait for the rewrite, where it will come for free.

[simonp0420]

I'm happy to wait for the rewrite. However, I'm also not seeing the advertised speedup on Float64:

Algorithm	n = 200	n = 500	n = 2000
Backslash	0.2	1.3	29
(Default)	0.2	1.2	29
LUFactorization	0.2	1.3	29
RFLUFactorization	0.2	1.2	48
Fast LUFactorization	0.3	2.1	53.7

n is the order of the matrix and the table entries are times in msec reported by @btime.

[ChrisRackauckas]

&& <= length(500)

https://github.com/SciML/LinearSolve.jl/blob/main/src/default.jl#L17

[simonp0420]

So the default algorithm was selected to be RFLUFactorization for the n = 200 and n=500 cases and it shouldn't be expected to be competitive for the n = 2000 case? Sounds good. But why isn't it any faster than backslash in the two smaller cases?

[rayegun]

I need to investigate what's happening in FastLUFactorization. This should be taking exactly as long as Backslash, until you do a re-solve where it should be faster. Should be easy to find. I wonder if this is also true for Float64.