Solver interfaces 1

Project.toml

@@ -4,6 +4,7 @@ authors = ["Santiago Badia <santiago.badia@monash.edu>", "Jordi Manyer <jordi.ma
 version = "0.2.0"
 
 [deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/GridapSolvers.jl

@@ -1,9 +1,11 @@
 module GridapSolvers
 
+  include("SolverInterfaces/SolverInterfaces.jl")
   include("MultilevelTools/MultilevelTools.jl")
   include("LinearSolvers/LinearSolvers.jl")
   include("PatchBasedSmoothers/PatchBasedSmoothers.jl")
 
+  using GridapSolvers.SolverInterfaces
   using GridapSolvers.MultilevelTools
   using GridapSolvers.LinearSolvers
   using GridapSolvers.PatchBasedSmoothers

src/LinearSolvers/Krylov/CGSolvers.jl

@@ -1,17 +1,18 @@
 
 struct CGSolver <: Gridap.Algebra.LinearSolver
   Pl       :: Gridap.Algebra.LinearSolver
-  maxiter  :: Int64
-  atol     :: Float64
-  rtol     :: Float64
+  log      :: ConvergenceLog{Float64}
   flexible :: Bool
-  verbose  :: Bool
 end
 
-function CGSolver(Pl;maxiter=10000,atol=1e-12,rtol=1.e-6,flexible=false,verbose=false)
-  return CGSolver(Pl,maxiter,atol,rtol,flexible,verbose)
+function CGSolver(Pl;maxiter=1000,atol=1e-12,rtol=1.e-6,flexible=false,verbose=0,name="CG")
+  tols = SolverTolerances{Float64}(;maxiter=maxiter,atol=atol,rtol=rtol)
+  log  = ConvergenceLog(name,tols;verbose=verbose)
+  return CGSolver(Pl,log,flexible)
 end
 
+AbstractTrees.children(s::CGSolver) = [s.Pl]
+
 struct CGSymbolicSetup <: Gridap.Algebra.SymbolicSetup
   solver
 end
@@ -49,27 +50,24 @@ end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup,b::AbstractVector)
   solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches
-  maxiter, atol, rtol = solver.maxiter, solver.atol, solver.rtol
-  flexible, verbose = solver.flexible, solver.verbose
+  flexible, log = solver.flexible, solver.log
   w,p,z,r = caches
-  verbose && println(" > Starting CG solver: ")
 
   # Initial residual
   mul!(w,A,x); r .= b .- w
   fill!(p,0.0); γ = 1.0
 
-  res  = norm(r); res_0 = res
-  iter = 0; converged = false
-  while !converged && (iter < maxiter)
-    verbose && println("   > Iteration ", iter," - Residual: ", res)
+  res  = norm(r)
+  done = init!(log,res)
+  while !done
 
     if !flexible # β = (zₖ₊₁ ⋅ rₖ₊₁)/(zₖ ⋅ rₖ)
       solve!(z, Pl, r)
       β = γ; γ = dot(z, r); β = γ / β
     else         # β = (zₖ₊₁ ⋅ (rₖ₊₁-rₖ))/(zₖ ⋅ rₖ)
-      β = γ; γ = dot(z, r)
+      δ = dot(z, r)
       solve!(z, Pl, r)
-      γ = dot(z, r) - γ; β = γ / β
+      β = γ; γ = dot(z, r); β = (γ-δ) / β
     end
     p .= z .+ β .* p
 
@@ -81,12 +79,10 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup,b::Abstrac
     x .+= α .* p
     r .-= α .* w
 
-    res = norm(r)
-    converged = (res < atol || res < rtol*res_0)
-    iter += 1
+    res  = norm(r)
+    done = update!(log,res)
   end
-  verbose && println("   > Num Iter: ", iter," - Final residual: ", res)
-  verbose && println("   Exiting CG solver.")
 
+  finalize!(log,res)
   return x
 end

src/LinearSolvers/Krylov/FGMRESSolvers.jl

@@ -1,18 +1,23 @@
 
 # FGMRES Solver
 struct FGMRESSolver <: Gridap.Algebra.LinearSolver
-  m       :: Int
-  Pr      :: Gridap.Algebra.LinearSolver
-  Pl      :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  atol    :: Float64
-  rtol    :: Float64
-  verbose :: Bool
+  m         :: Int
+  Pr        :: Gridap.Algebra.LinearSolver
+  Pl        :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  outer_log :: ConvergenceLog{Float64}
+  inner_log :: ConvergenceLog{Float64}
 end
 
-function FGMRESSolver(m,Pr;Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false)
-  return FGMRESSolver(m,Pr,Pl,atol,rtol,verbose)
+function FGMRESSolver(m,Pr;Pl=nothing,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="FGMRES")
+  outer_tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
+  outer_log  = ConvergenceLog(name,outer_tols,verbose=verbose)
+  inner_tols = SolverTolerances{Float64}(maxiter=m,atol=atol,rtol=rtol)
+  inner_log  = ConvergenceLog("$(name)_inner",inner_tols,verbose=verbose,nested=true)
+  return FGMRESSolver(m,Pr,Pl,outer_log,inner_log)
 end
 
+AbstractTrees.children(s::FGMRESSolver) = [s.Pr,s.Pl]
+
 struct FGMRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup
   solver
 end
@@ -30,11 +35,11 @@ mutable struct FGMRESNumericalSetup <: Gridap.Algebra.NumericalSetup
 end
 
 function get_solver_caches(solver::FGMRESSolver,A)
-  m = solver.m; Pl = solver.Pl
+  m = solver.m
 
   V  = [allocate_col_vector(A) for i in 1:m+1]
   Z  = [allocate_col_vector(A) for i in 1:m]
-  zl = !isa(Pl,Nothing) ? allocate_col_vector(A) : nothing
+  zl = allocate_col_vector(A)
 
   H = zeros(m+1,m)  # Hessenberg matrix
   g = zeros(m+1)    # Residual vector
@@ -61,26 +66,21 @@ end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::AbstractVector)
   solver, A, Pl, Pr, caches = ns.solver, ns.A, ns.Pl_ns, ns.Pr_ns, ns.caches
