Add an implementation of BLOCK-GMRES

JuliaSmoothOptimizers · Dec 7, 2023 · 59bee3f · 59bee3f
1 parent a3930dc
commit 59bee3f
Show file tree

Hide file tree

Showing 15 changed files with 648 additions and 134 deletions.
diff --git a/docs/make.jl b/docs/make.jl
@@ -21,6 +21,7 @@ makedocs(
                                 "Adjoint systems" => "solvers/as.md",
                                 "Saddle-point and Hermitian quasi-definite systems" => "solvers/sp_sqd.md",
                                 "Generalized saddle-point and non-Hermitian partitioned systems" => "solvers/gsp.md"],
+           "Block-Krylov methods" => "block_krylov.md",
            "In-place methods" => "inplace.md",
            "Preconditioners" => "preconditioners.md",
            "Storage requirements" => "storage.md",
@@ -39,6 +40,7 @@ makedocs(
                            "TriMR" => "examples/trimr.md",
                            "BICGSTAB" => "examples/bicgstab.md",
                            "DQGMRES" => "examples/dqgmres.md",
+                           "BLOCK-GMRES" => "examples/block_gmres.md",
                            "CGNE" => "examples/cgne.md",
                            "CRMR" => "examples/crmr.md",
                            "CRAIG" => "examples/craig.md",

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -16,6 +16,7 @@ Krylov.LsmrStats
 
 ```@docs
 KrylovSolver
+BlockKrylovSolver
 MinresSolver
 MinaresSolver
 CgSolver
@@ -51,6 +52,7 @@ CraigSolver
 CraigmrSolver
 GpmrSolver
 FgmresSolver
+BlockGmresSolver
 ```
 
 ## Utilities

diff --git a/docs/src/block_krylov.md b/docs/src/block_krylov.md
@@ -0,0 +1,6 @@
+## Block-GMRES
+
+```@docs
+block_gmres
+block_gmres!
+```
diff --git a/docs/src/block_processes.md b/docs/src/block_processes.md
@@ -99,6 +99,8 @@ H_{k+1,k} =
 
 The function [`arnoldi`](@ref arnoldi(::Any, ::AbstractMatrix{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\Gamma$, and $H_{k+1,k}$.
 
+Related method: [`BLOCK-GMRES`](@ref block_gmres).
+
 ```@docs
 arnoldi(::Any, ::AbstractMatrix{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)
 ```

diff --git a/docs/src/examples/block_gmres.md b/docs/src/examples/block_gmres.md
@@ -0,0 +1,19 @@
+```@example block_gmres
+using Krylov, MatrixMarket, SuiteSparseMatrixCollection
+using LinearAlgebra, Printf
+
+ssmc = ssmc_db(verbose=false)
+matrix = ssmc_matrices(ssmc, "HB", "sherman2")
+path = fetch_ssmc(matrix, format="MM")
+
+n = matrix.nrows[1]
+A = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1]).mtx"))
+B = Matrix{Float64}(I, n, 5)
+B_norm = norm(B)
+
+# Solve Ax = B
+X, stats = block_gmres(A, B)
+show(stats)
+R = B - A * X
+@printf("Relative residual: %8.1e\n", norm(R) / B_norm)
+```
diff --git a/docs/src/factorization-free.md b/docs/src/factorization-free.md
@@ -39,7 +39,7 @@ Some methods only require `A * v` products, whereas other ones also require `A'
 |:-----------------------------------------------:|:----------------------------------------:|
 | CG, CR, CAR                                     | CGLS, CRLS, CGNE, CRMR                   |
 | SYMMLQ, CG-LANCZOS, MINRES, MINRES-QLP, MINARES | LSLQ, LSQR, LSMR, LNLQ, CRAIG, CRAIGMR   |
-| DIOM, FOM, DQGMRES, GMRES, FGMRES               | BiLQ, QMR, BiLQR, USYMLQ, USYMQR, TriLQR |
+| DIOM, FOM, DQGMRES, GMRES, FGMRES, BLOCK-GMRES  | BiLQ, QMR, BiLQR, USYMLQ, USYMQR, TriLQR |
 | CGS, BICGSTAB                                   | TriCG, TriMR                             |
 
 !!! info

diff --git a/docs/src/inplace.md b/docs/src/inplace.md
@@ -15,7 +15,7 @@ Given an operator `A` and a right-hand side `b`, you can create a `KrylovSolver`
 For example, use `S = Vector{Float64}` if you want to solve linear systems in double precision on the CPU and `S = CuVector{Float32}` if you want to solve linear systems in single precision on an Nvidia GPU.
 
 !!! note
-    `DiomSolver`, `FomSolver`, `DqgmresSolver`, `GmresSolver`, `FgmresSolver`, `GpmrSolver` and `CgLanczosShiftSolver` require an additional argument (`memory` or `nshifts`).
+    `DiomSolver`, `FomSolver`, `DqgmresSolver`, `GmresSolver`, `BlockGmresSolver`, `FgmresSolver`, `GpmrSolver` and `CgLanczosShiftSolver` require an additional argument (`memory` or `nshifts`).
 
 The workspace is always the first argument of the in-place methods:
 

diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md
@@ -35,7 +35,7 @@ Krylov.jl supports both approaches thanks to the argument `ldiv` of the Krylov s
 
 ### Square non-Hermitian linear systems
 
-Methods concerned: [`CGS`](@ref cgs), [`BiCGSTAB`](@ref bicgstab), [`DQGMRES`](@ref dqgmres), [`GMRES`](@ref gmres), [`FGMRES`](@ref fgmres), [`DIOM`](@ref diom) and [`FOM`](@ref fom).
+Methods concerned: [`CGS`](@ref cgs), [`BiCGSTAB`](@ref bicgstab), [`DQGMRES`](@ref dqgmres), [`GMRES`](@ref gmres), [`BLOCK-GMRES`](@ref block_gmres), [`FGMRES`](@ref fgmres), [`DIOM`](@ref diom) and [`FOM`](@ref fom).
 
 A Krylov method dedicated to non-Hermitian linear systems allows the three variants of preconditioning.
 

diff --git a/src/Krylov.jl b/src/Krylov.jl
@@ -2,13 +2,17 @@ module Krylov
 
 using LinearAlgebra, SparseArrays, Printf
 
-include("krylov_utils.jl")
 include("krylov_stats.jl")
-include("krylov_solvers.jl")
 
-include("processes_utils.jl")
+include("krylov_utils.jl")
 include("krylov_processes.jl")
+include("krylov_solvers.jl")
+
+include("block_krylov_utils.jl")
 include("block_krylov_processes.jl")
+include("block_krylov_solvers.jl")
+
+include("block_gmres.jl")
 
 include("cg.jl")
 include("cr.jl")