Using KLU per default

Project.toml

@@ -6,6 +6,7 @@ version = "0.11.0"
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec"
@@ -30,6 +31,7 @@ Adapt = "3.3"
 CUDA = "4.1, 5"
 ForwardDiff = "0.10"
 KernelAbstractions = "0.9"
+KLU = "0.4"
 Krylov = "0.9"
 KrylovPreconditioners = "0.2"
 LazyArtifacts = "1.9"

examples/power_flow.jl

@@ -54,8 +54,8 @@ precond = LS.BlockJacobiPreconditioner(J, npartitions, localdevice, noverlap)
 #=
     Instantiate iterative linear solver
 =#
-lin_solver = ExaPF.LinearSolvers.KrylovBICGSTAB
-algo = LS.KrylovBICGSTAB(J; P=precond)
+lin_solver = ExaPF.LinearSolvers.Bicgstab
+algo = LS.Bicgstab(J; P=precond)
 
 #=
     Power flow resolution

ext/ExaPFAMDGPUExt.jl

@@ -21,10 +21,10 @@ const KP = KrylovPreconditioners
 LS.DirectSolver(J::ROCSparseMatrixCSR; options...) = ExaPF.LS.DirectSolver(nothing)
 LS.update!(solver::ExaPF.LS.AbstractIterativeLinearSolver, J::ROCSparseMatrixCSR) = KP.update!(solver.precond, J)
 LS._get_type(J::ROCSparseMatrixCSR) = ROCArray{Float64, 1, AMDGPU.Mem.HIPBuffer}
-LS.default_linear_solver(A::ROCSparseMatrixCSR, device::ROCBackend) = ExaPF.LS.KrylovBICGSTAB(A)
+LS.default_linear_solver(A::ROCSparseMatrixCSR, device::ROCBackend) = ExaPF.LS.Bicgstab(A)
 ExaPF._iscsr(::ROCSparseMatrixCSR) = true
 ExaPF._iscsc(::ROCSparseMatrixCSR) = false
-function LS.scaling!(::LS.KrylovBICGSTAB,A::ROCSparseMatrixCSR,b)
+function LS.scaling!(::LS.Bicgstab, A::ROCSparseMatrixCSR, b)
     KP.scaling_csr!(A,b)
 end
 
