From e907b27dae77fed27ec78225907949774de4cb76 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 28 Jul 2022 12:48:46 +0200 Subject: [PATCH] Add Hybrid solver. --- gen/solver_options.jl | 1 + src/solver_options.jl | 110 ++++++++++++++++++++++++++++++++++++++++++ src/solvers.jl | 36 ++++++++++++++ test/runtests.jl | 41 ++++++++++++++++ 4 files changed, 188 insertions(+) diff --git a/gen/solver_options.jl b/gen/solver_options.jl index 961e51a..aefd3ed 100644 --- a/gen/solver_options.jl +++ b/gen/solver_options.jl @@ -58,6 +58,7 @@ open(joinpath(@__DIR__, "..", "src", "solver_options.jl"), "w") do io generate_options(io, "FlexGMRES", "HYPRE_ParCSRFlexGMRESSet", "HYPRE_FlexGMRESSet") # generate_options(io, "FSAI", "HYPRE_FSAISet") generate_options(io, "GMRES", "HYPRE_ParCSRGMRESSet", "HYPRE_GMRESSet") + generate_options(io, "Hybrid", "HYPRE_ParCSRHybridSet") generate_options(io, "ILU", "HYPRE_ILUSet") generate_options(io, "ParaSails", "HYPRE_ParCSRParaSailsSet", "HYPRE_ParaSailsSet") generate_options(io, "PCG", "HYPRE_ParCSRPCGSet", "HYPRE_PCGSet") diff --git a/src/solver_options.jl b/src/solver_options.jl index a3405df..2885d59 100644 --- a/src/solver_options.jl +++ b/src/solver_options.jl @@ -353,6 +353,116 @@ function Internals.set_options(s::GMRES, kwargs) end end +function Internals.set_options(s::Hybrid, kwargs) + solver = s.solver + for (k, v) in kwargs + if k === :AbsoluteTol + @check HYPRE_ParCSRHybridSetAbsoluteTol(solver, v) + elseif k === :AggInterpType + @check HYPRE_ParCSRHybridSetAggInterpType(solver, v) + elseif k === :AggNumLevels + @check HYPRE_ParCSRHybridSetAggNumLevels(solver, v) + elseif k === :CoarsenType + @check HYPRE_ParCSRHybridSetCoarsenType(solver, v) + elseif k === :ConvergenceTol + @check HYPRE_ParCSRHybridSetConvergenceTol(solver, v) + elseif k === :CycleNumSweeps + @check HYPRE_ParCSRHybridSetCycleNumSweeps(solver, v...) + elseif k === :CycleRelaxType + @check HYPRE_ParCSRHybridSetCycleRelaxType(solver, v...) + elseif k === :CycleType + @check HYPRE_ParCSRHybridSetCycleType(solver, v) + elseif k === :DSCGMaxIter + @check HYPRE_ParCSRHybridSetDSCGMaxIter(solver, v) + elseif k === :DofFunc + @check HYPRE_ParCSRHybridSetDofFunc(solver, v) + elseif k === :GridRelaxPoints + @check HYPRE_ParCSRHybridSetGridRelaxPoints(solver, v) + elseif k === :GridRelaxType + @check HYPRE_ParCSRHybridSetGridRelaxType(solver, v) + elseif k === :InterpType + @check HYPRE_ParCSRHybridSetInterpType(solver, v) + elseif k === :KDim + @check HYPRE_ParCSRHybridSetKDim(solver, v) + elseif k === :KeepTranspose + @check HYPRE_ParCSRHybridSetKeepTranspose(solver, v) + elseif k === :LevelOuterWt + @check HYPRE_ParCSRHybridSetLevelOuterWt(solver, v...) + elseif k === :LevelRelaxWt + @check HYPRE_ParCSRHybridSetLevelRelaxWt(solver, v...) + elseif k === :Logging + @check HYPRE_ParCSRHybridSetLogging(solver, v) + elseif k === :MaxCoarseSize + @check HYPRE_ParCSRHybridSetMaxCoarseSize(solver, v) + elseif k === :MaxLevels + @check HYPRE_ParCSRHybridSetMaxLevels(solver, v) + elseif k === :MaxRowSum + @check HYPRE_ParCSRHybridSetMaxRowSum(solver, v) + elseif k === :MeasureType + @check HYPRE_ParCSRHybridSetMeasureType(solver, v) + elseif k === :MinCoarseSize + @check HYPRE_ParCSRHybridSetMinCoarseSize(solver, v) + elseif k === :Nodal + @check HYPRE_ParCSRHybridSetNodal(solver, v) + elseif k === :NonGalerkinTol + @check HYPRE_ParCSRHybridSetNonGalerkinTol(solver, v...) + elseif k === :NumFunctions + @check HYPRE_ParCSRHybridSetNumFunctions(solver, v) + elseif k === :NumGridSweeps + @check HYPRE_ParCSRHybridSetNumGridSweeps(solver, v) + elseif k === :NumPaths + @check HYPRE_ParCSRHybridSetNumPaths(solver, v) + elseif k === :NumSweeps + @check HYPRE_ParCSRHybridSetNumSweeps(solver, v) + elseif k === :Omega + @check HYPRE_ParCSRHybridSetOmega(solver, v) + elseif k === :OuterWt + @check HYPRE_ParCSRHybridSetOuterWt(solver, v) + elseif k === :PCGMaxIter + @check HYPRE_ParCSRHybridSetPCGMaxIter(solver, v) + elseif k === :PMaxElmts + @check HYPRE_ParCSRHybridSetPMaxElmts(solver, v) + elseif k === :Precond + Internals.set_precond_defaults(v) + Internals.set_precond(s, v) + elseif k === :PrintLevel + @check HYPRE_ParCSRHybridSetPrintLevel(solver, v) + elseif k === :RecomputeResidual + @check HYPRE_ParCSRHybridSetRecomputeResidual(solver, v) + elseif k === :RecomputeResidualP + @check HYPRE_ParCSRHybridSetRecomputeResidualP(solver, v) + elseif k === :RelChange + @check HYPRE_ParCSRHybridSetRelChange(solver, v) + elseif k === :RelaxOrder + @check HYPRE_ParCSRHybridSetRelaxOrder(solver, v) + elseif k === :RelaxType + @check HYPRE_ParCSRHybridSetRelaxType(solver, v) + elseif k === :RelaxWeight + @check HYPRE_ParCSRHybridSetRelaxWeight(solver, v) + elseif k === :RelaxWt + @check HYPRE_ParCSRHybridSetRelaxWt(solver, v) + elseif k === :SeqThreshold + @check HYPRE_ParCSRHybridSetSeqThreshold(solver, v) + elseif k === :SetupType + @check HYPRE_ParCSRHybridSetSetupType(solver, v) + elseif k === :SolverType + @check HYPRE_ParCSRHybridSetSolverType(solver, v) + elseif k === :StopCrit + @check HYPRE_ParCSRHybridSetStopCrit(solver, v) + elseif k === :StrongThreshold + @check HYPRE_ParCSRHybridSetStrongThreshold(solver, v) + elseif k === :Tol + @check HYPRE_ParCSRHybridSetTol(solver, v) + elseif k === :TruncFactor + @check HYPRE_ParCSRHybridSetTruncFactor(solver, v) + elseif k === :TwoNorm + @check HYPRE_ParCSRHybridSetTwoNorm(solver, v) + else + throw(ArgumentError("unknown option $k for HYPRE.Hybrid")) + end + end +end + function Internals.set_options(s::ILU, kwargs) solver = s.solver for (k, v) in kwargs diff --git a/src/solvers.jl b/src/solvers.jl index 5435292..0928f7a 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -265,6 +265,42 @@ function Internals.set_precond(gmres::GMRES, p::HYPRESolver) end +########## +# Hybrid # +########## + +mutable struct Hybrid <: HYPRESolver + solver::HYPRE_Solver + function Hybrid(; kwargs...) + solver = new(C_NULL) + solver_ref = Ref{HYPRE_Solver}(C_NULL) + @check HYPRE_ParCSRHybridCreate(solver_ref) + solver.solver = solver_ref[] + # Attach a finalizer + finalizer(x -> HYPRE_ParCSRHybridDestroy(x.solver), solver) + # Set the options + Internals.set_options(solver, kwargs) + return solver + end +end + +function solve!(hybrid::Hybrid, x::HYPREVector, A::HYPREMatrix, b::HYPREVector) + @check HYPRE_ParCSRHybridSetup(hybrid.solver, A.parmatrix, b.parvector, x.parvector) + @check HYPRE_ParCSRHybridSolve(hybrid.solver, A.parmatrix, b.parvector, x.parvector) + return x +end + +Internals.setup_func(::Hybrid) = HYPRE_ParCSRHybridSetup +Internals.solve_func(::Hybrid) = HYPRE_ParCSRHybridSolve + +function Internals.set_precond(hybrid::Hybrid, p::HYPRESolver) + solve_f = Internals.solve_func(p) + setup_f = Internals.setup_func(p) + @check HYPRE_ParCSRHybridSetPrecond(hybrid.solver, solve_f, setup_f, p.solver) + return nothing +end + + ####### # ILU # ####### diff --git a/test/runtests.jl b/test/runtests.jl index cbc3acf..93ddede 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -449,6 +449,47 @@ end @test x ≈ A \ b atol=tol end +@testset "Hybrid" begin + # Solver constructor and options + @test_throws( + ArgumentError("unknown option UnknownOption for HYPRE.Hybrid"), + HYPRE.Hybrid(; UnknownOption = 1) + ) + # Setup + A = sprand(100, 100, 0.05); A = A'A + 5I + b = rand(100) + x = zeros(100) + A_h = HYPREMatrix(A) + b_h = HYPREVector(b) + x_h = HYPREVector(x) + # Solve + tol = 1e-9 + hybrid = HYPRE.Hybrid(; Tol = tol) + HYPRE.solve!(hybrid, x_h, A_h, b_h) + copy!(x, x_h) + # Test result with direct solver + @test x ≈ A \ b atol=tol + # Test without passing initial guess + x_h = HYPRE.solve(hybrid, A_h, b_h) + copy!(x, x_h) + @test x ≈ A \ b atol=tol + + # Solve with given preconditioner + # XXX: https://github.com/hypre-space/hypre/issues/699 + precond = HYPRE.BoomerAMG() + hybrid = HYPRE.Hybrid(; Tol = tol, SolverType = 3, #= Precond = precond =#) + x_h = HYPREVector(zeros(100)) + HYPRE.solve!(hybrid, x_h, A_h, b_h) + copy!(x, x_h) + # Test result with direct solver + @test x ≈ A \ b atol=tol + # Test without passing initial guess + x_h = HYPRE.solve(hybrid, A_h, b_h) + copy!(x, x_h) + @test x ≈ A \ b atol=tol +end + + @testset "ILU" begin # Solver constructor and options @test_throws(