Update the API for workspaces (#975)

* Update the API for workspaces

* Fix the tests on GPUs

* More update for the solvers

* Use S::Type

Author: amontoison

docs/src/inplace.md

@@ -15,6 +15,10 @@ XyzSolver(kc::KrylovConstructor)
 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
@@ -24,20 +28,20 @@ See the section [custom workspaces](@ref custom_workspaces) for an example where
 
 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`, `BlockGmresSolver`, `FgmresSolver`, `GpmrSolver`, `CgLanczosShiftSolver` and `CglsLanczosShiftSolver` require an additional argument (`memory` or `nshifts`).
-
 The workspace is always the first argument of the in-place methods:
 
 ```@solvers
-minres_solver = MinresSolver(n, n, Vector{Float64})
+minres_solver = MinresSolver(m, n, Vector{Float64})
 minres!(minres_solver, A1, b1)
 
-dqgmres_solver = DqgmresSolver(n, n, memory, Vector{BigFloat})
-dqgmres!(dqgmres_solver, A2, b2)
+bicgstab_solver = BicgstabSolver(m, n, Vector{ComplexF64})
+bicgstab!(bicgstab_solver, A2, b2)
+
+gmres_solver = GmresSolver(m, n, Vector{BigFloat})
+gmres!(gmres_solver, A3, b3)
 
 lsqr_solver = LsqrSolver(m, n, CuVector{Float32})
-lsqr!(lsqr_solver, A3, b3)
+lsqr!(lsqr_solver, A4, b4)
 ```
 
 A generic function `solve!` is also available and dispatches to the appropriate Krylov method.
@@ -46,6 +50,9 @@ A generic function `solve!` is also available and dispatches to the appropriate
 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.

src/block_krylov_processes.jl

@@ -119,7 +119,7 @@ function nonhermitian_lanczos(A, B::AbstractMatrix{FC}, C::AbstractMatrix{FC}, k
   m, n = size(A)
   t, p = size(B)
   Aᴴ = A'
-  pivoting = VERSION < v"1.9" ? Val{false}() : NoPivot()
+  pivoting = NoPivot()
 
   nnzT = p*p + (k-1)*p*(2*p+1) + div(p*(p+1), 2)
   colptr = zeros(Int, p*k+1)

src/block_krylov_solvers.jl

@@ -17,7 +17,6 @@ The outer constructors
     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
@@ -42,7 +41,7 @@ mutable struct BlockMinresSolver{T,FC,SV,SM} <: BlockKrylovSolver{T,FC,SV,SM}
   stats      :: SimpleStats{T}
 end
 
-function BlockMinresSolver(m, n, p, SV, SM)
+function BlockMinresSolver(m::Integer, n::Integer, p::Integer, SV::Type, SM::Type)
   FC = eltype(SV)
   T  = real(FC)
   ΔX = SM(undef, 0, 0)
@@ -78,12 +77,11 @@ Type for storing the vectors required by the in-place version of BLOCK-GMRES.
 
 The outer constructors
 
-    solver = BlockGmresSolver(m, n, p, memory, SV, SM)
-    solver = BlockGmresSolver(A, B, memory = 5)
+    solver = BlockGmresSolver(m, n, p, SV, SM; memory = 5)
+    solver = BlockGmresSolver(A, B; memory = 5)
 
 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)`.
-`memory` is an optional argument in the second constructor.
 """
 mutable struct BlockGmresSolver{T,FC,SV,SM} <: BlockKrylovSolver{T,FC,SV,SM}
   m          :: Int
@@ -105,7 +103,7 @@ mutable struct BlockGmresSolver{T,FC,SV,SM} <: BlockKrylovSolver{T,FC,SV,SM}
   stats      :: SimpleStats{T}
 end
 
-function BlockGmresSolver(m, n, p, memory, SV, SM)
+function BlockGmresSolver(m::Integer, n::Integer, p::Integer, SV::Type, SM::Type; memory::Integer = 5)
   memory = min(div(n,p), memory)
   FC = eltype(SV)
   T  = real(FC)
@@ -126,12 +124,12 @@ function BlockGmresSolver(m, n, p, memory, SV, SM)
   return solver
 end
 
-function BlockGmresSolver(A, B, memory = 5)
+function BlockGmresSolver(A, B; memory::Integer = 5)
   m, n = size(A)
   s, p = size(B)
   SM = typeof(B)
   SV = matrix_to_vector(SM)
-  BlockGmresSolver(m, n, p, memory, SV, SM)
+  BlockGmresSolver(m, n, p, SV, SM; memory)
 end
 
 for (KS, fun, nsol, nA, nAt, warm_start) in [
@@ -220,11 +218,7 @@ function show(io :: IO, solver :: Union{KrylovSolver{T,FC,S}, BlockKrylovSolver{
     type_i = fieldtype(workspace, i)
     field_i = getfield(solver, name_i)
     size_i = ksizeof(field_i)
-    if (name_i::Symbol in [:w̅, :w̄, :d̅]) && (VERSION < v"1.8.0-DEV")
-      (size_i ≠ 0) && Printf.format(io, format2, string(name_i), type_i, format_bytes(size_i))
-    else
-      (size_i ≠ 0) && Printf.format(io, format, string(name_i), type_i, format_bytes(size_i))
-    end
+    (size_i ≠ 0) && Printf.format(io, format, string(name_i), type_i, format_bytes(size_i))
   end
   @printf(io, "└%s┴%s┴%s┘\n","─"^l1,"─"^l2,"─"^l3)
   if show_stats

src/krylov_solve.jl

@@ -18,7 +18,6 @@
     solve!(solver, args...; kwargs...)
 
 Generic function that dispatches to the appropriate in-place Krylov method based on the type of `solver`.
-This function is not exported to prevent potential conflicts with other Julia packages.
 """
 function solve! end
 
@@ -82,7 +81,7 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs
     @eval begin
       function $(krylov)($(def_args...); memory::Int=20, $(def_kwargs...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
         start_time = time_ns()
-        solver = $workspace(A, b, memory)
+        solver = $workspace(A, b; memory)
         elapsed_time = start_time |> ktimer
         timemax -= elapsed_time
         $(krylov!)(solver, $(args...); $(kwargs...))
@@ -93,7 +92,7 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs
       if !isempty($optargs)
         function $(krylov)($(def_args...), $(def_optargs...); memory::Int=20, $(def_kwargs...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
           start_time = time_ns()
-          solver = $workspace(A, b, memory)
+          solver = $workspace(A, b; memory)
           warm_start!(solver, $(optargs...))
           elapsed_time = start_time |> ktimer
           timemax -= elapsed_time
@@ -188,7 +187,7 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs
     @eval begin
       function $(krylov)($(def_args...); memory :: Int=20, $(def_kwargs...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
         start_time = time_ns()
-        solver = $workspace(A, B, memory)
+        solver = $workspace(A, B; memory)
         elapsed_time = ktimer(start_time)
         timemax -= elapsed_time
         $(krylov!)(solver, $(args...); $(kwargs...))
@@ -199,7 +198,7 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs
       if !isempty($optargs)
         function $(krylov)($(def_args...), $(def_optargs...); memory :: Int=20, $(def_kwargs...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
           start_time = time_ns()
-          solver = $workspace(A, B, memory)
+          solver = $workspace(A, B; memory)
           warm_start!(solver, $(optargs...))
           elapsed_time = ktimer(start_time)
           timemax -= elapsed_time