Added linear solver wrappers for lu and IterativeSolvers.cg

Author: fverdugo

PartitionedSolvers/src/PartitionedSolvers.jl

@@ -3,6 +3,7 @@ module PartitionedSolvers
 using PartitionedArrays
 using SparseArrays
 using LinearAlgebra
+using IterativeSolvers
 
 export setup
 export solve!

PartitionedSolvers/src/interfaces.jl

@@ -28,6 +28,10 @@ struct LinearSolver{A,B,C,D} <: AbstractLinearSolver
     finalize!::D
 end
 
+function linear_solver(s::LinearSolver)
+    s
+end
+
 struct Preconditioner{A,B}
     solver::A
     solver_setup::B

PartitionedSolvers/src/smoothers.jl

@@ -179,3 +179,35 @@ function additive_schwarz_correction(local_solver)
     linear_solver(;setup,update!,solve!,finalize!)
 end
 
+# Wrappers
+
+function linear_solver(::typeof(LinearAlgebra.lu))
+    lu_solver()
+end
+
+function linear_solver(::typeof(IterativeSolvers.cg);Pl,kwargs...)
+    function setup(x,A,b,options)
+        Pl_solver = linear_solver(Pl)
+        P = PartitionedSolvers.setup(Pl_solver,x,A,b;options...)
+        A_ref = Ref(A)
+        (;P,A_ref)
+    end
+    function update!(state,A,options)
+        (;P,A_ref) = state
+        A_ref[] = A
+        P = PartitionedSolvers.update!(P,A;options...)
+        state
+    end
+    function solve!(x,state,b,options)
+        (;P,A_ref) = state
+        A = A_ref[]
+        IterativeSolvers.cg!(x,A,b;Pl=P,kwargs...)
+    end
+    function finalize!(state,A,options)
+        (;P) = state
+        PartitionedSolvers.finalize!(P)
+        nothing
+    end
+    linear_solver(;setup,update!,solve!,finalize!)
+end
+

PartitionedSolvers/test/amg_tests.jl

@@ -6,6 +6,7 @@ using PartitionedSolvers
 using LinearAlgebra
 using Test
 using IterativeSolvers: cg!
+import IterativeSolvers
 
 # First with a sequential matrix
 nodes_per_dir = (100,100)
@@ -65,6 +66,12 @@ Pl = setup(amg(;fine_params),y,A,b;nullspace=B)
 y .= 0
 cg!(y,A,b;Pl,verbose=true)
 
+solver = linear_solver(IterativeSolvers.cg;Pl=amg(;fine_params),verbose=true)
+S = setup(solver,y,A,b)
+solve!(y,S,b)
+update!(S,2*A)
+solve!(y,S,b)
+
 # Now in parallel
 
 parts_per_dir = (2,2)

PartitionedSolvers/test/smoothers_tests.jl

@@ -26,6 +26,17 @@ solve!(y,S,b)
 @test norm(y-x/2)/norm(x/2) < tol
 finalize!(S)
 
+solver = linear_solver(LinearAlgebra.lu)
+y = similar(x)
+S = setup(solver,y,A,b)
+solve!(y,S,b)
+tol = 1.e-8
+@test norm(y-x)/norm(x) < tol
+update!(S,2*A)
+solve!(y,S,b)
+@test norm(y-x/2)/norm(x/2) < tol
+finalize!(S)
+
 solver = richardson(lu_solver(),iters=1)
 y = similar(x)
 y .= 0