-  m, atol, rtol, verbose = solver.m, solver.atol, solver.rtol, solver.verbose
+  log, ilog = solver.outer_log, solver.inner_log
   V, Z, zl, H, g, c, s = caches
-  verbose && println(" > Starting FGMRES solver: ")
 
   # Initial residual
   krylov_residual!(V[1],x,A,b,Pl,zl)
-
-  iter = 0
-  β    = norm(V[1]); β0 = β
-  converged = (β < atol || β < rtol*β0)
-  while !converged
-    verbose && println("   > Iteration ", iter," - Residual: ", β)
-    fill!(H,0.0)
-
+  β    = norm(V[1])
+  done = init!(log,β)
+  while !done
     # Arnoldi process
     j = 1
     V[1] ./= β
+    fill!(H,0.0)
     fill!(g,0.0); g[1] = β
-    while ( j < m+1 && !converged )
-      verbose && println("      > Inner iteration ", j," - Residual: ", β)
+    idone = init!(ilog,β)
+    while !idone
       # Arnoldi orthogonalization by Modified Gram-Schmidt
       krylov_mul!(V[j+1],A,V[j],Pr,Pl,Z[j],zl)
       for i in 1:j
@@ -102,8 +102,9 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::Abs
       H[j,j] = c[j]*H[j,j] + s[j]*H[j+1,j]; H[j+1,j] = 0.0
       g[j+1] = -s[j]*g[j]; g[j] = c[j]*g[j]
 
-      β  = abs(g[j+1]); converged = (β < atol || β < rtol*β0)
+      β  = abs(g[j+1])
       j += 1
+      idone = update!(ilog,β)
     end
     j = j-1
 
@@ -117,11 +118,9 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::Abs
       x .+= g[i] .* Z[i]
     end
     krylov_residual!(V[1],x,A,b,Pl,zl)
-
-    iter += 1
+    done = update!(log,β)
   end
-  verbose && println("   > Num Iter: ", iter," - Final residual: ", β)
-  verbose && println("   Exiting FGMRES solver.")
 
