LinearSolve.jl integration and remove some dead code.

Author: termi-official

Project.toml

@@ -5,6 +5,7 @@ version = "0.5.1"
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

src/AlgebraicMultigrid.jl

@@ -2,18 +2,13 @@ module AlgebraicMultigrid
 
 using Reexport
 using LinearAlgebra
+using LinearSolve
 using SparseArrays, Printf
-using Base.Threads
 @reexport import CommonSolve: solve, solve!, init
 using Reexport
 
 using LinearAlgebra: rmul!
 
-# const mul! = A_mul_B!
-
-const MT = false
-const AMG = AlgebraicMultigrid
-
 include("utils.jl")
 export approximate_spectral_radius
 

src/multilevel.jl

@@ -52,6 +52,12 @@ function coarse_b!(m::MultiLevelWorkspace{TX, bs}, n) where {TX, bs}
 end
 
 abstract type CoarseSolver end
+
+"""
+    Pinv{T} <: CoarseSolver
+
+Moore-Penrose pseudo inverse coarse solver. Calls `pinv`
+"""
 struct Pinv{T} <: CoarseSolver
     pinvA::Matrix{T}
     Pinv{T}(A) where T = new{T}(pinv(Matrix(A)))
@@ -61,6 +67,43 @@ Base.show(io::IO, p::Pinv) = print(io, "Pinv")
 
 (p::Pinv)(x, b) = mul!(x, p.pinvA, b)
 
+# This one is used internally.
+"""
+    LinearSolveWrapperInternal <: CoarseSolver
+
+Helper to allow the usage of LinearSolve.jl solvers for the coarse-level solve. Constructed via `LinearSolveWrapper`.
+"""
+struct LinearSolveWrapperInternal{LC <: LinearSolve.LinearCache} <: CoarseSolver
+    linsolve::LC
+    function LinearSolveWrapperInternal(A, alg::LinearSolve.SciMLLinearSolveAlgorithm)
+        rhs_tmp = zeros(eltype(A), size(A,1))
+        u_tmp   = zeros(eltype(A), size(A,2))
+        linprob = LinearProblem(A, rhs_tmp; u0 = u_tmp, alias_A = false, alias_b = false)
+        linsolve = init(linprob, alg)
+        new{typeof(linsolve)}(linsolve)
+    end
+end
+
+function (p::LinearSolveWrapperInternal{LC})(x, b) where {LC <: LinearSolve.LinearCache}
+    for i ∈ 1:size(b, 2)
+        # Update right hand side
+        p.linsolve.b = b[:, i]
+        # Solve for x and update
+        x[:, i] = solve!(p.linsolve).u
+    end
+end
+
+# This one simplifies passing of LinearSolve.jl algorithms into AlgebraicMultigrid.jl as coarse solvers.
+"""
+    LinearSolveWrapper <: CoarseSolver
+
+Helper to allow the usage of LinearSolve.jl solvers for the coarse-level solve.
+"""
+struct LinearSolveWrapper <: CoarseSolver
+    alg::LinearSolve.SciMLLinearSolveAlgorithm
+end
+(p::LinearSolveWrapper)(A::AbstractMatrix) = LinearSolveWrapperInternal(A, p.alg)
+
 Base.length(ml::MultiLevel) = length(ml.levels) + 1
 
 function Base.show(io::IO, ml::MultiLevel)

src/utils.jl

@@ -99,7 +99,6 @@ end
 find_breakdown(::Type{Float64}) = eps(Float64) * 10^6
 find_breakdown(::Type{Float32}) = eps(Float64) * 10^3
 
-using Base.Threads
 #=function mul!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)
     A.n == size(B, 1) || throw(DimensionMismatch())
     A.m == size(C, 1) || throw(DimensionMismatch())