@@ -33,7 +33,7 @@ end
 
 List all linear solvers available solving the power flow on an NVIDIA GPU.
 """
-ExaPF.list_solvers(::ROCBackend) = [LS.BICGSTAB, LS.DQGMRES, LS.EigenBICGSTAB, LS.KrylovBICGSTAB]
+ExaPF.list_solvers(::ROCBackend) = [LS.Dqgmres, LS.Bicgstab]
 
 include("amdgpu_wrapper.jl")
 end

ext/ExaPFCUDAExt.jl

@@ -24,15 +24,15 @@ LS._get_type(J::CuSparseMatrixCSR) = CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}
 LS.default_linear_solver(A::CuSparseMatrixCSR, device::CUDABackend) = ExaPF.LS.DirectSolver(A)
 ExaPF._iscsr(::CuSparseMatrixCSR) = true
 ExaPF._iscsc(::CuSparseMatrixCSR) = false
-function LS.scaling!(::LS.KrylovBICGSTAB,A::CuSparseMatrixCSR,b)
+function LS.scaling!(::LS.Bicgstab, A::CuSparseMatrixCSR, b)
     KP.scaling_csr!(A,b)
 end
 """
     list_solvers(::CUDABackend)
 
 List all linear solvers available solving the power flow on an NVIDIA GPU.
 """
-ExaPF.list_solvers(::CUDABackend) = [LS.DirectSolver, LS.BICGSTAB, LS.DQGMRES, LS.EigenBICGSTAB, LS.KrylovBICGSTAB]
+ExaPF.list_solvers(::CUDABackend) = [LS.DirectSolver, LS.Dqgmres, LS.Bicgstab]
 
 include("cuda_wrapper.jl")
 end

scripts/linear_solvers/benchmark.jl

@@ -65,7 +65,7 @@ function benchmark_gpu_krylov(datafile, pf_solver; ntrials=3)
     n_partitions = div(n_states, n_blocks)
     jac_gpu = instance.jacobian.J
     precond = BlockJacobiPreconditioner(jac_gpu, n_partitions, CUDABackend(), 0)
-    krylov_solver = ExaPF.KrylovBICGSTAB(
+    krylov_solver = ExaPF.Bicgstab(
         jac_gpu; P=precond, ldiv=false, scaling=true,
         rtol=1e-7, atol=1e-7, verbose=0,
     )

scripts/linear_solvers/benchmark_block.jl

@@ -85,7 +85,7 @@ function benchmark_gpu_krylov(datafile, n_scenarios, pf_solver; ntrials=3, magni
     n_partitions = div(n_states, n_blocks)
     jac_gpu = instance.jacobian.J
     precond = BlockJacobiPreconditioner(jac_gpu, n_partitions, CUDABackend(), 0)
-    krylov_solver = ExaPF.KrylovBICGSTAB(
+    krylov_solver = ExaPF.Bicgstab(
         jac_gpu; P=precond, ldiv=false, scaling=true,
         rtol=1e-7, atol=1e-7, verbose=0,
     )

src/LinearSolvers/LinearSolvers.jl

@@ -16,14 +16,16 @@ import Metis
 import ..ExaPF: xnorm
 
 import Base.size, Base.sizeof, Base.format_bytes
+
+using KLU
 import Krylov: KrylovStats, allocate_if, ksizeof, FloatOrComplex, ktimer, matrix_to_vector, kdisplay, mulorldiv!
 using KrylovPreconditioners
 
 const KA = KernelAbstractions
 const KP = KrylovPreconditioners
 
 export bicgstab, list_solvers, default_linear_solver
-export DirectSolver, BICGSTAB, EigenBICGSTAB, KrylovBICGSTAB
+export DirectSolver, BICGSTAB, EigenBICGSTAB, Bicgstab
 export do_scaling, scaling!
 
 @enum(
@@ -35,9 +37,6 @@ export do_scaling, scaling!
     Diverged,
 )
 
-include("bicgstab.jl")
-include("bicgstab_eigen.jl")
-
 include("utils.jl")
 include("block_gmres.jl")
 
@@ -102,11 +101,11 @@ struct DirectSolver{Fac<:Union{Nothing, LinearAlgebra.Factorization}} <: Abstrac
     factorization::Fac
 end
 
-DirectSolver(J; options...) = DirectSolver(lu(J))
+DirectSolver(J; options...) = DirectSolver(klu(J))
 DirectSolver() = DirectSolver(nothing)
 
 function update!(s::DirectSolver, J::AbstractMatrix)
-    lu!(s.factorization, J) # Update factorization inplace
+    klu!(s.factorization, J) # Update factorization inplace
 end
 
 # Reuse factorization in update
@@ -138,89 +137,28 @@ end
 update!(solver::AbstractIterativeLinearSolver, J::SparseMatrixCSC) = KP.update!(solver.precond, J)
 
 """
-    BICGSTAB <: AbstractIterativeLinearSolver
-    BICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)
-
-Custom BICGSTAB implementation to solve iteratively the linear system
-``A  x = y``.
-"""
-struct BICGSTAB <: AbstractIterativeLinearSolver
-    precond::AbstractKrylovPreconditioner
-    maxiter::Int
-    tol::Float64
-    verbose::Bool
-end
-function BICGSTAB(J::AbstractSparseMatrix;
-    P=BlockJacobiPreconditioner(J), maxiter=2_000, tol=1e-8, verbose=false
-)
-    return BICGSTAB(P, maxiter, tol, verbose)
-end
-
-function ldiv!(solver::BICGSTAB,
-    y::AbstractVector, J::AbstractMatrix, x::AbstractVector; options...
-)
-    y[:], n_iters, status = bicgstab(J, x, solver.precond, y; maxiter=solver.maxiter,
-                                     verbose=solver.verbose, tol=solver.tol)
-    if status != Converged
-        @warn("BICGSTAB failed to converge. Final status is $(status)")
-    end
-    return n_iters
-end
-
-"""
-    EigenBICGSTAB <: AbstractIterativeLinearSolver
-    EigenBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)
-
-Julia's port of Eigen's BICGSTAB to solve iteratively the linear system
-``A x = y``.
-"""
-struct EigenBICGSTAB <: AbstractIterativeLinearSolver
-    precond::AbstractKrylovPreconditioner
-    maxiter::Int
-    tol::Float64
-    verbose::Bool
-end
-function EigenBICGSTAB(J::AbstractSparseMatrix;
-    P=BlockJacobiPreconditioner(J), maxiter=2_000, tol=1e-8, verbose=false
-)
-    return EigenBICGSTAB(P, maxiter, tol, verbose)
-end
-
-function ldiv!(solver::EigenBICGSTAB,
-    y::AbstractVector, J::AbstractMatrix, x::AbstractVector; options...
-)
-    y[:], n_iters, status = bicgstab_eigen(J, x, solver.precond, y; maxiter=solver.maxiter,
-                                           verbose=solver.verbose, tol=solver.tol)
-    if status != Converged
-        @warn("EigenBICGSTAB failed to converge. Final status is $(status)")
-    end
-
-    return n_iters
-end
-
-"""
-    DQGMRES <: AbstractIterativeLinearSolver
-    DQGMRES(precond; verbose=false, memory=4)
+    Dqgmres <: AbstractIterativeLinearSolver
+    Dqgmres(precond; verbose=false, memory=4)
 