+  finalize!(log,β)
   return x
 end

src/LinearSolvers/Krylov/GMRESSolvers.jl

@@ -1,17 +1,22 @@
 # GMRES Solver
 struct GMRESSolver <: Gridap.Algebra.LinearSolver
-  m       :: Int
-  Pr      :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  Pl      :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  atol    :: Float64
-  rtol    :: Float64
-  verbose :: Bool
+  m         :: Int
+  Pr        :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  Pl        :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  outer_log :: ConvergenceLog{Float64}
+  inner_log :: ConvergenceLog{Float64}
 end
 
-function GMRESSolver(m;Pr=nothing,Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false)
-  return GMRESSolver(m,Pr,Pl,atol,rtol,verbose)
+function GMRESSolver(m;Pr=nothing,Pl=nothing,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="GMRES")
+  outer_tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
+  outer_log  = ConvergenceLog(name,outer_tols,verbose=verbose)
+  inner_tols = SolverTolerances{Float64}(maxiter=m,atol=atol,rtol=rtol)
+  inner_log  = ConvergenceLog("$(name)_inner",inner_tols,verbose=verbose,nested=true)
+  return GMRESSolver(m,Pr,Pl,outer_log,inner_log)
 end
 
+AbstractTrees.children(s::GMRESSolver) = [s.Pr,s.Pl]
+
 struct GMRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup
   solver
 end
@@ -33,7 +38,7 @@ function get_solver_caches(solver::GMRESSolver,A)
 
   V  = [allocate_col_vector(A) for i in 1:m+1]
   zr = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing
-  zl = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing
+  zl = allocate_col_vector(A)
 
   H = zeros(m+1,m)  # Hessenberg matrix
   g = zeros(m+1)    # Residual vector
@@ -62,25 +67,21 @@ end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::AbstractVector)
   solver, A, Pl, Pr, caches = ns.solver, ns.A, ns.Pl_ns, ns.Pr_ns, ns.caches
-  m, atol, rtol, verbose = solver.m, solver.atol, solver.rtol, solver.verbose
+  log, ilog = solver.outer_log, solver.inner_log
   V, zr, zl, H, g, c, s = caches
-  verbose && println(" > Starting GMRES solver: ")
 
   # Initial residual
   krylov_residual!(V[1],x,A,b,Pl,zl)
-  β    = norm(V[1]); β0 = β
-  iter = 0
-  converged = (β < atol || β < rtol*β0)
-  while !converged
-    verbose && println("   > Iteration ", iter," - Residual: ", β)
-    fill!(H,0.0)
-
+  β    = norm(V[1])
+  done = init!(log,β)
+  while !done
     # Arnoldi process
     j = 1
     V[1] ./= β
+    fill!(H,0.0)
     fill!(g,0.0); g[1] = β
-    while ( j < m+1 && !converged )
-      verbose && println("      > Inner iteration ", j," - Residual: ", β)
+    idone = init!(ilog,β)
+    while !idone
       # Arnoldi orthogonalization by Modified Gram-Schmidt
       krylov_mul!(V[j+1],A,V[j],Pr,Pl,zr,zl)
       for i in 1:j
@@ -102,8 +103,9 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst
       H[j,j] = c[j]*H[j,j] + s[j]*H[j+1,j]; H[j+1,j] = 0.0
       g[j+1] = -s[j]*g[j]; g[j] = c[j]*g[j]
 
-      β  = abs(g[j+1]); converged = (β < atol || β < rtol*β0)
+      β  = abs(g[j+1])
       j += 1
+      idone = update!(ilog,β)
     end
     j = j-1
 
@@ -126,11 +128,8 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst
       x .+= zr
     end
     krylov_residual!(V[1],x,A,b,Pl,zl)
-
-    iter += 1
+    done = update!(log,β)
   end
-  verbose && println("   > Num Iter: ", iter," - Final residual: ", β)
-  verbose && println("   Exiting GMRES solver.")
-
+  finalize!(log,β)
   return x
-end
+end

src/LinearSolvers/Krylov/KrylovUtils.jl

@@ -31,7 +31,7 @@ function krylov_residual!(r,x,A,b,Pl,w)
   w .= b .- w
   solve!(r,Pl,w)
 end
