Merge pull request #175 from fverdugo/matrix_attributes

Release PartitionedSolvers 0.2.2

Author: fverdugo

PartitionedSolvers/CHANGELOG.md

@@ -6,7 +6,16 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.2.1] - 2025-07-26
+## [0.2.2] - 2024-10-03
+
+### Added
+
+- Support for user-defined near nullspace in AMG solver.
+- Flag `history` to return solve history in `solve!`.
+- Expose iterations for iterative solvers with function `iterations!`.
+- Solver traits `uses_nullspace` and `uses_initial_guess`. These allow the user to build the nullspace or reset the initial guess to zero only when needed by the solver.
+
+## [0.2.1] - 2024-07-26
 
 ### Added
 - Support for PartitionedArrays v0.5.

PartitionedSolvers/Project.toml

@@ -1,7 +1,7 @@
 name = "PartitionedSolvers"
 uuid = "11b65f7f-80ac-401b-9ef2-3db765482d62"
 authors = ["Francesc Verdugo <[email protected]>"]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"

PartitionedSolvers/src/PartitionedSolvers.jl

@@ -1,6 +1,7 @@
 module PartitionedSolvers
 
 using PartitionedArrays
+using PartitionedArrays: val_parameter
 using SparseArrays
 using LinearAlgebra
 using IterativeSolvers
@@ -13,6 +14,9 @@ export AbstractLinearSolver
 export linear_solver
 export default_nullspace
 export nullspace
+export uses_nullspace
+export uses_initial_guess
+export iterations!
 include("interfaces.jl")
 
 export lu_solver

PartitionedSolvers/src/amg.jl

