Use commonsolve interface

Author: ranjanan

Project.toml

@@ -3,12 +3,14 @@ uuid = "2169fc97-5a83-5252-b627-83903c6c433c"
 version = "0.4.2"
 
 [deps]
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 CompatHelper = "aa819f21-2bde-4658-8897-bab36330d9b7"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

src/AlgebraicMultigrid.jl

@@ -1,9 +1,12 @@
 module AlgebraicMultigrid
 
+using Reexport
 import IterativeSolvers: gauss_seidel!
 using LinearAlgebra
 using SparseArrays, Printf
 using Base.Threads
+@reexport import CommonSolve: solve, solve!, init
+using Reexport
 
 using LinearAlgebra: rmul!
 
@@ -29,7 +32,7 @@ export GaussSeidel, SymmetricSweep, ForwardSweep, BackwardSweep,
         JacobiProlongation
 
 include("multilevel.jl")
-export solve
+export RugeStubenAMG, SmoothedAggregationAMG
 
 include("classical.jl")
 export ruge_stuben

src/multilevel.jl

@@ -126,7 +126,7 @@ struct F <: Cycle
 end
 
 """
-    solve(ml::MultiLevel, b::AbstractArray, cycle, kwargs...)
+    _solve(ml::MultiLevel, b::AbstractArray, cycle, kwargs...)
 
 Execute multigrid cycling.
 
@@ -145,13 +145,13 @@ Keyword Arguments
 * log::Bool - return vector of residuals along with solution
 
 """
-function solve(ml::MultiLevel, b::AbstractArray, args...; kwargs...)
+function _solve(ml::MultiLevel, b::AbstractArray, args...; kwargs...)
     n = length(ml) == 1 ? size(ml.final_A, 1) : size(ml.levels[1].A, 1)
     V = promote_type(eltype(ml.workspace), eltype(b))
     x = zeros(V, size(b))
-    return solve!(x, ml, b, args...; kwargs...)
+    return _solve!(x, ml, b, args...; kwargs...)
 end
-function solve!(x, ml::MultiLevel, b::AbstractArray{T},
+function _solve!(x, ml::MultiLevel, b::AbstractArray{T},
                                     cycle::Cycle = V();
                                     maxiter::Int = 100,
                                     abstol::Real = zero(real(eltype(b))),
@@ -233,3 +233,29 @@ function __solve!(x, ml, cycle::Cycle, b, lvl)
 
     x
 end
+
+
+### CommonSolve.jl spec
+struct AMGSolver{T}
+    ml::MultiLevel
+    b::Vector{T}
+end
+
+abstract type AMGAlg end
+
+struct RugeStubenAMG  <: AMGAlg end
+struct SmoothedAggregationAMG  <: AMGAlg end
+
+function solve(A::AbstractMatrix, b::Vector, s::AMGAlg, args...; kwargs...)
+    solt = init(s, A, b, args...; kwargs...)
+    solve!(solt, args...; kwargs...)
+end
+function init(::RugeStubenAMG, A, b, args...; kwargs...)
+    AMGSolver(ruge_stuben(A; kwargs...), b)
+end
+function init(::SmoothedAggregationAMG, A, b; kwargs...)
+    AMGSolver(smoothed_aggregation(A; kwargs...), b)
+end
+function solve!(solt::AMGSolver, args...; kwargs...) 
+    _solve(solt.ml, solt.b, args...; kwargs...)   
+end

test/cycle_tests.jl

@@ -14,7 +14,7 @@ function test_cycles()
         ml = method(A)
 
         for cycle in [AlgebraicMultigrid.V(),AlgebraicMultigrid.W(),AlgebraicMultigrid.F()]
-            x,convhist = solve(ml, b, cycle; reltol = reltol, log = true)
+            x,convhist = AlgebraicMultigrid._solve(ml, b, cycle; reltol = reltol, log = true)
 
             @debug "number of iterations for $cycle using $method: $(length(convhist))"
             @test norm(b - A*x) < reltol * norm(b)

test/runtests.jl

@@ -327,4 +327,4 @@ end
 X = poisson(27_000)+24.0*I
 ml = ruge_stuben(X)
 b = rand(27_000)
-@test solve(ml, b, reltol = 1e-10) ≈ X \ b rtol = 1e-10
+@test AlgebraicMultigrid._solve(ml, b, reltol = 1e-10) ≈ X \ b rtol = 1e-10