Preallocate workspace for solvers

[dpo]

Krylov methods

• CG
• CR
• SYMMLQ
• CG-LANCZOS
• CG-LANCZOS-SHIFT-SEQ
• MINRES
• MINRES-QLP
• DQGMRES
• DIOM
• USYMLQ
• USYMQR
• TRICG
• TRIMR
• TRILQR
• CGS
• BICGSTAB
• BILQ
• QMR
• BILQR
• CGLS
• CRLS
• CGNE
• CRMR
• LSLQ
• LSQR
• LSMR
• LNLQ
• CRAIG
• CRAIGMR

For MINRES, this could take the form (not tested)

function minres(A, b :: AbstractVector{T}; kwargs...) where T <: AbstractFloat
  x = similar(b)
  r1 = similar(b)
  window = kwargs.get(:window, 5)
  err_vec = Vector{T}(undef, window)
  minres!(A, b, x, r1, err_vec; kwargs...)
end

function minres!(A, b, x, r1, err_vec;
                M=opEye(), λ :: T=zero(T), atol :: T=√eps(T)/100,
                rtol :: T=√eps(T)/100, etol :: T=√eps(T),
                window :: Int=5, itmax :: Int=0, conlim :: T=1/√eps(T),
                verbose :: Int=0) where T <: AbstractFloat
  ...
  x .= 0
  r1 .= b
  err_vec .= 0
  ...
end

The objective is for minres!() to not allocate at all. We could also define

abstract type KrylovSolver end
mutable struct MinresSolver <: KrylovSolver
  x
  r1
  err_vec
  # maybe default parameter values...
  function MinresSolver(A, b)
    new(...)
  end
end

and pass an instance of MinresSolver around.

@amontoison @geoffroyleconte

[amontoison]

Speak of the devil ! We were talking about it this afternoon with Geoffroy.

[dpo]

😈

[amontoison]

I have something in mind, I will open a pull request when it's ready. It should not break GPU support, which a little bit tricky.

[amontoison]

@geoffroyleconte
I let you try to add this feature for MINRES.
To not break GPU support you should parameterize your solver :

abstract type KrylovSolver{S} end
mutable struct MinresSolver{S} <: KrylovSolver{S}
  x  :: S
  r1 :: S
  ...
end

where S is Vector{Float64}, CuVector{Float64}, ... It should not be Float32, Float64, ...
You can find the variable names of MinresSolver here : https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/master/test/test_alloc.jl#L41

Like Dominique proposed you can create a function like that after :

function minres(A, b :: AbstractVector{T}; kwargs...) where T <: AbstractFloat
  solver = ... # create your MinresSolver 
  minres!(A, b, solver; kwargs...)
end

In minres!, you can add the following lines at the beginning to reuse all the actual code of minres :

x = solver.x
r1 = solver.r1
....

You only need to update the initialization of the variables and add a check to verify that eltype(S) == T.

[amontoison]

fixed by #271