-Wrap `Krylov.jl`'s DQGMRES algorithm to solve iteratively the linear system
+Wrap `Krylov.jl`'s Dqgmres algorithm to solve iteratively the linear system
 ``A x = y``.
 """
-struct DQGMRES <: AbstractIterativeLinearSolver
+struct Dqgmres <: AbstractIterativeLinearSolver
     inner::Krylov.DqgmresSolver
     precond::AbstractKrylovPreconditioner
     memory::Int
     verbose::Bool
 end
-function DQGMRES(J::AbstractSparseMatrix;
+function Dqgmres(J::AbstractSparseMatrix;
     P=BlockJacobiPreconditioner(J), memory=4, verbose=false
 )
     n, m = size(J)
     S = _get_type(J)
     solver = Krylov.DqgmresSolver(n, m, memory, S)
-    return DQGMRES(solver, P, memory, verbose)
+    return Dqgmres(solver, P, memory, verbose)
 end
 
-function ldiv!(solver::DQGMRES,
+function ldiv!(solver::Dqgmres,
     y::AbstractVector, J::AbstractMatrix, x::AbstractVector; options...
 )
     Krylov.dqgmres!(solver.inner, J, x; N=solver.precond)
@@ -229,13 +167,13 @@ function ldiv!(solver::DQGMRES,
 end
 
 """
-    KrylovBICGSTAB <: AbstractIterativeLinearSolver
-    KrylovBICGSTAB(precond; verbose=0, rtol=1e-10, atol=1e-10)
+    Bicgstab <: AbstractIterativeLinearSolver
+    Bicgstab(precond; verbose=0, rtol=1e-10, atol=1e-10)
 
 Wrap `Krylov.jl`'s BICGSTAB algorithm to solve iteratively the linear system
 ``A x = y``.
 """
-struct KrylovBICGSTAB <: AbstractIterativeLinearSolver
+struct Bicgstab <: AbstractIterativeLinearSolver
     inner::Krylov.BicgstabSolver
     precond::AbstractKrylovPreconditioner
     verbose::Int
@@ -245,17 +183,17 @@ struct KrylovBICGSTAB <: AbstractIterativeLinearSolver
     scaling::Bool
     maxiter::Int64
 end
-do_scaling(linear_solver::KrylovBICGSTAB) = linear_solver.scaling
-function KrylovBICGSTAB(J::AbstractSparseMatrix;
+do_scaling(linear_solver::Bicgstab) = linear_solver.scaling
+function Bicgstab(J::AbstractSparseMatrix;
     P=BlockJacobiPreconditioner(J), verbose=0, rtol=1e-10, atol=1e-10, ldiv=false, scaling=false, maxiter=size(J,1)
 )
     n, m = size(J)
     S = _get_type(J)
     solver = Krylov.BicgstabSolver(n, m, S)
-    return KrylovBICGSTAB(solver, P, verbose, atol, rtol, ldiv, scaling, maxiter)
+    return Bicgstab(solver, P, verbose, atol, rtol, ldiv, scaling, maxiter)
 end
 
-function ldiv!(solver::KrylovBICGSTAB,
+function ldiv!(solver::Bicgstab,
     y::AbstractVector, J::AbstractMatrix, x::AbstractVector;
     max_atol = solver.atol, max_rtol = solver.rtol, options...
 )
@@ -284,7 +222,7 @@ end
 
 List all linear solvers available solving the power flow on the CPU.
 """
-list_solvers(::KA.CPU) = [DirectSolver, DQGMRES, BICGSTAB, EigenBICGSTAB, KrylovBICGSTAB]
+list_solvers(::KA.CPU) = [DirectSolver, Dqgmres, Bicgstab]
 
 """
     default_linear_solver(A, ::KA.CPU)