From 9eba04500ecd384cfee937007fbb1be01ead08f5 Mon Sep 17 00:00:00 2001
From: Fredrik Ekre <ekrefredrik@gmail.com>
Date: Thu, 28 Jul 2022 03:03:14 +0200
Subject: [PATCH] Add ParaSails preconditioner.

---
 gen/solver_options.jl |  1 +
 src/solver_options.jl | 21 +++++++++++++++++++++
 src/solvers.jl        | 27 +++++++++++++++++++++++++++
 test/runtests.jl      | 24 ++++++++++++++++++++++++
 4 files changed, 73 insertions(+)

diff --git a/gen/solver_options.jl b/gen/solver_options.jl
index 89b3c19..b30fe32 100644
--- a/gen/solver_options.jl
+++ b/gen/solver_options.jl
@@ -57,5 +57,6 @@ open(joinpath(@__DIR__, "..", "src", "solver_options.jl"), "w") do io
     generate_options(io, "BoomerAMG", "HYPRE_BoomerAMGSet")
     # generate_options(io, "FSAI", "HYPRE_FSAISet")
     generate_options(io, "GMRES", "HYPRE_ParCSRGMRESSet", "HYPRE_GMRESSet")
+    generate_options(io, "ParaSails", "HYPRE_ParCSRParaSailsSet", "HYPRE_ParaSailsSet")
     generate_options(io, "PCG", "HYPRE_ParCSRPCGSet", "HYPRE_PCGSet")
 end
diff --git a/src/solver_options.jl b/src/solver_options.jl
index 9fc7287..789026d 100644
--- a/src/solver_options.jl
+++ b/src/solver_options.jl
@@ -323,6 +323,27 @@ function Internals.set_options(s::GMRES, kwargs)
     end
 end
 
+function Internals.set_options(s::ParaSails, kwargs)
+    solver = s.solver
+    for (k, v) in kwargs
+        if k === :Filter
+            @check HYPRE_ParCSRParaSailsSetFilter(solver, v)
+        elseif k === :Loadbal
+            @check HYPRE_ParCSRParaSailsSetLoadbal(solver, v)
+        elseif k === :Logging
+            @check HYPRE_ParCSRParaSailsSetLogging(solver, v)
+        elseif k === :Params
+            @check HYPRE_ParCSRParaSailsSetParams(solver, v...)
+        elseif k === :Reuse
+            @check HYPRE_ParCSRParaSailsSetReuse(solver, v)
+        elseif k === :Sym
+            @check HYPRE_ParCSRParaSailsSetSym(solver, v)
+        else
+            throw(ArgumentError("unknown option $k for HYPRE.ParaSails"))
+        end
+    end
+end
+
 function Internals.set_options(s::PCG, kwargs)
     solver = s.solver
     for (k, v) in kwargs
diff --git a/src/solvers.jl b/src/solvers.jl
index bcb28d4..b5ceeee 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -227,6 +227,33 @@ function Internals.set_precond(gmres::GMRES, p::HYPRESolver)
 end
 
 
+#####################
+# (ParCSR)ParaSails #
+#####################
+
+mutable struct ParaSails <: HYPRESolver
+    comm::MPI.Comm
+    solver::HYPRE_Solver
+    function ParaSails(comm::MPI.Comm=MPI.COMM_WORLD; kwargs...)
+        # Note: comm is used in this solver so default to COMM_WORLD
+        solver = new(comm, C_NULL)
+        solver_ref = Ref{HYPRE_Solver}(C_NULL)
+        @check HYPRE_ParCSRParaSailsCreate(comm, solver_ref)
+        solver.solver = solver_ref[]
+        # Attach a finalizer
+        finalizer(x -> HYPRE_ParCSRParaSailsDestroy(x.solver), solver)
+        # Set the options
+        Internals.set_options(solver, kwargs)
+        return solver
+    end
+end
+
+const ParCSRParaSails = ParaSails
+
+Internals.setup_func(::ParaSails) = HYPRE_ParCSRParaSailsSetup
+Internals.solve_func(::ParaSails) = HYPRE_ParCSRParaSailsSolve
+
+
 ###############
 # (ParCSR)PCG #
 ###############
diff --git a/test/runtests.jl b/test/runtests.jl
index 179cc5e..3d7bde8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -409,6 +409,30 @@ end
     @test x ≈ A \ b atol=tol
 end
 
+@testset "(ParCSR)ParaSails" begin
+    # Solver constructor and options
+    @test_throws(
+        ArgumentError("unknown option UnknownOption for HYPRE.ParaSails"),
+        HYPRE.ParaSails(; UnknownOption = 1)
+    )
+    # Setup
+    A = sprand(100, 100, 0.05); A = A'A + 5I
+    b = rand(100)
+    x = zeros(100)
+    ilower, iupper = 1, size(A, 1)
+    A_h = HYPREMatrix(A, ilower, iupper)
+    b_h = HYPREVector(b, ilower, iupper)
+    x_h = HYPREVector(b, ilower, iupper)
+    # Solve with ParaSails as preconditioner
+    tol = 1e-9
+    parasails = HYPRE.ParaSails()
+    pcg = HYPRE.PCG(; Tol = tol, Precond = parasails)
+    HYPRE.solve!(pcg, x_h, A_h, b_h)
+    copy!(x, x_h)
+    # Test result with direct solver
+    @test x ≈ A \ b atol=tol
+end
+
 @testset "(ParCSR)PCG" begin
     # Solver constructor and options
     @test_throws(