Add an implementation of BLOCK-MINRES

JuliaSmoothOptimizers · Oct 31, 2024 · 3a3c03e · 3a3c03e
1 parent 783ed93
commit 3a3c03e
Show file tree

Hide file tree

Showing 14 changed files with 411 additions and 70 deletions.
diff --git a/docs/src/api.md b/docs/src/api.md
@@ -12,11 +12,10 @@ Krylov.LSLQStats
 Krylov.LsmrStats
 ```
 
-## Solver Types
+## Workspace of Krylov methods
 
 ```@docs
 KrylovSolver
-BlockKrylovSolver
 MinresSolver
 MinaresSolver
 CgSolver
@@ -53,6 +52,13 @@ CraigSolver
 CraigmrSolver
 GpmrSolver
 FgmresSolver
+```
+
+## Workspace of block-Krylov methods
+
+```@docs
+BlockKrylovSolver
+BlockMinresSolver
 BlockGmresSolver
 ```
 

diff --git a/docs/src/block_krylov.md b/docs/src/block_krylov.md
@@ -1,10 +1,7 @@
-## Block-GMRES
-
 !!! note
-    `block_gmres` works on GPUs
-     with Julia 1.11.
+    `block_minres` and `block_gmres` work on GPUs with Julia 1.11.
 
-If you want to use `block_gmres` on previous Julia versions, you can overload the function `Krylov.copy_triangle` with the following code:
+If you want to use `block_minres` and `block_gmres` on previous Julia versions, you can overload the function `Krylov.copy_triangle` with the following code:
 ```julia
 using KernelAbstractions, Krylov
 
