Merge pull request #87 from JuliaLinearAlgebra/RA/commonsolve

Add 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"
 

README.md

@@ -9,6 +9,19 @@ This package lets you solve sparse linear systems using Algebraic Multigrid (AMG
 
 ## Usage
 
+### Using the CommonSolve interface
+
+This is highest level API. It internally creates the multilevel object 
+and calls the multigrid cycling `_solve`. 
+
+```julia
+A = poisson(100); 
+b = rand(100);
+solve(A, b, RugeStubenAMG(), maxiter = 1, abstol = 1e-6)
+```
+
+### Multigrid cycling
+
 ```julia
 using AlgebraicMultigrid
 
@@ -32,7 +45,7 @@ ml = ruge_stuben(A) # Construct a Ruge-Stuben solver
 #     8            7           19 [ 0.32%]
 
 
-solve(ml, A * ones(1000)) # should return ones(1000)
+AlgebraicMultigrid._solve(ml, A * ones(1000)) # should return ones(1000)
 ```
 
 ### As a Preconditioner

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/aggregation.jl

@@ -11,7 +11,7 @@ function smoothed_aggregation(A::TA,
                         max_coarse = 10,
                         diagonal_dominance = false,
                         keep = false,
-                        coarse_solver = Pinv) where {T,V,bs,TA<:SparseMatrixCSC{T,V}}
+                        coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}}
 
     n = size(A, 1)
     # B = kron(ones(n, 1), eye(1))

src/classical.jl

@@ -15,7 +15,7 @@ function ruge_stuben(_A::Union{TA, Symmetric{Ti, TA}, Hermitian{Ti, TA}},
                 postsmoother = GaussSeidel(),
                 max_levels = 10,
                 max_coarse = 10,
-                coarse_solver = Pinv) where {Ti,Tv,bs,TA<:SparseMatrixCSC{Ti,Tv}}
+                coarse_solver = Pinv, kwargs...) where {Ti,Tv,bs,TA<:SparseMatrixCSC{Ti,Tv}}
 
     s = Solver(strength, CF, presmoother,
                 postsmoother, max_levels, max_levels)

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,20 +145,20 @@ 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))),
                                     reltol::Real = sqrt(eps(real(eltype(b)))),
                                     verbose::Bool = false,
                                     log::Bool = false,
-                                    calculate_residual = true) where {T}
+                                    calculate_residual = true, kwargs...) where {T}
 
     A = length(ml) == 1 ? ml.final_A : ml.levels[1].A
     V = promote_type(eltype(A), eltype(b))
@@ -233,3 +233,28 @@ 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

src/preconditioner.jl

@@ -15,7 +15,7 @@ function ldiv!(x, p::Preconditioner, b)
     else
         x .= b
     end
-    solve!(x, p.ml, b, p.cycle, maxiter = 1, calculate_residual = false)
+    _solve!(x, p.ml, b, p.cycle, maxiter = 1, calculate_residual = false)
 end
 mul!(b, p::Preconditioner, x) = mul!(b, p.ml.levels[1].A, x)
 

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

@@ -120,29 +120,26 @@ fsmoother = GaussSeidel(ForwardSweep())
 A = poisson(1000)
 A = float.(A)
 ml = ruge_stuben(A)
-x = solve(ml, A * ones(1000))
+x = AlgebraicMultigrid._solve(ml, A * ones(1000))
 @test sum(abs2, x - ones(1000)) < 1e-8
 
 ml = ruge_stuben(A, presmoother = fsmoother,
                      postsmoother = fsmoother)
-x = solve(ml, A * ones(1000))
+x = AlgebraicMultigrid._solve(ml, A * ones(1000))
 @test sum(abs2, x - ones(1000)) < 1e-8
 
 
 A = include("randlap.jl")
 
 ml = ruge_stuben(A, presmoother = fsmoother,
                     postsmoother = fsmoother)
