Solving Sparse Linear Systems fast

General Usage Performance
,
9

You don’t need it for MKLPardiso.

MWE with LinearSolve.jl:

using LinearSolve, LinearSolvePardiso, SparseArrays, Random

A = sparse([1.0 0 -2 3
    0 5 1 2
    -2 1 4 -7
    3 2 -7 5])
b = rand(4)
prob = LinearProblem(A, b)
u = solve(prob, MKLPardisoFactorize(); cache_kwargs...).u

should work and automatically install any binaries needed.

This has a regular sparsity pattern, so it should be a case where UMFPACK-type algorithms do well because they should be able to exploit the parallelism of the BLAS3 kernel. I’d try:

using LinearSolve, LinearSolvePardiso, SparseArrays
A = # your matrix
b = # test b

prob = LinearProblem(A, b)

for alg in (
    MKLPardisoFactorize(),
    MKLPardisoIterate(),
    UMFPACKFactorization(),
    KLUFactorization())

    u = solve(prob, alg).u
end

using AlgebraicMultigrid
ml = ruge_stuben(A) # Construct a Ruge-Stuben solver
pl = aspreconditioner(ml)
solve(prob1, KrylovJL_GMRES(), Pl = pl).u

using IncompleteLU
pl = ilu(A, τ = 0.1) # τ needs to be tuned per problem
solve(prob1, KrylovJL_GMRES(), Pl = pl).u

I should make that into a tutorial. You can see this in action in https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/

6 Likes