From a43a913974f628c5f870608f56747e0f4b8ae01a Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 26 Sep 2023 13:20:04 +1000 Subject: [PATCH 1/5] Started implementing CG --- src/LinearSolvers/CGSolvers.jl | 127 +++++++++++++++++++++++++++++ src/LinearSolvers/MINRESSolvers.jl | 0 2 files changed, 127 insertions(+) create mode 100644 src/LinearSolvers/CGSolvers.jl create mode 100644 src/LinearSolvers/MINRESSolvers.jl diff --git a/src/LinearSolvers/CGSolvers.jl b/src/LinearSolvers/CGSolvers.jl new file mode 100644 index 00000000..d8d19bf4 --- /dev/null +++ b/src/LinearSolvers/CGSolvers.jl @@ -0,0 +1,127 @@ + +struct CGSolver <: Gridap.Algebra.LinearSolver + Pl ::Gridap.Algebra.LinearSolver + maxiter:: Int64 + atol ::Float64 + rtol ::Float64 + variant::Symbol + verbose::Bool +end + +function CGSolver(Pl;maxiter=10000,atol=1e-12,rtol=1.e-6,flexible=false,verbose=false) + variant = flexible ? :flexible : :standard + return CGSolver(Pl,maxiter,atol,rtol,variant,verbose) +end + +struct CGSymbolicSetup <: Gridap.Algebra.SymbolicSetup + solver +end + +function Gridap.Algebra.symbolic_setup(solver::CGSolver, A::AbstractMatrix) + return CGSymbolicSetup(solver) +end + +mutable struct CGNumericalSetup{T} <: Gridap.Algebra.NumericalSetup + solver + A + Pl_ns + caches +end + +function get_cg_caches(A) + w = allocate_col_vector(A) + p = allocate_col_vector(A) + r = allocate_col_vector(A) + return (w,p,r) +end + +function Gridap.Algebra.numerical_setup(ss::CGSymbolicSetup, A::AbstractMatrix) + solver = ss.solver + Pl_ns = numerical_setup(symbolic_setup(solver.Pl,A),A) + caches = get_cg_caches(A) + return CGNumericalSetup{solver.variant}(solver,A,Pl_ns,caches) +end + +function Gridap.Algebra.numerical_setup!(ns::CGNumericalSetup, A::AbstractMatrix) + numerical_setup!(ns.Pl_ns,A) + ns.A = A +end + +function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup{:standard},b::AbstractVector) + solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches + maxiter, atol, rtol, verbose = solver.maxiter, solver.atol, solver.rtol, solver.verbose + w,p,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) + + # Apply left preconditioner + solve!(z, Pl, r) + + # p := z + β⋅p , β = (zₖ₊₁ ⋅ rₖ₊₁)/(zₖ ⋅ rₖ) + β = γ; γ = dot(z, r); β = γ / β + p .= z .+ β .* p + + # w = A⋅p + mul!(w,A,p) + α = γ / dot(p, w) + + # Update solution and residual + x .+= α .* p + r .-= α .* w + + res = norm(r) + converged = (res < atol || res < rtol*res_0) + iter += 1 + end + verbose && println(" > Num Iter: ", iter," - Final residual: ", res) + verbose && println(" Exiting CG solver.") + + return x +end + +function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup{:flexible},b::AbstractVector) + solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches + maxiter, atol, rtol, verbose = solver.maxiter, solver.atol, solver.rtol, solver.verbose + w,p,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) + + # p := z + β⋅p , β = (zₖ₊₁ ⋅ (rₖ₊₁-rₖ))/(zₖ ⋅ rₖ) + β = γ; γ = dot(z, r) + solve!(z, Pl, r) + γ = dot(z, r) - γ; β = γ / β + p .= z .+ β .* p + + # w = A⋅p + mul!(w,A,p) + α = γ / dot(p, w) + + # Update solution and residual + x .+= α .* p + r .-= α .* w + + res = norm(r) + converged = (res < atol || res < rtol*res_0) + iter += 1 + end + verbose && println(" > Num Iter: ", iter," - Final residual: ", res) + verbose && println(" Exiting CG solver.") + + return x +end diff --git a/src/LinearSolvers/MINRESSolvers.jl b/src/LinearSolvers/MINRESSolvers.jl new file mode 100644 index 00000000..e69de29b From 2b907e940879b6f450d360472e290bec43f7fc49 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 29 Sep 2023 16:23:14 +1000 Subject: [PATCH 2/5] Added PCG and FPCG solvers --- src/LinearSolvers/CGSolvers.jl | 19 +++++----- src/LinearSolvers/LinearSolvers.jl | 1 + test/runtests.jl | 2 +- ...SSolversTests.jl => KrylovSolversTests.jl} | 36 ++++++++++++------- 4 files changed, 35 insertions(+), 23 deletions(-) rename test/seq/{GMRESSolversTests.jl => KrylovSolversTests.jl} (78%) diff --git a/src/LinearSolvers/CGSolvers.jl b/src/LinearSolvers/CGSolvers.jl index d8d19bf4..aaab94e3 100644 --- a/src/LinearSolvers/CGSolvers.jl +++ b/src/LinearSolvers/CGSolvers.jl @@ -1,11 +1,11 @@ struct CGSolver <: Gridap.Algebra.LinearSolver - Pl ::Gridap.Algebra.LinearSolver - maxiter:: Int64 - atol ::Float64 - rtol ::Float64 - variant::Symbol - verbose::Bool + Pl :: Gridap.Algebra.LinearSolver + maxiter :: Int64 + atol :: Float64 + rtol :: Float64 + variant :: Symbol + verbose :: Bool end function CGSolver(Pl;maxiter=10000,atol=1e-12,rtol=1.e-6,flexible=false,verbose=false) @@ -31,8 +31,9 @@ end function get_cg_caches(A) w = allocate_col_vector(A) p = allocate_col_vector(A) + z = allocate_col_vector(A) r = allocate_col_vector(A) - return (w,p,r) + return (w,p,z,r) end function Gridap.Algebra.numerical_setup(ss::CGSymbolicSetup, A::AbstractMatrix) @@ -50,7 +51,7 @@ end function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup{:standard},b::AbstractVector) solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches maxiter, atol, rtol, verbose = solver.maxiter, solver.atol, solver.rtol, solver.verbose - w,p,r = caches + w,p,z,r = caches verbose && println(" > Starting CG solver: ") # Initial residual @@ -90,7 +91,7 @@ end function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup{:flexible},b::AbstractVector) solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches maxiter, atol, rtol, verbose = solver.maxiter, solver.atol, solver.rtol, solver.verbose - w,p,r = caches + w,p,z,r = caches verbose && println(" > Starting CG solver: ") # Initial residual diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl index 14477214..ec7c547c 100644 --- a/src/LinearSolvers/LinearSolvers.jl +++ b/src/LinearSolvers/LinearSolvers.jl @@ -41,6 +41,7 @@ include("GMGLinearSolvers.jl") include("BlockDiagonalSmoothers.jl") include("IterativeLinearSolvers.jl") include("GMRESSolvers.jl") +include("CGSolvers.jl") include("SchurComplementSolvers.jl") end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index fb4835f7..4fd2e586 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,7 +65,7 @@ run_tests(joinpath(@__DIR__, "mpi")) # Sequential tests @time @testset "BlockDiagonalSmoothersTests" begin include("seq/BlockDiagonalSmoothersTests.jl") end @time @testset "DistributedPatchFESpacesTests" begin include("seq/DistributedPatchFESpacesTests.jl") end -@time @testset "GMRESSolversTests" begin include("seq/GMRESSolversTests.jl") end +@time @testset "KrylovSolversTests" begin include("seq/KrylovSolversTests.jl") end @time @testset "IterativeSolversTests" begin include("seq/IterativeSolversTests.jl") end @time @testset "PatchLinearSolverTests" begin include("seq/PatchLinearSolverTests.jl") end @time @testset "SymGaussSeidelSmoothersTests" begin include("seq/SymGaussSeidelSmoothersTests.jl") end diff --git a/test/seq/GMRESSolversTests.jl b/test/seq/KrylovSolversTests.jl similarity index 78% rename from test/seq/GMRESSolversTests.jl rename to test/seq/KrylovSolversTests.jl index 146d5e17..86547248 100644 --- a/test/seq/GMRESSolversTests.jl +++ b/test/seq/KrylovSolversTests.jl @@ -12,6 +12,20 @@ using GridapSolvers.LinearSolvers sol(x) = x[1] + x[2] f(x) = -Δ(sol)(x) +function test_solver(solver,op,Uh,dΩ) + A, b = get_matrix(op), get_vector(op); + ns = numerical_setup(symbolic_setup(solver,A),A) + + x = LinearSolvers.allocate_col_vector(A) + solve!(x,ns,b) + + u = interpolate(sol,Uh) + uh = FEFunction(Uh,x) + eh = uh - u + E = sum(∫(eh*eh)*dΩ) + @test E < 1.e-8 +end + function main(model) order = 1 qorder = order*2 + 1 @@ -24,29 +38,25 @@ function main(model) dΩ = Measure(Ω,qorder) a(u,v) = ∫(∇(v)⋅∇(u))*dΩ l(v) = ∫(v⋅f)*dΩ - op = AffineFEOperator(a,l,Uh,Vh) - A, b = get_matrix(op), get_vector(op); - + Pl = JacobiLinearSolver() - solver = LinearSolvers.GMRESSolver(40,Pl,1.e-8) - ns = numerical_setup(symbolic_setup(solver,A),A) - x = LinearSolvers.allocate_col_vector(A) - solve!(x,ns,b) + gmres = LinearSolvers.GMRESSolver(40,Pl;rtol=1.e-8,verbose=true) + test_solver(gmres,op,Uh,dΩ) - u = interpolate(sol,Uh) - uh = FEFunction(Uh,x) - eh = uh - u - E = sum(∫(eh*eh)*dΩ) - return E < 1.e-8 + pcg = LinearSolvers.CGSolver(Pl;verbose=true) + test_solver(pcg,op,Uh,dΩ) + + fpcg = LinearSolvers.CGSolver(Pl;flexible=true,verbose=true) + test_solver(fpcg,op,Uh,dΩ) end # Completely serial mesh_partition = (10,10) domain = (0,1,0,1) model = CartesianDiscreteModel(domain,mesh_partition) -@test main(model) +main(model) # Sequential num_ranks = (1,2) From 2f05929360209a7e2a2ed6ba0568300dfb6aa36e Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 30 Sep 2023 14:37:51 +1000 Subject: [PATCH 3/5] Added left preconditioning for GMRES --- src/LinearSolvers/CGSolvers.jl | 128 ------------------ src/LinearSolvers/Krylov/CGSolvers.jl | 92 +++++++++++++ src/LinearSolvers/Krylov/FGMRESSolvers.jl | 127 +++++++++++++++++ .../{ => Krylov}/GMRESSolvers.jl | 86 +++++++----- src/LinearSolvers/Krylov/KrylovUtils.jl | 37 +++++ src/LinearSolvers/Krylov/MINRESSolvers.jl | 121 +++++++++++++++++ src/LinearSolvers/LinearSolvers.jl | 14 +- src/LinearSolvers/MINRESSolvers.jl | 0 test/seq/KrylovSolversTests.jl | 11 +- 9 files changed, 446 insertions(+), 170 deletions(-) delete mode 100644 src/LinearSolvers/CGSolvers.jl create mode 100644 src/LinearSolvers/Krylov/CGSolvers.jl create mode 100644 src/LinearSolvers/Krylov/FGMRESSolvers.jl rename src/LinearSolvers/{ => Krylov}/GMRESSolvers.jl (57%) create mode 100644 src/LinearSolvers/Krylov/KrylovUtils.jl create mode 100644 src/LinearSolvers/Krylov/MINRESSolvers.jl delete mode 100644 src/LinearSolvers/MINRESSolvers.jl diff --git a/src/LinearSolvers/CGSolvers.jl b/src/LinearSolvers/CGSolvers.jl deleted file mode 100644 index aaab94e3..00000000 --- a/src/LinearSolvers/CGSolvers.jl +++ /dev/null @@ -1,128 +0,0 @@ - -struct CGSolver <: Gridap.Algebra.LinearSolver - Pl :: Gridap.Algebra.LinearSolver - maxiter :: Int64 - atol :: Float64 - rtol :: Float64 - variant :: Symbol - verbose :: Bool -end - -function CGSolver(Pl;maxiter=10000,atol=1e-12,rtol=1.e-6,flexible=false,verbose=false) - variant = flexible ? :flexible : :standard - return CGSolver(Pl,maxiter,atol,rtol,variant,verbose) -end - -struct CGSymbolicSetup <: Gridap.Algebra.SymbolicSetup - solver -end - -function Gridap.Algebra.symbolic_setup(solver::CGSolver, A::AbstractMatrix) - return CGSymbolicSetup(solver) -end - -mutable struct CGNumericalSetup{T} <: Gridap.Algebra.NumericalSetup - solver - A - Pl_ns - caches -end - -function get_cg_caches(A) - w = allocate_col_vector(A) - p = allocate_col_vector(A) - z = allocate_col_vector(A) - r = allocate_col_vector(A) - return (w,p,z,r) -end - -function Gridap.Algebra.numerical_setup(ss::CGSymbolicSetup, A::AbstractMatrix) - solver = ss.solver - Pl_ns = numerical_setup(symbolic_setup(solver.Pl,A),A) - caches = get_cg_caches(A) - return CGNumericalSetup{solver.variant}(solver,A,Pl_ns,caches) -end - -function Gridap.Algebra.numerical_setup!(ns::CGNumericalSetup, A::AbstractMatrix) - numerical_setup!(ns.Pl_ns,A) - ns.A = A -end - -function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup{:standard},b::AbstractVector) - solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches - maxiter, atol, rtol, verbose = solver.maxiter, solver.atol, solver.rtol, solver.verbose - 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) - - # Apply left preconditioner - solve!(z, Pl, r) - - # p := z + β⋅p , β = (zₖ₊₁ ⋅ rₖ₊₁)/(zₖ ⋅ rₖ) - β = γ; γ = dot(z, r); β = γ / β - p .= z .+ β .* p - - # w = A⋅p - mul!(w,A,p) - α = γ / dot(p, w) - - # Update solution and residual - x .+= α .* p - r .-= α .* w - - res = norm(r) - converged = (res < atol || res < rtol*res_0) - iter += 1 - end - verbose && println(" > Num Iter: ", iter," - Final residual: ", res) - verbose && println(" Exiting CG solver.") - - return x -end - -function Gridap.Algebra.solve!(x::AbstractVector,ns::CGNumericalSetup{:flexible},b::AbstractVector) - solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches - maxiter, atol, rtol, verbose = solver.maxiter, solver.atol, solver.rtol, solver.verbose - 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) - - # p := z + β⋅p , β = (zₖ₊₁ ⋅ (rₖ₊₁-rₖ))/(zₖ ⋅ rₖ) - β = γ; γ = dot(z, r) - solve!(z, Pl, r) - γ = dot(z, r) - γ; β = γ / β - p .= z .+ β .* p - - # w = A⋅p - mul!(w,A,p) - α = γ / dot(p, w) - - # Update solution and residual - x .+= α .* p - r .-= α .* w - - res = norm(r) - converged = (res < atol || res < rtol*res_0) - iter += 1 - end - verbose && println(" > Num Iter: ", iter," - Final residual: ", res) - verbose && println(" Exiting CG solver.") - - return x -end diff --git a/src/LinearSolvers/Krylov/CGSolvers.jl b/src/LinearSolvers/Krylov/CGSolvers.jl new file mode 100644 index 00000000..69767ddd --- /dev/null +++ b/src/LinearSolvers/Krylov/CGSolvers.jl @@ -0,0 +1,92 @@ + +struct CGSolver <: Gridap.Algebra.LinearSolver + Pl :: Gridap.Algebra.LinearSolver + maxiter :: Int64 + atol :: Float64 + rtol :: 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) +end + +struct CGSymbolicSetup <: Gridap.Algebra.SymbolicSetup + solver +end + +function Gridap.Algebra.symbolic_setup(solver::CGSolver, A::AbstractMatrix) + return CGSymbolicSetup(solver) +end + +mutable struct CGNumericalSetup <: Gridap.Algebra.NumericalSetup + solver + A + Pl_ns + caches +end + +function get_solver_caches(solver::CGSolver,A) + w = allocate_col_vector(A) + p = allocate_col_vector(A) + z = allocate_col_vector(A) + r = allocate_col_vector(A) + return (w,p,z,r) +end + +function Gridap.Algebra.numerical_setup(ss::CGSymbolicSetup, A::AbstractMatrix) + solver = ss.solver + Pl_ns = numerical_setup(symbolic_setup(solver.Pl,A),A) + caches = get_solver_caches(solver,A) + return CGNumericalSetup(solver,A,Pl_ns,caches) +end + +function Gridap.Algebra.numerical_setup!(ns::CGNumericalSetup, A::AbstractMatrix) + numerical_setup!(ns.Pl_ns,A) + ns.A = A +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 + 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) + + if !flexible # β = (zₖ₊₁ ⋅ rₖ₊₁)/(zₖ ⋅ rₖ) + solve!(z, Pl, r) + β = γ; γ = dot(z, r); β = γ / β + else # β = (zₖ₊₁ ⋅ (rₖ₊₁-rₖ))/(zₖ ⋅ rₖ) + β = γ; γ = dot(z, r) + solve!(z, Pl, r) + γ = dot(z, r) - γ; β = γ / β + end + p .= z .+ β .* p + + # w = A⋅p + mul!(w,A,p) + α = γ / dot(p, w) + + # Update solution and residual + x .+= α .* p + r .-= α .* w + + res = norm(r) + converged = (res < atol || res < rtol*res_0) + iter += 1 + end + verbose && println(" > Num Iter: ", iter," - Final residual: ", res) + verbose && println(" Exiting CG solver.") + + return x +end diff --git a/src/LinearSolvers/Krylov/FGMRESSolvers.jl b/src/LinearSolvers/Krylov/FGMRESSolvers.jl new file mode 100644 index 00000000..4aeae87b --- /dev/null +++ b/src/LinearSolvers/Krylov/FGMRESSolvers.jl @@ -0,0 +1,127 @@ + +# 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 +end + +function FGMRESSolver(m,Pr;Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false) + return FGMRESSolver(m,Pr,Pl,atol,rtol,verbose) +end + +struct FGMRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup + solver +end + +function Gridap.Algebra.symbolic_setup(solver::FGMRESSolver, A::AbstractMatrix) + return FGMRESSymbolicSetup(solver) +end + +mutable struct FGMRESNumericalSetup <: Gridap.Algebra.NumericalSetup + solver + A + Pr_ns + Pl_ns + caches +end + +function get_solver_caches(solver::FGMRESSolver,A) + m = solver.m; Pl = solver.Pl + + 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 + + H = zeros(m+1,m) # Hessenberg matrix + g = zeros(m+1) # Residual vector + c = zeros(m) # Gibens rotation cosines + s = zeros(m) # Gibens rotation sines + return (V,Z,zl,H,g,c,s) +end + +function Gridap.Algebra.numerical_setup(ss::FGMRESSymbolicSetup, A::AbstractMatrix) + solver = ss.solver + Pr_ns = numerical_setup(symbolic_setup(solver.Pr,A),A) + Pl_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A) + caches = get_solver_caches(solver,A) + return FGMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches) +end + +function Gridap.Algebra.numerical_setup!(ns::FGMRESNumericalSetup, A::AbstractMatrix) + numerical_setup!(ns.Pr_ns,A) + if !isa(ns.Pl_ns,Nothing) + numerical_setup!(ns.Pl_ns,A) + end + ns.A = A +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 + 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) + + # Arnoldi process + j = 1 + V[1] ./= β + fill!(g,0.0); g[1] = β + while ( j < m+1 && !converged ) + verbose && println(" > Inner iteration ", j," - Residual: ", β) + # Arnoldi orthogonalization by Modified Gram-Schmidt + krylov_mul!(V[j+1],A,V[j],Pr,Pl,Z[j],zl) + for i in 1:j + H[i,j] = dot(V[j+1],V[i]) + V[j+1] .= V[j+1] .- H[i,j] .* V[i] + end + H[j+1,j] = norm(V[j+1]) + V[j+1] ./= H[j+1,j] + + # Update QR + for i in 1:j-1 + γ = c[i]*H[i,j] + s[i]*H[i+1,j] + H[i+1,j] = -s[i]*H[i,j] + c[i]*H[i+1,j] + H[i,j] = γ + end + + # New Givens rotation, update QR and residual + c[j], s[j], _ = LinearAlgebra.givensAlgorithm(H[j,j],H[j+1,j]) + 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) + j += 1 + end + j = j-1 + + # Solve least squares problem Hy = g by backward substitution + for i in j:-1:1 + g[i] = (g[i] - dot(H[i,i+1:j],g[i+1:j])) / H[i,i] + end + + # Update solution & residual + for i in 1:j + x .+= g[i] .* Z[i] + end + krylov_residual!(V[1],x,A,b,Pl,zl) + + iter += 1 + end + verbose && println(" > Num Iter: ", iter," - Final residual: ", β) + verbose && println(" Exiting FGMRES solver.") + + return x +end diff --git a/src/LinearSolvers/GMRESSolvers.jl b/src/LinearSolvers/Krylov/GMRESSolvers.jl similarity index 57% rename from src/LinearSolvers/GMRESSolvers.jl rename to src/LinearSolvers/Krylov/GMRESSolvers.jl index fc2de310..6412eab8 100644 --- a/src/LinearSolvers/GMRESSolvers.jl +++ b/src/LinearSolvers/Krylov/GMRESSolvers.jl @@ -1,15 +1,15 @@ - # GMRES Solver struct GMRESSolver <: Gridap.Algebra.LinearSolver - m ::Int - Pl ::Gridap.Algebra.LinearSolver - atol::Float64 - rtol::Float64 - verbose::Bool + m :: Int + Pr :: Union{Gridap.Algebra.LinearSolver,Nothing} + Pl :: Union{Gridap.Algebra.LinearSolver,Nothing} + atol :: Float64 + rtol :: Float64 + verbose :: Bool end -function GMRESSolver(m,Pl;atol=1e-12,rtol=1.e-6,verbose=false) - return GMRESSolver(m,Pl,atol,rtol,verbose) +function GMRESSolver(m;Pr=nothing,Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false) + return GMRESSolver(m,Pr,Pl,atol,rtol,verbose) end struct GMRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup @@ -23,65 +23,72 @@ end mutable struct GMRESNumericalSetup <: Gridap.Algebra.NumericalSetup solver A + Pr_ns Pl_ns caches end -function get_gmres_caches(m,A) - w = allocate_col_vector(A) - V = [allocate_col_vector(A) for i in 1:m+1] - Z = [allocate_col_vector(A) for i in 1:m] +function get_solver_caches(solver::GMRESSolver,A) + m, Pl, Pr = solver.m, solver.Pl, solver.Pr + + V = [allocate_col_vector(A) for i in 1:m+1] + zr = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing + zl = !isa(Pl,Nothing) ? allocate_col_vector(A) : nothing H = zeros(m+1,m) # Hessenberg matrix g = zeros(m+1) # Residual vector c = zeros(m) # Gibens rotation cosines s = zeros(m) # Gibens rotation sines - return (w,V,Z,H,g,c,s) + return (V,zr,zl,H,g,c,s) end function Gridap.Algebra.numerical_setup(ss::GMRESSymbolicSetup, A::AbstractMatrix) solver = ss.solver - Pl_ns = numerical_setup(symbolic_setup(solver.Pl,A),A) - caches = get_gmres_caches(solver.m,A) - return GMRESNumericalSetup(solver,A,Pl_ns,caches) + Pr_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A) + Pl_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A) + caches = get_solver_caches(solver,A) + return GMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches) end function Gridap.Algebra.numerical_setup!(ns::GMRESNumericalSetup, A::AbstractMatrix) - numerical_setup!(ns.Pl_ns,A) + if !isa(ns.Pr_ns,Nothing) + numerical_setup!(ns.Pr_ns,A) + end + if !isa(ns.Pl_ns,Nothing) + numerical_setup!(ns.Pl_ns,A) + end ns.A = A end function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::AbstractVector) - solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches + 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 - w, V, Z, H, g, c, s = caches + V, zr, zl, H, g, c, s = caches verbose && println(" > Starting GMRES solver: ") # Initial residual - mul!(w,A,x); w .= b .- w - - β = norm(w); β0 = β - converged = (β < atol || β < rtol*β0) + 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) # Arnoldi process - fill!(g,0.0); g[1] = β - V[1] .= w ./ β j = 1 + V[1] ./= β + fill!(g,0.0); g[1] = β while ( j < m+1 && !converged ) verbose && println(" > Inner iteration ", j," - Residual: ", β) # Arnoldi orthogonalization by Modified Gram-Schmidt - solve!(Z[j],Pl,V[j]) - mul!(w,A,Z[j]) + krylov_mul!(V[j+1],A,V[j],Pr,Pl,zr,zl) for i in 1:j - H[i,j] = dot(w,V[i]) - w .= w .- H[i,j] .* V[i] + H[i,j] = dot(V[j+1],V[i]) + V[j+1] .= V[j+1] .- H[i,j] .* V[i] end - H[j+1,j] = norm(w) - V[j+1] = w ./ H[j+1,j] + H[j+1,j] = norm(V[j+1]) + V[j+1] ./= H[j+1,j] # Update QR for i in 1:j-1 @@ -106,10 +113,19 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst end # Update solution & residual - for i in 1:j - x .+= g[i] .* Z[i] + if isa(Pr,Nothing) + for i in 1:j + x .+= g[i] .* V[i] + end + else + fill!(zl,0.0) + for i in 1:j + zl .+= g[i] .* V[i] + end + solve!(zr,Pr,zl) + x .+= zr end - mul!(w,A,x); w .= b .- w + krylov_residual!(V[1],x,A,b,Pl,zl) iter += 1 end @@ -117,4 +133,4 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::GMRESNumericalSetup,b::Abst verbose && println(" Exiting GMRES solver.") return x -end +end \ No newline at end of file diff --git a/src/LinearSolvers/Krylov/KrylovUtils.jl b/src/LinearSolvers/Krylov/KrylovUtils.jl new file mode 100644 index 00000000..c8cde945 --- /dev/null +++ b/src/LinearSolvers/Krylov/KrylovUtils.jl @@ -0,0 +1,37 @@ + +""" + Computes the Krylov matrix-vector product y = Pl⁻¹⋅A⋅Pr⁻¹⋅x + by solving: + Pr⋅wr = x + wl = A⋅wr + Pl⋅y = wl +""" +function krylov_mul!(y,A,x,Pr,Pl,wr,wl) + solve!(wr,Pr,x) + mul!(wl,A,wr) + solve!(y,Pl,wl) +end +function krylov_mul!(y,A,x,Pr,Pl::Nothing,wr,wl) + solve!(wr,Pr,x) + mul!(y,A,wr) +end +function krylov_mul!(y,A,x,Pr::Nothing,Pl,wr,wl) + mul!(wl,A,x) + solve!(y,Pl,wl) +end + +""" + Computes the Krylov residual r = Pl⁻¹(A⋅x - b). + by solving: + w = A⋅x - b + Pl⋅r = w +""" +function krylov_residual!(r,x,A,b,Pl,w) + mul!(w,A,x) + w .= b .- w + solve!(r,Pl,w) +end +function krylov_residual!(r,x,A,b,Pl::Nothing,w::Nothing) + mul!(r,A,x) + r .= b .- r +end diff --git a/src/LinearSolvers/Krylov/MINRESSolvers.jl b/src/LinearSolvers/Krylov/MINRESSolvers.jl new file mode 100644 index 00000000..697bf0b2 --- /dev/null +++ b/src/LinearSolvers/Krylov/MINRESSolvers.jl @@ -0,0 +1,121 @@ + + +# MINRES Solver +struct MINRESSolver <: Gridap.Algebra.LinearSolver + m ::Int + Pl ::Gridap.Algebra.LinearSolver + atol::Float64 + rtol::Float64 + verbose::Bool +end + +function MINRESSolver(m,Pl;atol=1e-12,rtol=1.e-6,verbose=false) + return MINRESSolver(m,Pl,atol,rtol,verbose) +end + +struct MINRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup + solver +end + +function Gridap.Algebra.symbolic_setup(solver::MINRESSolver, A::AbstractMatrix) + return MINRESSymbolicSetup(solver) +end + +mutable struct MINRESNumericalSetup <: Gridap.Algebra.NumericalSetup + solver + A + Pl_ns + caches +end + +function get_MINRES_caches(m,A) + w = allocate_col_vector(A) + V = [allocate_col_vector(A) for i in 1:3] + Z = [allocate_col_vector(A) for i in 1:3] + + H = zeros(m+1,m) # Hessenberg matrix + g = zeros(m+1) # Residual vector + c = zeros(m) # Gibens rotation cosines + s = zeros(m) # Gibens rotation sines + return (w,V,Z,H,g,c,s) +end + +function Gridap.Algebra.numerical_setup(ss::MINRESSymbolicSetup, A::AbstractMatrix) + solver = ss.solver + Pl_ns = numerical_setup(symbolic_setup(solver.Pl,A),A) + caches = get_MINRES_caches(solver.m,A) + return MINRESNumericalSetup(solver,A,Pl_ns,caches) +end + +function Gridap.Algebra.numerical_setup!(ns::MINRESNumericalSetup, A::AbstractMatrix) + numerical_setup!(ns.Pl_ns,A) + ns.A = A +end + +function Gridap.Algebra.solve!(x::AbstractVector,ns::MINRESNumericalSetup,b::AbstractVector) + solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches + m, atol, rtol, verbose = solver.m, solver.atol, solver.rtol, solver.verbose + w, V, Z, H, g, c, s = caches + verbose && println(" > Starting MINRES solver: ") + + # Initial residual + mul!(w,A,x); w .= b .- w + + β = norm(w); β0 = β + converged = (β < atol || β < rtol*β0) + iter = 0 + while !converged + verbose && println(" > Iteration ", iter," - Residual: ", β) + fill!(H,0.0) + + # Arnoldi process + fill!(g,0.0); g[1] = β + V[1] .= w ./ β + j = 1 + + # Arnoldi orthogonalization by Modified Gram-Schmidt + solve!(Z[j],Pl,V[j]) + mul!(w,A,Z[j]) + for i in 1:j + H[i,j] = dot(w,V[i]) + w .= w .- H[i,j] .* V[i] + end + H[j+1,j] = norm(w) + V[j+1] = w ./ H[j+1,j] + + # Update QR + for i in 1:j-1 + γ = c[i]*H[i,j] + s[i]*H[i+1,j] + H[i+1,j] = -s[i]*H[i,j] + c[i]*H[i+1,j] + H[i,j] = γ + end + + # New Givens rotation, update QR and residual + c[j], s[j], _ = LinearAlgebra.givensAlgorithm(H[j,j],H[j+1,j]) + 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) + + # Solve least squares problem Hy = g by backward substitution + for i in j:-1:1 + g[i] = (g[i] - dot(H[i,i+1:j],g[i+1:j])) / H[i,i] + end + + # Update solution & residual + for i in 1:j + x .+= g[i] .* Z[i] + end + mul!(w,A,x); w .= b .- w + + iter += 1 + end + verbose && println(" > Num Iter: ", iter," - Final residual: ", β) + verbose && println(" Exiting MINRES solver.") + + return x +end + + + + diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl index ec7c547c..0db43168 100644 --- a/src/LinearSolvers/LinearSolvers.jl +++ b/src/LinearSolvers/LinearSolvers.jl @@ -23,6 +23,7 @@ export RichardsonSmoother export SymGaussSeidelSmoother export GMGLinearSolver export BlockDiagonalSmoother +export SchurComplementSolver # Wrappers for IterativeSolvers.jl export IS_ConjugateGradientSolver @@ -30,8 +31,17 @@ export IS_GMRESSolver export IS_MINRESSolver export IS_SSORSolver +# Krylov solvers +export CGSolver export GMRESSolver -export SchurComplementSolver +export FGMRESSolver +export MINRESSolver + +include("Krylov/KrylovUtils.jl") +include("Krylov/CGSolvers.jl") +include("Krylov/GMRESSolvers.jl") +include("Krylov/FGMRESSolvers.jl") +include("Krylov/MINRESSolvers.jl") include("IdentityLinearSolvers.jl") include("JacobiLinearSolvers.jl") @@ -40,8 +50,6 @@ include("SymGaussSeidelSmoothers.jl") include("GMGLinearSolvers.jl") include("BlockDiagonalSmoothers.jl") include("IterativeLinearSolvers.jl") -include("GMRESSolvers.jl") -include("CGSolvers.jl") include("SchurComplementSolvers.jl") end \ No newline at end of file diff --git a/src/LinearSolvers/MINRESSolvers.jl b/src/LinearSolvers/MINRESSolvers.jl deleted file mode 100644 index e69de29b..00000000 diff --git a/test/seq/KrylovSolversTests.jl b/test/seq/KrylovSolversTests.jl index 86547248..382df157 100644 --- a/test/seq/KrylovSolversTests.jl +++ b/test/seq/KrylovSolversTests.jl @@ -40,15 +40,18 @@ function main(model) l(v) = ∫(v⋅f)*dΩ op = AffineFEOperator(a,l,Uh,Vh) - Pl = JacobiLinearSolver() + P = JacobiLinearSolver() - gmres = LinearSolvers.GMRESSolver(40,Pl;rtol=1.e-8,verbose=true) + gmres = LinearSolvers.GMRESSolver(40;Pr=P,Pl=P,rtol=1.e-8,verbose=true) test_solver(gmres,op,Uh,dΩ) - pcg = LinearSolvers.CGSolver(Pl;verbose=true) + fgmres = LinearSolvers.FGMRESSolver(40,P;rtol=1.e-8,verbose=true) + test_solver(fgmres,op,Uh,dΩ) + + pcg = LinearSolvers.CGSolver(P;verbose=true) test_solver(pcg,op,Uh,dΩ) - fpcg = LinearSolvers.CGSolver(Pl;flexible=true,verbose=true) + fpcg = LinearSolvers.CGSolver(P;flexible=true,verbose=true) test_solver(fpcg,op,Uh,dΩ) end From ec564911fd0ccf0cca1a53c12788f4e7fd429099 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 2 Oct 2023 11:38:14 +1100 Subject: [PATCH 4/5] Non-preconditioned MINRES solver --- src/LinearSolvers/Krylov/MINRESSolvers.jl | 135 ++++++++++------------ test/seq/KrylovSolversTests.jl | 9 +- 2 files changed, 69 insertions(+), 75 deletions(-) diff --git a/src/LinearSolvers/Krylov/MINRESSolvers.jl b/src/LinearSolvers/Krylov/MINRESSolvers.jl index 697bf0b2..8a328377 100644 --- a/src/LinearSolvers/Krylov/MINRESSolvers.jl +++ b/src/LinearSolvers/Krylov/MINRESSolvers.jl @@ -1,16 +1,14 @@ - - # MINRES Solver struct MINRESSolver <: Gridap.Algebra.LinearSolver - m ::Int - Pl ::Gridap.Algebra.LinearSolver - atol::Float64 - rtol::Float64 - verbose::Bool + Pr :: Union{Gridap.Algebra.LinearSolver,Nothing} + Pl :: Union{Gridap.Algebra.LinearSolver,Nothing} + atol :: Float64 + rtol :: Float64 + verbose :: Bool end -function MINRESSolver(m,Pl;atol=1e-12,rtol=1.e-6,verbose=false) - return MINRESSolver(m,Pl,atol,rtol,verbose) +function MINRESSolver(;Pr=nothing,Pl=nothing,atol=1e-12,rtol=1.e-6,verbose=false) + return MINRESSolver(Pr,Pl,atol,rtol,verbose) end struct MINRESSymbolicSetup <: Gridap.Algebra.SymbolicSetup @@ -24,98 +22,91 @@ end mutable struct MINRESNumericalSetup <: Gridap.Algebra.NumericalSetup solver A + Pr_ns Pl_ns caches end -function get_MINRES_caches(m,A) - w = allocate_col_vector(A) - V = [allocate_col_vector(A) for i in 1:3] - Z = [allocate_col_vector(A) for i in 1:3] +function get_solver_caches(solver::MINRESSolver,A) + Pl, Pr = solver.Pl, solver.Pr + + V = [allocate_col_vector(A) for i in 1:3] + Z = [allocate_col_vector(A) for i in 1:3] + zr = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing + zl = !isa(Pl,Nothing) ? allocate_col_vector(A) : nothing - H = zeros(m+1,m) # Hessenberg matrix - g = zeros(m+1) # Residual vector - c = zeros(m) # Gibens rotation cosines - s = zeros(m) # Gibens rotation sines - return (w,V,Z,H,g,c,s) + H = zeros(4) # Hessenberg matrix + g = zeros(2) # Residual vector + c = zeros(2) # Gibens rotation cosines + s = zeros(2) # Gibens rotation sines + return (V,Z,zr,zl,H,g,c,s) end function Gridap.Algebra.numerical_setup(ss::MINRESSymbolicSetup, A::AbstractMatrix) solver = ss.solver - Pl_ns = numerical_setup(symbolic_setup(solver.Pl,A),A) - caches = get_MINRES_caches(solver.m,A) - return MINRESNumericalSetup(solver,A,Pl_ns,caches) + Pr_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A) + Pl_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A) + caches = get_solver_caches(solver,A) + return MINRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches) end function Gridap.Algebra.numerical_setup!(ns::MINRESNumericalSetup, A::AbstractMatrix) - numerical_setup!(ns.Pl_ns,A) + if !isa(ns.Pr_ns,Nothing) + numerical_setup!(ns.Pr_ns,A) + end + if !isa(ns.Pl_ns,Nothing) + numerical_setup!(ns.Pl_ns,A) + end ns.A = A end function Gridap.Algebra.solve!(x::AbstractVector,ns::MINRESNumericalSetup,b::AbstractVector) - solver, A, Pl, caches = ns.solver, ns.A, ns.Pl_ns, ns.caches - m, atol, rtol, verbose = solver.m, solver.atol, solver.rtol, solver.verbose - w, V, Z, H, g, c, s = caches + 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, Z, zr, zl, H, g, c, s = caches verbose && println(" > Starting MINRES solver: ") - # Initial residual - mul!(w,A,x); w .= b .- w + Vjm1, Vj, Vjp1 = V + Zjm1, Zj, Zjp1 = Z - β = norm(w); β0 = β - converged = (β < atol || β < rtol*β0) + fill!(Vjm1,0.0); fill!(Vjp1,0.0); copy!(Vj,b) + fill!(H,0.0), fill!(c,1.0); fill!(s,0.0); fill!(g,0.0) + + mul!(Vj,A,x,-1.0,1.0) + β = norm(Vj); β0 = β; Vj ./= β; g[1] = β iter = 0 + converged = (β < atol || β < rtol*β0) while !converged verbose && println(" > Iteration ", iter," - Residual: ", β) - fill!(H,0.0) - - # Arnoldi process - fill!(g,0.0); g[1] = β - V[1] .= w ./ β - j = 1 - - # Arnoldi orthogonalization by Modified Gram-Schmidt - solve!(Z[j],Pl,V[j]) - mul!(w,A,Z[j]) - for i in 1:j - H[i,j] = dot(w,V[i]) - w .= w .- H[i,j] .* V[i] - end - H[j+1,j] = norm(w) - V[j+1] = w ./ H[j+1,j] + + mul!(Vjp1,A,Vj) + H[3] = dot(Vjp1,Vj) + Vjp1 .= Vjp1 .- H[3] .* Vj .- H[2] .* Vjm1 + H[4] = norm(Vjp1) + Vjp1 ./= H[4] # Update QR - for i in 1:j-1 - γ = c[i]*H[i,j] + s[i]*H[i+1,j] - H[i+1,j] = -s[i]*H[i,j] + c[i]*H[i+1,j] - H[i,j] = γ - end + H[1] = s[1]*H[2]; H[2] = c[1]*H[2] + γ = c[2]*H[2] + s[2]*H[3]; H[3] = -s[2]*H[2] + c[2]*H[3]; H[2] = γ # New Givens rotation, update QR and residual - c[j], s[j], _ = LinearAlgebra.givensAlgorithm(H[j,j],H[j+1,j]) - 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) - - # Solve least squares problem Hy = g by backward substitution - for i in j:-1:1 - g[i] = (g[i] - dot(H[i,i+1:j],g[i+1:j])) / H[i,i] - end - - # Update solution & residual - for i in 1:j - x .+= g[i] .* Z[i] - end - mul!(w,A,x); w .= b .- w - + c[1], s[1] = c[2], s[2] + c[2], s[2], H[3] = LinearAlgebra.givensAlgorithm(H[3],H[4]) + g[2] = -s[2]*g[1]; g[1] = c[2]*g[1] + + # Update solution + Zjp1 .= Vj .- H[2] .* Zj .- H[1] .* Zjm1 + Zjp1 ./= H[3] + x .+= g[1] .* Zjp1 + + β = abs(g[2]); converged = (β < atol || β < rtol*β0) + Vjm1, Vj, Vjp1 = Vj, Vjp1, Vjm1 + Zjm1, Zj, Zjp1 = Zj, Zjp1, Zjm1 + g[1] = g[2]; H[2] = H[4]; iter += 1 end verbose && println(" > Num Iter: ", iter," - Final residual: ", β) verbose && println(" Exiting MINRES solver.") return x -end - - - - +end \ No newline at end of file diff --git a/test/seq/KrylovSolversTests.jl b/test/seq/KrylovSolversTests.jl index 382df157..cded6be9 100644 --- a/test/seq/KrylovSolversTests.jl +++ b/test/seq/KrylovSolversTests.jl @@ -48,11 +48,14 @@ function main(model) fgmres = LinearSolvers.FGMRESSolver(40,P;rtol=1.e-8,verbose=true) test_solver(fgmres,op,Uh,dΩ) - pcg = LinearSolvers.CGSolver(P;verbose=true) + pcg = LinearSolvers.CGSolver(P;rtol=1.e-8,verbose=true) test_solver(pcg,op,Uh,dΩ) - fpcg = LinearSolvers.CGSolver(P;flexible=true,verbose=true) + fpcg = LinearSolvers.CGSolver(P;flexible=true,rtol=1.e-8,verbose=true) test_solver(fpcg,op,Uh,dΩ) + + minres = LinearSolvers.MINRESSolver(;rtol=1.e-8,verbose=true) + test_solver(minres,op,Uh,dΩ) end # Completely serial @@ -68,6 +71,6 @@ parts = with_debug() do distribute end model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition) -@test main(model) +main(model) end \ No newline at end of file From 1ed6a392b0eb2b16934d2dfa42218ee608c6f0a8 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 2 Oct 2023 11:53:33 +1100 Subject: [PATCH 5/5] Added preconditioning to MINRES --- src/LinearSolvers/Krylov/GMRESSolvers.jl | 4 ++-- src/LinearSolvers/Krylov/MINRESSolvers.jl | 28 ++++++++++++++--------- test/seq/KrylovSolversTests.jl | 2 +- 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/LinearSolvers/Krylov/GMRESSolvers.jl b/src/LinearSolvers/Krylov/GMRESSolvers.jl index 6412eab8..92234193 100644 --- a/src/LinearSolvers/Krylov/GMRESSolvers.jl +++ b/src/LinearSolvers/Krylov/GMRESSolvers.jl @@ -33,7 +33,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(Pl,Nothing) ? allocate_col_vector(A) : nothing + zl = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing H = zeros(m+1,m) # Hessenberg matrix g = zeros(m+1) # Residual vector @@ -44,7 +44,7 @@ end function Gridap.Algebra.numerical_setup(ss::GMRESSymbolicSetup, A::AbstractMatrix) solver = ss.solver - Pr_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A) + Pr_ns = isa(solver.Pr,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A) Pl_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A) caches = get_solver_caches(solver,A) return GMRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches) diff --git a/src/LinearSolvers/Krylov/MINRESSolvers.jl b/src/LinearSolvers/Krylov/MINRESSolvers.jl index 8a328377..0f010aa1 100644 --- a/src/LinearSolvers/Krylov/MINRESSolvers.jl +++ b/src/LinearSolvers/Krylov/MINRESSolvers.jl @@ -31,7 +31,7 @@ function get_solver_caches(solver::MINRESSolver,A) Pl, Pr = solver.Pl, solver.Pr V = [allocate_col_vector(A) for i in 1:3] - Z = [allocate_col_vector(A) for i in 1:3] + W = [allocate_col_vector(A) for i in 1:3] zr = !isa(Pr,Nothing) ? allocate_col_vector(A) : nothing zl = !isa(Pl,Nothing) ? allocate_col_vector(A) : nothing @@ -39,12 +39,12 @@ function get_solver_caches(solver::MINRESSolver,A) g = zeros(2) # Residual vector c = zeros(2) # Gibens rotation cosines s = zeros(2) # Gibens rotation sines - return (V,Z,zr,zl,H,g,c,s) + return (V,W,zr,zl,H,g,c,s) end function Gridap.Algebra.numerical_setup(ss::MINRESSymbolicSetup, A::AbstractMatrix) solver = ss.solver - Pr_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A) + Pr_ns = isa(solver.Pr,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pr,A),A) Pl_ns = isa(solver.Pl,Nothing) ? nothing : numerical_setup(symbolic_setup(solver.Pl,A),A) caches = get_solver_caches(solver,A) return MINRESNumericalSetup(solver,A,Pr_ns,Pl_ns,caches) @@ -63,23 +63,24 @@ 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, Z, zr, zl, H, g, c, s = caches + V, W, zr, zl, H, g, c, s = caches verbose && println(" > Starting MINRES solver: ") Vjm1, Vj, Vjp1 = V - Zjm1, Zj, Zjp1 = Z + Wjm1, Wj, Wjp1 = W fill!(Vjm1,0.0); fill!(Vjp1,0.0); copy!(Vj,b) + fill!(Wjm1,0.0); fill!(Wjp1,0.0); fill!(Wj,0.0) fill!(H,0.0), fill!(c,1.0); fill!(s,0.0); fill!(g,0.0) - mul!(Vj,A,x,-1.0,1.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: ", β) - mul!(Vjp1,A,Vj) + krylov_mul!(Vjp1,A,Vj,Pr,Pl,zr,zl) H[3] = dot(Vjp1,Vj) Vjp1 .= Vjp1 .- H[3] .* Vj .- H[2] .* Vjm1 H[4] = norm(Vjp1) @@ -95,13 +96,18 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::MINRESNumericalSetup,b::Abs g[2] = -s[2]*g[1]; g[1] = c[2]*g[1] # Update solution - Zjp1 .= Vj .- H[2] .* Zj .- H[1] .* Zjm1 - Zjp1 ./= H[3] - x .+= g[1] .* Zjp1 + Wjp1 .= Vj .- H[2] .* Wj .- H[1] .* Wjm1 + Wjp1 ./= H[3] + if isa(Pr,Nothing) + x .+= g[1] .* Wjp1 + else + solve!(zr,Pr,Wjp1) + x .+= g[1] .* zr + end β = abs(g[2]); converged = (β < atol || β < rtol*β0) Vjm1, Vj, Vjp1 = Vj, Vjp1, Vjm1 - Zjm1, Zj, Zjp1 = Zj, Zjp1, Zjm1 + Wjm1, Wj, Wjp1 = Wj, Wjp1, Wjm1 g[1] = g[2]; H[2] = H[4]; iter += 1 end diff --git a/test/seq/KrylovSolversTests.jl b/test/seq/KrylovSolversTests.jl index cded6be9..d5452065 100644 --- a/test/seq/KrylovSolversTests.jl +++ b/test/seq/KrylovSolversTests.jl @@ -54,7 +54,7 @@ function main(model) fpcg = LinearSolvers.CGSolver(P;flexible=true,rtol=1.e-8,verbose=true) test_solver(fpcg,op,Uh,dΩ) - minres = LinearSolvers.MINRESSolver(;rtol=1.e-8,verbose=true) + minres = LinearSolvers.MINRESSolver(;Pl=P,Pr=P,rtol=1.e-8,verbose=true) test_solver(minres,op,Uh,dΩ) end