Add a generic interface for Krylov.jl (#983)

Author: amontoison

docs/make.jl

@@ -23,7 +23,8 @@ makedocs(
                                 "Generalized saddle-point and non-Hermitian partitioned systems" => "solvers/gsp.md"],
            "Block-Krylov methods" => "block_krylov.md",
            "In-place methods" => "inplace.md",
-            "Storage requirements" => "storage.md",
+           "Generic interface" => "generic_interface.md",
+           "Storage requirements" => "storage.md",
            "Preconditioners" => "preconditioners.md",
            "GPU support" => "gpu.md",
            "Warm-start" => "warm-start.md",

docs/src/generic_interface.md

@@ -0,0 +1,73 @@
+## [Generic interface](@id generic-interface)
+
+Krylov.jl provides a generic interface for solving linear systems using (block) Krylov methods.
+The interface is designed to be common for all methods and contains three routines [`krylov_workspace`](@ref krylov_workspace), [`krylov_solve`](@ref krylov_solve) and [`krylov_solve!`](@ref krylov_solve!).
+They allow to build Krylov workspaces and call both in-place and out-of-place variants of the solvers with a unified API.
+
+```@docs
+krylov_workspace
+krylov_solve!
+krylov_solve
+```
+
+In-place solvers update the workspace, from which solutions and statistics can be retrieved.
+The following functions are available for post-solve analysis:
+
+```@docs
+Krylov.results
+Krylov.solution
+Krylov.statistics
+Krylov.nsolution
+Krylov.issolved
+```
+
+## Examples
+
+```julia
+using Krylov, Random
+
+# Define a symmetric positive definite matrix A and a right-hand side vector b
+n = 1000
+A = sprandn(n, n, 0.005)
+A = A * A' + I
+b = randn(n)
+
+# Out-of-place interface
+for method in (:cg, :cr, :car)
+    x, stats = krylov_solve(Val{method}(), A, b)
+    r = b - A * x
+    println("Residual norm for $(method): ", norm(r))
+end
+```
+
+```julia
+using Krylov, Random
+
+# Define a square nonsymmetric matrix A and a right-hand side vector b
+n = 100
+A = sprand(n, n, 0.05) + I
+b = rand(n)
+
+# In-place interface
+for method in (:bicgstab, :gmres)
+    # Create a workspace for the Krylov method
+    solver = krylov_workspace(Val(method), A, b)
+
+    # Solve the system in-place
+    krylov_solve!(solver, A, b)
+
+    # Get the statistics
+    stats = statistics(solver)
+
+    # Retrieve the solution
+    x = solution(solver)
+
+    # Check if the solver converged
+    solved = issolved(solver)
+    println("Converged $method: ", solved)
+
+    # Display the number of iterations
+    niter = niterations(solver)
+    println("Number of iterations for $method: ", niter)
+end
+```

docs/src/inplace.md

@@ -3,22 +3,22 @@
 All solvers in Krylov.jl have an in-place variant implemented in a method whose name ends with `!`.
 A workspace (`KrylovSolver`), which contains the storage needed by a Krylov method, can be used to solve multiple linear systems with the same dimensions and the same floating-point precision.
 The section [storage requirements](@ref storage-requirements) specifies the memory needed for each Krylov method.
-Each `KrylovSolver` has three constructors:
+
+Each `KrylovSolver` has three constructors with consistent argument patterns:
 
 ```@constructors
 XyzSolver(A, b)
 XyzSolver(m, n, S)
 XyzSolver(kc::KrylovConstructor)
 ```
+The only exceptions are `CgLanczosShiftSolver` and `CglsLanczosShiftSolver`, which require an additional argument `nshifts`.
+Additionally, some constructors accept keyword arguments.
 
-`Xyz` represents the name of the Krylov method, written in lowercase except for its first letter (e.g., `Cg`, `Minres`, `Lsmr`, `Bicgstab`, etc.).
+`Xyz` represents the name of the Krylov method, written in lowercase except for its first letter (such as `Cg`, `Minres`, `Lsmr` or `Bicgstab`).
 If the name of the Krylov method contains an underscore (e.g., `minres_qlp` or `cgls_lanczos_shift`), the workspace constructor transforms it by capitalizing each word and removing underscores, resulting in names like `MinresQlpSolver` or `CglsLanczosShiftSolver`.
 
 Given an operator `A` and a right-hand side `b`, you can create a `KrylovSolver` based on the size of `A` and the type of `b`, or explicitly provide the dimensions `(m, n)` and the storage type `S`.
 
-!!! note
-    The constructors of `CgLanczosShiftSolver` and `CglsLanczosShiftSolver` require an additional argument `nshifts`.
-
 We assume that `S(undef, 0)`, `S(undef, n)`, and `S(undef, m)` are well-defined for the storage type `S`.
 For more advanced vector types, workspaces can also be created with the help of a `KrylovConstructor`.
 ```@docs
@@ -44,27 +44,6 @@ lsqr_solver = LsqrSolver(m, n, CuVector{Float32})
 lsqr!(lsqr_solver, A4, b4)
 ```
 
-A generic function `solve!` is also available and dispatches to the appropriate Krylov method.
-
-```@docs
-Krylov.solve!
-```
-
-!!! note
-    The function `solve!` is not exported to prevent potential conflicts with other Julia packages.
-
-In-place methods return an updated `solver` workspace.
-Solutions and statistics can be recovered via `solver.x`, `solver.y` and `solver.stats`.
-Functions `solution`, `statistics` and `results` can be also used.
-
-```@docs
-Krylov.nsolution
-Krylov.issolved
-Krylov.solution
-Krylov.statistics
-Krylov.results
-```
-
 ## Examples
 
 We illustrate the use of in-place Krylov solvers with two well-known optimization methods.

src/block_gmres.jl

@@ -7,7 +7,7 @@ export block_gmres, block_gmres!
 
 """
     (X, stats) = block_gmres(A, b::AbstractMatrix{FC};
-                             memory::Int=20, M=I, N=I, ldiv::Bool=false,
+                             memory::Int=5, M=I, N=I, ldiv::Bool=false,
                              restart::Bool=false, reorthogonalization::Bool=false,
                              atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0,
                              timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
@@ -61,6 +61,10 @@ function block_gmres end
 
 where `kwargs` are keyword arguments of [`block_gmres`](@ref).
 
+The keyword argument `memory` is the only exception.
+It is only supported by `block_gmres` and is required to create a `BlockGmresSolver`.
+It cannot be changed later.
+
 See [`BlockGmresSolver`](@ref) for more details about the `solver`.
 """
 function block_gmres! end

src/block_krylov_solvers.jl

@@ -1,11 +1,6 @@
-export BlockKrylovSolver
+export BlockKrylovSolver, BlockMinresSolver, BlockGmresSolver
 
-export BlockMinresSolver, BlockGmresSolver
-
-const BLOCK_KRYLOV_SOLVERS = Dict(:block_minres => :BlockMinresSolver,
-                                  :block_gmres  => :BlockGmresSolver )
-
-"Abstract type for using block Krylov solvers in-place"
+"Abstract type for using block Krylov solvers in-place."
 abstract type BlockKrylovSolver{T,FC,SV,SM} end
 
 """

src/diom.jl

@@ -77,8 +77,9 @@ function diom end
 
 where `kwargs` are keyword arguments of [`diom`](@ref).
 
-Note that the `memory` keyword argument is the only exception.
-It's required to create a `DiomSolver` and can't be changed later.
+The keyword argument `memory` is the only exception.
+It is only supported by `diom` and is required to create a `DiomSolver`.
+It cannot be changed later.
 
 See [`DiomSolver`](@ref) for more details about the `solver`.
 """

src/dqgmres.jl

@@ -77,8 +77,9 @@ function dqgmres end
 
 where `kwargs` are keyword arguments of [`dqgmres`](@ref).
 
-Note that the `memory` keyword argument is the only exception.
-It's required to create a `DqgmresSolver` and can't be changed later.
+The keyword argument `memory` is the only exception.
+It is only supported by [`dqgmres`](@ref) and is required to create a `DqgmresSolver`.
+It cannot be changed later.
 
 See [`DqgmresSolver`](@ref) for more details about the `solver`.
 """

src/fgmres.jl

@@ -79,8 +79,9 @@ function fgmres end
 
 where `kwargs` are keyword arguments of [`fgmres`](@ref).
 
-Note that the `memory` keyword argument is the only exception.
-It's required to create a `FgmresSolver` and can't be changed later.
+The keyword argument `memory` is the only exception.
+It is only supported by [`fgmres`](@ref) and is required to create a `FgmresSolver`.
+It cannot be changed later.
 
 See [`FgmresSolver`](@ref) for more details about the `solver`.
 """

src/fom.jl

@@ -72,8 +72,9 @@ function fom end
 
 where `kwargs` are keyword arguments of [`fom`](@ref).
 
-Note that the `memory` keyword argument is the only exception.
-It's required to create a `FomSolver` and can't be changed later.
+The keyword argument `memory` is the only exception.
+It is only supported by [`fom`](@ref) and is required to create a `FomSolver`.
+It cannot be changed later.
 
 See [`FomSolver`](@ref) for more details about the `solver`.
 """

src/gmres.jl

@@ -72,8 +72,9 @@ function gmres end
 
 where `kwargs` are keyword arguments of [`gmres`](@ref).
 
-Note that the `memory` keyword argument is the only exception.
-It's required to create a `GmresSolver` and can't be changed later.
+The keyword argument `memory` is the only exception.
+It is only supported by [`gmres`](@ref) and is required to create a `GmresSolver`.
+It cannot be changed later.
 
 See [`GmresSolver`](@ref) for more details about the `solver`.
 """