-x = solve(ml, A * ones(100))
+x = AlgebraicMultigrid._solve(ml, A * ones(100))
 @test sum(abs2, x - zeros(100)) < 1e-8
 
 ml = ruge_stuben(A)
-x = solve(ml, A * ones(100))
+x = AlgebraicMultigrid._solve(ml, A * ones(100))
 @test sum(abs2, x - zeros(100)) < 1e-6
 
-
-
-
 end
 
 @testset "Preconditioning" begin
@@ -155,7 +152,7 @@ p = aspreconditioner(ml)
 b = zeros(n)
 b[1] = 1
 b[2] = -1
-x = solve(p.ml, A * ones(n), maxiter = 1, abstol = 1e-12)
+x = AlgebraicMultigrid._solve(ml, A * ones(n), maxiter = 1, abstol = 1e-12)
 diff = x - [  1.88664780e-16,   2.34982727e-16,   2.33917697e-16,
          8.77869044e-17,   7.16783490e-17,   1.43415460e-16,
          3.69199021e-17,   9.70950385e-17,   4.77034895e-17,
@@ -173,7 +170,9 @@ diff = x - [  1.88664780e-16,   2.34982727e-16,   2.33917697e-16,
         -6.76965535e-16,  -7.00643227e-16,  -6.23581397e-16,
         -7.03016682e-16]
 @test sum(abs2, diff) < 1e-8
-x = solve(p.ml, b, maxiter = 1, abstol = 1e-12)
+x = solve(A, b, RugeStubenAMG(); presmoother = smoother, 
+                                postsmoother = smoother,
+                                maxiter = 1, abstol = 1e-12)
 diff = x - [ 0.76347046, -0.5498286 , -0.2705487 , -0.15047352, -0.10248021,
         0.60292674, -0.11497073, -0.08460548, -0.06931461,  0.38230708,
        -0.055664  , -0.04854558, -0.04577031,  0.09964325,  0.01825624,
@@ -214,7 +213,7 @@ diff = x - [0.823762, -0.537478, -0.306212, -0.19359, -0.147621, 0.685002,
             0.0511691, 0.0502043, 0.0498349, 0.0498134]
 @test sum(abs2, diff) < 1e-8
 
-x = solve(ml, b, maxiter = 1, reltol = 1e-12)
+x = AlgebraicMultigrid._solve(ml, b, maxiter = 1, reltol = 1e-12)
 diff = x - [0.775725, -0.571202, -0.290989, -0.157001, -0.106981, 0.622652,
             -0.122318, -0.0891874, -0.0709834, 0.392621, -0.055544, -0.0507485,
             -0.0466376, 0.107175, 0.0267468, -0.0200843, -0.0282827, -0.0299929,
@@ -240,7 +239,7 @@ for (T,V) in ((Float64, Float64), (Float32,Float32),
         ml = smoothed_aggregation(a)
         b = V.(b)
         c = cg(a, b, maxiter = 10)
-        @test eltype(solve(ml, b)) == eltype(c)
+        @test eltype(AlgebraicMultigrid._solve(ml, b)) == eltype(c)
 end
 
 end
@@ -308,15 +307,16 @@ for sz in [10, 5, 2]
 end
 
 # Issue #46
-for f in (smoothed_aggregation, ruge_stuben)
+for f in ((smoothed_aggregation, SmoothedAggregationAMG), 
+            (ruge_stuben, RugeStubenAMG))
 
     a = load("bug.jld2")["G"]
-    ml = f(a)
+    ml = f[1](a)
     p = aspreconditioner(ml)
     b = zeros(size(a,1))
     b[1] = 1
     b[2] = -1
-    @test sum(abs2, a * solve(ml, b) - b) < 1e-10
+    @test sum(abs2, a * solve(a, b, f[2]()) - b) < 1e-10
     @test sum(abs2, a * cg(a, b, Pl = p, maxiter = 1000) - b) < 1e-10
 
 end
@@ -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