@@ -23,6 +20,15 @@ function Krylov.copy_triangle(Q::AbstractMatrix{FC}, R::AbstractMatrix{FC}, k::I
 end
 ```
 
+## Block-MINRES
+
+```@docs
+block_minres
+block_minres!
+```
+
+## Block-GMRES
+
 ```@docs
 block_gmres
 block_gmres!

diff --git a/docs/src/block_processes.md b/docs/src/block_processes.md
@@ -27,6 +27,8 @@ T_{k+1,k} =
 
 The function [`hermitian_lanczos`](@ref hermitian_lanczos(::Any, ::AbstractMatrix{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\Psi_1$ and $T_{k+1,k}$.
 
+Related method: [`BLOCK-MINRES`](@ref block_minres).
+
 ```@docs
 hermitian_lanczos(::Any, ::AbstractMatrix{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)
 ```

diff --git a/docs/src/gpu.md b/docs/src/gpu.md
@@ -107,7 +107,7 @@ if CUDA.functional()
   symmetric = hermitian = true
   opM = LinearOperator(T, n, n, symmetric, hermitian, (y, x) -> ldiv_ic0!(P, x, y, z))
 
-  # Solve an Hermitian positive definite system with an IC(0) preconditioner on GPU
+  # Solve a Hermitian positive definite system with an IC(0) preconditioner on GPU
   x, stats = cg(A_gpu, b_gpu, M=opM)
 end
 ```

diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md
@@ -46,7 +46,7 @@ A Krylov method dedicated to non-Hermitian linear systems allows the three varia
 
 ### Hermitian linear systems
 
-Methods concerned: [`SYMMLQ`](@ref symmlq), [`CG`](@ref cg), [`CG-LANCZOS`](@ref cg_lanczos), [`CG-LANCZOS-SHIFT`](@ref cg_lanczos_shift), [`CR`](@ref cr), [`CAR`](@ref car), [`MINRES`](@ref minres), [`MINRES-QLP`](@ref minres_qlp) and [`MINARES`](@ref minares).
+Methods concerned: [`SYMMLQ`](@ref symmlq), [`CG`](@ref cg), [`CG-LANCZOS`](@ref cg_lanczos), [`CG-LANCZOS-SHIFT`](@ref cg_lanczos_shift), [`CR`](@ref cr), [`CAR`](@ref car), [`MINRES`](@ref minres), [`BLOCK-MINRES`](@ref block_minres), [`MINRES-QLP`](@ref minres_qlp) and [`MINARES`](@ref minares).
 
 When $A$ is Hermitian, we can only use centered preconditioning $L^{-1}AL^{-H}y = L^{-1}b$ with $x = L^{-H}y$.
 Centered preconditioning is a special case of two-sided preconditioning with $P_{\ell} = L = P_r^H$ that maintains hermicity.

diff --git a/src/Krylov.jl b/src/Krylov.jl
@@ -12,6 +12,7 @@ include("block_krylov_utils.jl")
 include("block_krylov_processes.jl")
 include("block_krylov_solvers.jl")
 
+include("block_minres.jl")
 include("block_gmres.jl")
 
 include("cg.jl")

diff --git a/src/block_gmres.jl b/src/block_gmres.jl
@@ -1,6 +1,6 @@
 # An implementation of block-GMRES for the solution of the square linear system AX = B.
 #
-# Alexis Montoison, <alexis.montoison@polymtl.ca>
+# Alexis Montoison, <alexis.montoison@polymtl.ca> -- <amontoison@anl.gov>
 # Argonne National Laboratory -- Chicago, October 2023.
 
 export block_gmres, block_gmres!

diff --git a/src/block_krylov_processes.jl b/src/block_krylov_processes.jl
@@ -3,7 +3,7 @@
 
 #### Input arguments
 
-* `A`: a linear operator that models an Hermitian matrix of dimension `n`;
+* `A`: a linear operator that models a Hermitian matrix of dimension `n`;
 * `B`: a matrix of size `n × p`;
 * `k`: the number of iterations of the block Hermitian Lanczos process.
 

diff --git a/src/block_krylov_solvers.jl b/src/block_krylov_solvers.jl
@@ -1,12 +1,64 @@
 export BlockKrylovSolver
 
-export BlockGmresSolver
+export BlockMinresSolver, BlockGmresSolver
 
-const BLOCK_KRYLOV_SOLVERS = Dict(:block_gmres => :BlockGmresSolver)
+const BLOCK_KRYLOV_SOLVERS = Dict(:block_minres => :BlockMinresSolver,
+                                  :block_gmres  => :BlockGmresSolver )
 
 "Abstract type for using block Krylov solvers in-place"
 abstract type BlockKrylovSolver{T,FC,SV,SM} end
 
+"""
+Type for storing the vectors required by the in-place version of BLOCK-MINRES.
+
+The outer constructors
+
+    solver = BlockMinresSolver(m, n, p, SV, SM)
+    solver = BlockMinresSolver(A, B)
+
+may be used in order to create these vectors.
+`memory` is set to `div(n,p)` if the value given is larger than `div(n,p)`.
+"""
+mutable struct BlockMinresSolver{T,FC,SV,SM} <: BlockKrylovSolver{T,FC,SV,SM}
+  m          :: Int
+  n          :: Int
+  p          :: Int
+  ΔX         :: SM
+  X          :: SM
+  W          :: SM
+  P          :: SM
+  Q          :: SM
+  C          :: SM
+  D          :: SM
+  τ          :: SV
+  warm_start :: Bool
+  stats      :: SimpleStats{T}
+end
+
+function BlockMinresSolver(m, n, p, SV, SM)
+  FC = eltype(SV)
+  T  = real(FC)
+  ΔX = SM(undef, 0, 0)
+  X  = SM(undef, n, p)
+  W  = SM(undef, n, p)
+  P  = SM(undef, 0, 0)
+  Q  = SM(undef, 0, 0)
+  C  = SM(undef, p, p)
+  D  = SM(undef, 2p, p)
+  τ  = SV(undef, p)
+  stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown")
+  solver = BlockMinresSolver{T,FC,SV,SM}(m, n, p, ΔX, X, W, P, Q, C, D, τ, false, stats)
+  return solver
+end
+
+function BlockMinresSolver(A, B)
+  m, n = size(A)
+  s, p = size(B)
+  SM = typeof(B)
+  SV = matrix_to_vector(SM)
+  BlockMinresSolver(m, n, p, SV, SM)
+end
+
 """
 Type for storing the vectors required by the in-place version of BLOCK-GMRES.
 
@@ -71,7 +123,8 @@ function BlockGmresSolver(A, B, memory = 5)
 end
 
 for (KS, fun, nsol, nA, nAt, warm_start) in [
-  (:BlockGmresSolver, :block_gmres!, 1, 1, 0, true)
+  (:BlockMinresSolver, :block_minres!, 1, 1, 0, true)
+  (:BlockGmresSolver , :block_gmres! , 1, 1, 0, true)
 ]
   @eval begin
     size(solver :: $KS) = solver.m, solver.n