@@ -789,7 +789,8 @@ function amg(;
     update! = amg_update!
     solve! = amg_solve!
     finalize! = amg_finalize!
-    linear_solver(;setup,update!,solve!,finalize!)
+    uses_nullspace = Val(true)
+    linear_solver(;setup,update!,solve!,finalize!,uses_nullspace)
 end
 
 function amg_setup(x,A,b,::Nothing,amg_params)
@@ -939,4 +940,4 @@ function amg_finalize!(setup)
     end
     coarse_solver_setup = setup.coarse_level.coarse_solver_setup
     finalize!(coarse_solver_setup)
-end
+end

PartitionedSolvers/src/interfaces.jl

@@ -1,4 +1,9 @@
 
+abstract type AbstractType end
+function Base.show(io::IO,data::AbstractType)
+    print(io,"PartitionedSolvers.$(nameof(typeof(data)))(…)")
+end
+
 function default_nullspace(A)
     T = eltype(A)
     [ones(T,size(A,2))]
@@ -10,29 +15,51 @@ function default_nullspace(A::PSparseMatrix)
     [ pones(T,col_partition) ]
 end
 
-abstract type AbstractLinearSolver end
+abstract type AbstractLinearSolver <: AbstractType end
 
 function linear_solver(;
         setup,
         solve!,
         update!,
         finalize! = ls_setup->nothing,
+        step! = nothing,
+        uses_nullspace = Val(false),
+        uses_initial_guess = Val(true),
+        returns_history = Val(false),
     )
-    LinearSolver(setup,solve!,update!,finalize!)
+    if step! === nothing
+        step! = (x,ls_setup,b,options,step=0) -> begin
+            if step !=0
+                return nothing
+            end
+            x = solve!(x,ls_setup,b,options)
+            x,step+1
+        end
+    end
+    traits = LinearSolverTraits(uses_nullspace,uses_initial_guess,returns_history)
+    LinearSolver(setup,solve!,update!,finalize!,step!,traits)
 end
 
-struct LinearSolver{A,B,C,D} <: AbstractLinearSolver
+struct LinearSolverTraits{A,B,C} <: AbstractType
+    uses_nullspace::A
+    uses_initial_guess::B
+    returns_history::C
+end
+
+struct LinearSolver{A,B,C,D,E,F} <: AbstractLinearSolver
     setup::A
     solve!::B
     update!::C
     finalize!::D
+    step!::E
+    traits::F
 end
 
 function linear_solver(s::LinearSolver)
     s
 end
 
-struct Preconditioner{A,B}
+struct Preconditioner{A,B} <: AbstractType
     solver::A
     solver_setup::B
 end
@@ -45,6 +72,18 @@ function nullspace(options)
     options.nullspace
 end
 
+function uses_nullspace(a::LinearSolver)
+    val_parameter(a.traits.uses_nullspace)
+end
+
+function uses_initial_guess(a::LinearSolver)
+    val_parameter(a.traits.uses_initial_guess)
+end
+
+function returns_history(a::LinearSolver)
+    val_parameter(a.traits.returns_history)
+end
+
 function setup(solver::LinearSolver,x,A,b;kwargs...)
     options = setup_options(;kwargs...)
     solver_setup = solver.setup(x,A,b,options)
@@ -57,18 +96,30 @@ function update!(P::Preconditioner,A;kwargs...)
     P
 end
 
-function solve_options(;zero_guess=false)
-    options = (;zero_guess)
+function solve_options(;zero_guess=false,history=Val(false))
+    options = (;zero_guess,history)
 end
 
 function solve!(x,P::Preconditioner,b;kwargs...)
     options = solve_options(;kwargs...)
-    P.solver.solve!(x,P.solver_setup,b,options)
-    x
+    next = P.solver.solve!(x,P.solver_setup,b,options)
+    if returns_history(P.solver)
+        x,log = next
+    else
+        x = next
+        log = nothing
+    end
+    if val_parameter(options.history) == true
+        return x, log
+    else
+        return x
+    end
 end
 
 function LinearAlgebra.ldiv!(x,P::Preconditioner,b)
-    fill!(x,zero(eltype(x)))
+    if uses_initial_guess(P.solver)
+        fill!(x,zero(eltype(x)))
+    end
     solve!(x,P,b;zero_guess=true)
     x
 end
@@ -77,3 +128,31 @@ function finalize!(P::Preconditioner)
     P.solver.finalize!(P.solver_setup)
 end
 
+function iterations!(x,P::Preconditioner,b;kwargs...)
+    options = solve_options(;kwargs...)
+    params = (;options,x,P,b)
+    LinearSolverIterator(params)
+end
+
+struct LinearSolverIterator{A} <: AbstractType
+    params::A
+end
+
+function Base.iterate(a::LinearSolverIterator)
+    P = a.params.P
+    options = a.params.options
+    b = a.params.b
+    x = a.params.x
+    next = P.solver.step!(x,P.solver_setup,b,options)
+    next
+end
+
+function Base.iterate(a::LinearSolverIterator,state)
+    P = a.params.P
+    options = a.params.options
+    b = a.params.b
+    x = a.params.x
+    next = P.solver.step!(x,P.solver_setup,b,options,state)
+    next
+end
+

PartitionedSolvers/src/smoothers.jl

@@ -3,16 +3,19 @@ function lu_solver()
     setup(x,op,b,options) = lu(op)
     update!(state,op,options) = lu!(state,op)
     solve!(x,P,b,options) = ldiv!(x,P,b)
-    linear_solver(;setup,solve!,update!)
+    uses_initial_guess = Val(false)
+    linear_solver(;setup,solve!,update!,uses_initial_guess)
 end
 
 function jacobi_correction()
     setup(x,op,b,options) = dense_diag!(similar(b),op)
     update!(state,op,options) = dense_diag!(state,op)
     function solve!(x,state,b,options)
         x .= state .\ b
+        x
     end
-    linear_solver(;setup,update!,solve!)
+    uses_initial_guess = Val(false)
+    linear_solver(;setup,update!,solve!,uses_initial_guess)
 end
 
 function richardson(solver;iters,omega=1)
@@ -39,13 +42,27 @@ function richardson(solver;iters,omega=1)
             ldiv!(dx,P,r)
             x .-= omega .* dx
         end
-        (;iters)
+        x, (;iters)
+    end
+    function step!(x,state,b,options,step=0)
+        if step == iters
+            return nothing
+        end
+        (r,dx,P,A_ref) = state
+        A = A_ref[]
+        dx .= x
+        mul!(r,A,dx)
+        r .-= b
+        ldiv!(dx,P,r)
+        x .-= omega .* dx
+        x,step+1
     end
     function finalize!(state)
         (r,dx,P,A_ref) = state
         PartitionedSolvers.finalize!(P)
     end
-    linear_solver(;setup,update!,solve!,finalize!)
+    returns_history = Val(true)
+    linear_solver(;setup,update!,solve!,finalize!,step!,returns_history)
 end
 
 function jacobi(;kwargs...)
@@ -122,7 +139,7 @@ function local_solver_options(A,options)
     end
 end
 
-struct AdditiveSchwarzSetup{A}
+struct AdditiveSchwarzSetup{A} <: AbstractType
     local_setups::A
 end
 

PartitionedSolvers/test/smoothers_tests.jl

@@ -24,6 +24,18 @@ tol = 1.e-8
 update!(S,2*A)
 solve!(y,S,b)
 @test norm(y-x/2)/norm(x/2) < tol
+
+istep = Ref(0)
+for y in iterations!(y,S,b)
+    isa(y,AbstractVector)
+    istep[] += 1
+end
+@test istep[] == 1
+
+y2,hist = solve!(y,S,b;history=true)
+@test y2 === y
+@test hist === nothing
+
 finalize!(S)
 
 solver = linear_solver(LinearAlgebra.lu)
@@ -59,6 +71,18 @@ tol = 1.e-8
 update!(S,2*A)
 solve!(y,S,b)
 @test norm(y-x/2)/norm(x/2) < tol
+
+y,hist = solve!(y,S,b;history=true)
+@test hist.iters == 1000
+
+istep = Ref(0)
+y .= 0
+for y in iterations!(y,S,b)
+    isa(y,AbstractVector)
+    istep[] += 1
+end
+@test norm(y-x/2)/norm(x/2) < tol
+@test istep[] == 1000
 finalize!(S)
 
 solver = additive_schwarz(lu_solver())