-function krylov_residual!(r,x,A,b,Pl::Nothing,w::Nothing)
+function krylov_residual!(r,x,A,b,Pl::Nothing,w)
   mul!(r,A,x)
   r .= b .- r
 end

src/LinearSolvers/Krylov/MINRESSolvers.jl

@@ -1,16 +1,18 @@
 # MINRES Solver
 struct MINRESSolver <: Gridap.Algebra.LinearSolver
-  Pr      :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  Pl      :: Union{Gridap.Algebra.LinearSolver,Nothing}
-  atol    :: Float64
-  rtol    :: Float64
-  verbose :: Bool
+  Pr  :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  Pl  :: Union{Gridap.Algebra.LinearSolver,Nothing}
+  log :: ConvergenceLog{Float64}
 end
 
-function MINRESSolver(;Pr=nothing,Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false)
-  return MINRESSolver(Pr,Pl,atol,rtol,verbose)
+function MINRESSolver(;Pr=nothing,Pl=nothing,maxiter=1000,atol=1e-12,rtol=1.e-6,verbose=false,name="MINRES")
+  tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
+  log  = ConvergenceLog(name,tols,verbose=verbose)
+  return MINRESSolver(Pr,Pl,log)
 end
 
+AbstractTrees.children(s::MINRESSolver) = [s.Pr,s.Pl]
+
 struct MINRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup
   solver
 end
@@ -62,9 +64,8 @@ end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::MINRESNumericalSetup,b::AbstractVector)
   solver, A, Pl, Pr, caches = ns.solver, ns.A, ns.Pl_ns, ns.Pr_ns, ns.caches
-  atol, rtol, verbose = solver.atol, solver.rtol, solver.verbose
   V, W, zr, zl, H, g, c, s = caches
-  verbose && println(" > Starting MINRES solver: ")
+  log = solver.log
 
   Vjm1, Vj, Vjp1 = V
   Wjm1, Wj, Wjp1 = W
@@ -74,12 +75,10 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::MINRESNumericalSetup,b::Abs
   fill!(H,0.0), fill!(c,1.0); fill!(s,0.0); fill!(g,0.0)
 
   krylov_residual!(Vj,x,A,b,Pl,zl)
-  β    = norm(Vj); β0 = β; Vj ./= β; g[1] = β
-  iter = 0
-  converged = (β < atol || β < rtol*β0)
-  while !converged
-    verbose && println("   > Iteration ", iter," - Residual: ", β)
-
+  β    = norm(Vj); Vj ./= β; g[1] = β
+  done = init!(log,β)
+  while !done
+    # Lanczos process
     krylov_mul!(Vjp1,A,Vj,Pr,Pl,zr,zl)
     H[3] = dot(Vjp1,Vj)
     Vjp1 .= Vjp1 .- H[3] .* Vj .- H[2] .* Vjm1
@@ -105,14 +104,13 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::MINRESNumericalSetup,b::Abs
       x .+= g[1] .* zr
     end
 
-    β  = abs(g[2]); converged = (β < atol || β < rtol*β0)
+    β  = abs(g[2])
     Vjm1, Vj, Vjp1 = Vj, Vjp1, Vjm1
     Wjm1, Wj, Wjp1 = Wj, Wjp1, Wjm1
     g[1] = g[2]; H[2] = H[4];
-    iter += 1
+    done = update!(log,β)
   end
-  verbose && println("   > Num Iter: ", iter," - Final residual: ", β)
-  verbose && println("   Exiting MINRES solver.")
-
+
+  finalize!(log,β)
   return x
 end

src/LinearSolvers/LinearSolvers.jl

@@ -1,6 +1,7 @@
 module LinearSolvers
 
 using Printf
+using AbstractTrees
 using LinearAlgebra
 using SparseArrays
 using SparseMatricesCSR
@@ -17,6 +18,7 @@ using GridapPETSc
 
 using GridapDistributed
 using GridapSolvers.MultilevelTools
+using GridapSolvers.SolverInterfaces
 
 export JacobiLinearSolver
 export RichardsonSmoother