Skip to content

Commit

Permalink
Add an implementation of BLOCK-MINRES
Browse files Browse the repository at this point in the history
  • Loading branch information
amontoison committed Nov 1, 2024
1 parent 9536ef7 commit 29f9b7d
Show file tree
Hide file tree
Showing 19 changed files with 491 additions and 49 deletions.
10 changes: 8 additions & 2 deletions docs/src/api.md
Original file line number Diff line number Diff line change
Expand Up @@ -12,11 +12,10 @@ Krylov.LSLQStats
Krylov.LsmrStats
```

## Solver Types
## Workspace of Krylov methods

```@docs
KrylovSolver
BlockKrylovSolver
MinresSolver
MinaresSolver
CgSolver
Expand Down Expand Up @@ -53,6 +52,13 @@ CraigSolver
CraigmrSolver
GpmrSolver
FgmresSolver
```

## Workspace of block-Krylov methods

```@docs
BlockKrylovSolver
BlockMinresSolver
BlockGmresSolver
```

Expand Down
16 changes: 11 additions & 5 deletions docs/src/block_krylov.md
Original file line number Diff line number Diff line change
@@ -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

Expand All @@ -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!
Expand Down
2 changes: 2 additions & 0 deletions docs/src/block_processes.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
```
Expand Down
2 changes: 1 addition & 1 deletion docs/src/gpu.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
```
Expand Down
2 changes: 1 addition & 1 deletion docs/src/preconditioners.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
1 change: 1 addition & 0 deletions src/Krylov.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
2 changes: 1 addition & 1 deletion src/block_gmres.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# An implementation of block-GMRES for the solution of the square linear system AX = B.
#
# Alexis Montoison, <[email protected]>
# Alexis Montoison, <[email protected]> -- <[email protected]>
# Argonne National Laboratory -- Chicago, October 2023.

export block_gmres, block_gmres!
Expand Down
2 changes: 1 addition & 1 deletion src/block_krylov_processes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
75 changes: 72 additions & 3 deletions src/block_krylov_solvers.jl
Original file line number Diff line number Diff line change
@@ -1,12 +1,80 @@
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
P :: SM
Q :: SM
C :: SM
D :: SM
Φ :: SM
tmp :: SM
Vₖ₋₁ :: SM
Vₖ :: SM
wₖ₋₂ :: SM
wₖ₋₁ :: SM
Hₖ₋₂ :: SM
Hₖ₋₁ :: SM
τₖ₋₂ :: SV
τₖ₋₁ :: 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)
P = SM(undef, 0, 0)
Q = SM(undef, n, p)
C = SM(undef, p, p)
D = SM(undef, 2p, p)
Φ = SM(undef, p, p)
tmp = C isa Matrix ? SM(undef, 0, 0) : SM(undef, p, p)
Vₖ₋₁ = SM(undef, n, p)
Vₖ = SM(undef, n, p)
wₖ₋₂ = SM(undef, n, p)
wₖ₋₁ = SM(undef, n, p)
Hₖ₋₂ = SM(undef, 2p, p)
Hₖ₋₁ = SM(undef, 2p, p)
τₖ₋₂ = SV(undef, 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, P, Q, C, D, Φ, tmp, Vₖ₋₁, Vₖ, wₖ₋₂, wₖ₋₁, Hₖ₋₂, Hₖ₋₁, τₖ₋₂, τₖ₋₁, 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.
Expand Down Expand Up @@ -71,7 +139,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
Expand Down
Loading

0 comments on commit 29f9b7d

Please sign in to comment.