Skip to content

Commit ae0b5e0

Browse files
format
1 parent 3beafc8 commit ae0b5e0

File tree

4 files changed

+51
-44
lines changed

4 files changed

+51
-44
lines changed

docs/pages.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@
22

33
pages = ["index.md",
44
"Tutorials" => Any["tutorials/linear.md",
5-
"tutorials/caching_interface.md",
6-
"tutorials/accelerating_choices.md"],
5+
"tutorials/caching_interface.md",
6+
"tutorials/accelerating_choices.md"],
77
"Basics" => Any["basics/LinearProblem.md",
88
"basics/common_solver_opts.md",
99
"basics/OperatorAssumptions.md",

docs/src/solvers/solvers.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -129,7 +129,7 @@ LinearSolve.jl provides a wrapper to these routines in a way where an initialize
129129
has a non-allocating LU factorization. In theory, this post-initialized solve should always
130130
be faster than the Base.LinearAlgebra version. In practice, with the way we wrap the solvers,
131131
we do not see a performance benefit and in fact benchmarks tend to show this inhibits
132-
performance.
132+
performance.
133133

134134
!!! note
135135

Lines changed: 36 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
# Accelerating your Linear Solves
22

33
!!! note
4+
45
This section is essential if you wish to achieve maximum performance with
56
LinearSolve.jl, especially on v7 and above. Please ensure the tips of this
67
section are adhered to when optimizing code and benchmarking.
@@ -16,7 +17,7 @@ scenarios, so let's dive in.
1617
## Understanding Performance of Dense Linear Solves
1718

1819
The performance of dense linear solvers is highly dependent on the size of the matrix
19-
and the chosen architecture to run on, i.e. the CPU.
20+
and the chosen architecture to run on, i.e. the CPU.
2021
[This issue](https://github.com/SciML/LinearSolve.jl/issues/357) gathered benchmark data
2122
from many different users and is summarized in the following graphs:
2223

@@ -25,33 +26,34 @@ from many different users and is summarized in the following graphs:
2526
Now one thing that is immediate is for example that AppleAccelerate generally does well
2627
on Apple M-series chips, MKL generally does well on Intel, etc. And we know this in
2728
LinearSolve.jl, in fact we automatically default to different BLASes based on the CPU
28-
architecture already as part of the design! So that covers most of the variation, but
29+
architecture already as part of the design! So that covers most of the variation, but
2930
there are a few major tips to note when fine tuning the results to your system:
3031

31-
1. One of the best methods for size 150x150 matrices and below is RecursiveFactorization.jl.
32-
This is a pure Julia BLAS system, but it has a high load time overhead, and thus as of
33-
v7 it's no longer loaded by default! Thus if your matrices are in this range and you would
34-
value better run times at the cost of compile and load times, it is recommended you add
35-
`using RecursiveFactorization`. The defaulting algorithm will then consider it in its list
36-
and will automatically (in an architecture-specific way) insert it as it feels necesssary.
37-
2. One of the major factors that can inhibit BLAS performance on LU factorization is multithreading.
38-
In many of these plots you can see a giant dip in GFLOPs (higher is better) when a certain size
39-
threshold is hit. This is because, for the number of chosen threads, there was not enough work
40-
and thus when the threading threshold is hit you get a hit to the performance due to the added
41-
overhead. The threading performance can be a per-system thing, and it can be greatly influenced
42-
by the number of cores on your system and the number of threads you allow. Thus for example,
43-
OpenBLAS' LU factorization seems to generally be really bad at guessing the thread switch point
44-
for CPUs with really high core/thread counts. If this is the case, you may want to investigate
45-
decreasing your number of BLAS threads, i.e. via `BLAS.set_num_threads(i)`. Note that
46-
RecursiveFactorization.jl uses your Julia thread pool instead of the BLAS threads.
47-
3. The switch points between algorithms can be fairly inexact. LinearSolve.jl tried to keep a tab
48-
on where they are per platform and keep updated, but it can be a moving battle. You may be
49-
able to eek out some performance by testing between the various options on your platform, i.e.
50-
RFLUFactorization vs LUFactorization vs AppleAccelerateLUFactorization (M-series) vs
51-
MKLFactorization (X86) and hardcoding the choice for your problem if the default did not make
52-
the right guess.
32+
1. One of the best methods for size 150x150 matrices and below is RecursiveFactorization.jl.
33+
This is a pure Julia BLAS system, but it has a high load time overhead, and thus as of
34+
v7 it's no longer loaded by default! Thus if your matrices are in this range and you would
35+
value better run times at the cost of compile and load times, it is recommended you add
36+
`using RecursiveFactorization`. The defaulting algorithm will then consider it in its list
37+
and will automatically (in an architecture-specific way) insert it as it feels necesssary.
38+
2. One of the major factors that can inhibit BLAS performance on LU factorization is multithreading.
39+
In many of these plots you can see a giant dip in GFLOPs (higher is better) when a certain size
40+
threshold is hit. This is because, for the number of chosen threads, there was not enough work
41+
and thus when the threading threshold is hit you get a hit to the performance due to the added
42+
overhead. The threading performance can be a per-system thing, and it can be greatly influenced
43+
by the number of cores on your system and the number of threads you allow. Thus for example,
44+
OpenBLAS' LU factorization seems to generally be really bad at guessing the thread switch point
45+
for CPUs with really high core/thread counts. If this is the case, you may want to investigate
46+
decreasing your number of BLAS threads, i.e. via `BLAS.set_num_threads(i)`. Note that
47+
RecursiveFactorization.jl uses your Julia thread pool instead of the BLAS threads.
48+
3. The switch points between algorithms can be fairly inexact. LinearSolve.jl tried to keep a tab
49+
on where they are per platform and keep updated, but it can be a moving battle. You may be
50+
able to eek out some performance by testing between the various options on your platform, i.e.
51+
RFLUFactorization vs LUFactorization vs AppleAccelerateLUFactorization (M-series) vs
52+
MKLFactorization (X86) and hardcoding the choice for your problem if the default did not make
53+
the right guess.
5354

5455
!!! warn
56+
5557
As noted, RecursiveFactorization.jl is one of the fastest linear solvers for smaller dense
5658
matrices but requires `using RecursiveFactorization` in order to be used in the default
5759
solver setups! Thus it's recommended that any optimized code or benchmarks sets this up.
@@ -69,17 +71,18 @@ LinearSolve.jl just uses a very simple "if small then use KLU and if large use U
6971
is validated by this plot, but leaves a lot to be desired. In particular, the following rules
7072
should be thought about:
7173

72-
1. Pardiso is a great solver, you should try `using Pardiso` and using `MKLPardiso()` in many
73-
scenarios.
74-
2. The more structured a sparsity pattern is, the worse KLU is in comparison to the other
75-
algorithms.
76-
3. A Krylov subspace method with proper preconditioning will be better than direct solvers
77-
when the matrices get large enough. You could always precondition a sparse matrix with
78-
iLU as an easy choice, though the tolerance would need to be tuned in a problem-specific
79-
way.
74+
1. Pardiso is a great solver, you should try `using Pardiso` and using `MKLPardiso()` in many
75+
scenarios.
76+
2. The more structured a sparsity pattern is, the worse KLU is in comparison to the other
77+
algorithms.
78+
3. A Krylov subspace method with proper preconditioning will be better than direct solvers
79+
when the matrices get large enough. You could always precondition a sparse matrix with
80+
iLU as an easy choice, though the tolerance would need to be tuned in a problem-specific
81+
way.
8082

8183
!!! note
84+
8285
UMFPACK does better when the BLAS is not OpenBLAS. Try `using MKL` on Intel and AMD Ryzen
8386
platforms and UMPACK will be faster! LinearSolve.jl cannot default to this as this changes
8487
global settings and thus only defaults to MKL locally, and thus cannot change the setting
85-
within UMFPACK.
88+
within UMFPACK.

ext/LinearSolveFastLapackInterfaceExt.jl

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -12,17 +12,20 @@ function LinearSolve.init_cacheval(::FastLUFactorization, A, b, u, Pl, Pr,
1212
maxiters::Int, abstol, reltol, verbose::Bool,
1313
assumptions::OperatorAssumptions)
1414
ws = LUWs(A)
15-
return WorkspaceAndFactors(ws, LinearSolve.ArrayInterface.lu_instance(convert(AbstractMatrix, A)))
15+
return WorkspaceAndFactors(
16+
ws, LinearSolve.ArrayInterface.lu_instance(convert(AbstractMatrix, A)))
1617
end
1718

18-
function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::FastLUFactorization; kwargs...)
19+
function SciMLBase.solve!(
20+
cache::LinearSolve.LinearCache, alg::FastLUFactorization; kwargs...)
1921
A = cache.A
2022
A = convert(AbstractMatrix, A)
2123
ws_and_fact = LinearSolve.@get_cacheval(cache, :FastLUFactorization)
2224
if cache.isfresh
2325
# we will fail here if A is a different *size* than in a previous version of the same cache.
2426
# it may instead be desirable to resize the workspace.
25-
LinearSolve.@set! ws_and_fact.factors = LinearAlgebra.LU(LAPACK.getrf!(ws_and_fact.workspace,
27+
LinearSolve.@set! ws_and_fact.factors = LinearAlgebra.LU(LAPACK.getrf!(
28+
ws_and_fact.workspace,
2629
A)...)
2730
cache.cacheval = ws_and_fact
2831
cache.isfresh = false
@@ -31,19 +34,21 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::FastLUFactorizati
3134
SciMLBase.build_linear_solution(alg, y, nothing, cache)
3235
end
3336

34-
function LinearSolve.init_cacheval(alg::FastQRFactorization{NoPivot}, A::AbstractMatrix, b, u, Pl, Pr,
37+
function LinearSolve.init_cacheval(
38+
alg::FastQRFactorization{NoPivot}, A::AbstractMatrix, b, u, Pl, Pr,
3539
maxiters::Int, abstol, reltol, verbose::Bool,
3640
assumptions::OperatorAssumptions)
3741
ws = QRWYWs(A; blocksize = alg.blocksize)
3842
return WorkspaceAndFactors(ws,
3943
LinearSolve.ArrayInterface.qr_instance(convert(AbstractMatrix, A)))
4044
end
41-
function LinearSolve.init_cacheval(::FastQRFactorization{ColumnNorm}, A::AbstractMatrix, b, u, Pl, Pr,
45+
function LinearSolve.init_cacheval(
46+
::FastQRFactorization{ColumnNorm}, A::AbstractMatrix, b, u, Pl, Pr,
4247
maxiters::Int, abstol, reltol, verbose::Bool,
4348
assumptions::OperatorAssumptions)
4449
ws = QRpWs(A)
4550
return WorkspaceAndFactors(ws,
46-
LinearSolve.ArrayInterface.qr_instance(convert(AbstractMatrix, A)))
51+
LinearSolve.ArrayInterface.qr_instance(convert(AbstractMatrix, A)))
4752
end
4853

4954
function LinearSolve.init_cacheval(alg::FastQRFactorization, A, b, u, Pl, Pr,
@@ -78,5 +83,4 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::FastQRFactorizati
7883
SciMLBase.build_linear_solution(alg, y, nothing, cache)
7984
end
8085

81-
82-
end
86+
end

0 commit comments

Comments
 (0)