Merge pull request #116 from j-fu/LinearSolvePrecs

Provide interface to the precs API of LinearSolve

Author: ranjanan

README.md

@@ -56,6 +56,26 @@ p = aspreconditioner(ml)
 c = cg(A, A*ones(1000), Pl = p)
 ```
 
+
+### As a preconditioner with LinearSolve.jl
+
+`RugeStubenPreconBuilder` and `SmoothedAggregationPreconBuilder` work with the 
+[`precs` API](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#Specifying-Preconditioners)
+of LinearSolve.jl
+
+```julia
+A = poisson( (100,100) )
+u0= rand(size(A,1))
+b=A*u0
+
+prob = LinearProblem(A, b)
+strategy = KrylovJL_CG(precs = RugeStubenPreconBuilder())
+sol = solve(prob, strategy, atol=1.0e-14)
+
+strategy = KrylovJL_CG(precs = SmoothedAggregationPreconBuilder())
+sol = solve(prob, strategy, atol=1.0e-14)
+```
+
 ## Features and Roadmap
 
 This package currently supports:

src/AlgebraicMultigrid.jl

@@ -3,7 +3,9 @@ module AlgebraicMultigrid
 using Reexport
 using LinearAlgebra
 using LinearSolve
-using SparseArrays, Printf
+using SparseArrays
+using SparseArrays: AbstractSparseMatrixCSC
+using Printf
 @reexport import CommonSolve: solve, solve!, init
 using Reexport
 
@@ -40,4 +42,7 @@ export fit_candidates, smoothed_aggregation
 include("preconditioner.jl")
 export aspreconditioner
 
+include("precs.jl")
+export SmoothedAggregationPreconBuilder, RugeStubenPreconBuilder
+
 end # module

src/precs.jl

@@ -0,0 +1,38 @@
+"""
+    SmoothedAggregationPreconBuilder(;blocksize=1, kwargs...)
+
+Return callable object constructing a left smoothed aggregation algebraic multigrid preconditioner
+to be used with the `precs` API of LinearSolve.
+"""
+struct SmoothedAggregationPreconBuilder{Tk}
+    blocksize::Int
+    kwargs::Tk
+end
+
+function SmoothedAggregationPreconBuilder(; blocksize = 1, kwargs...)
+    return SmoothedAggregationPreconBuilder(blocksize, kwargs)
+end
+
+function (b::SmoothedAggregationPreconBuilder)(A::AbstractSparseMatrixCSC, p)
+    return (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)), I)
+end
+
+
+"""
+   RugeStubenPreconBuilder(;blocksize=1, kwargs...)
+
+Return callable object constructing a  left algebraic multigrid preconditioner after Ruge & Stüben
+to be used with the `precs` API of LinearSolve.
+"""
+struct RugeStubenPreconBuilder{Tk}
+    blocksize::Int
+    kwargs::Tk
+end
+
+function RugeStubenPreconBuilder(; blocksize = 1, kwargs...)
+    return RugeStubenPreconBuilder(blocksize, kwargs)
+end
+
+function (b::RugeStubenPreconBuilder)(A::AbstractSparseMatrixCSC, p)
+    return (aspreconditioner(ruge_stuben(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)), I)
+end

test/runtests.jl

@@ -330,4 +330,27 @@ end
 X = poisson(27_000)+24.0*I
 ml = ruge_stuben(X)
 b = rand(27_000)
-@test AlgebraicMultigrid._solve(ml, b, reltol = 1e-10) ≈ X \ b rtol = 1e-10
+@test AlgebraicMultigrid._solve(ml, b, reltol = 1e-10) ≈ X \ b rtol = 1e-10
+
+
+
+
+# LinearSolve precs interface
+@testset "LinearSolvePrecs" begin
+
+for sz in [ (10,10), (20,20), (50,50)]
+    A = poisson(sz)
+    u0= ones(size(A,1))
+    b=A*u0
+
+    prob = LinearProblem(A, b)
+    strategy = KrylovJL_CG(precs = RugeStubenPreconBuilder())
+    @test solve(prob, strategy, atol=1.0e-14) ≈ u0 rtol = 1.0e-8
+
+
+    strategy = KrylovJL_CG(precs = SmoothedAggregationPreconBuilder())
+    @test solve(prob, strategy, atol=1.0e-14) ≈ u0 rtol = 1.0e-8
+    
+end
+
+end