High memory allocation when using preconditioners

[mtanneau]

I have noticed some higher-than-expected memory allocations when using a preconditioner other than opEye().

I only tested minres and cg so far, and the culprit seems to be a matrix-vector product with M, namely, in cg and minres.
If I replace, e.g., v = M * r2 in minres by mul!(v, M, r2), the allocations are more reasonable (see example below).

using Krylov, LinearAlgebra, SparseArrays, Random
using BenchmarkTools

Random.seed!(0)
m = 2^10
A = sprand(m, m, 0.01)
S = A+A'
b = ones(m)

M1 = Krylov.opEye()
M2 = 1.0I
M3 = Diagonal(ones(m))

# Warm-up will be done during calibration
@btime minres($S, $b)         # No preconditioner
@btime minres($S, $b, M=$M1)  # M = opEye
@btime minres($S, $b, M=$M2)  # M = I (uniform scaling)
@btime minres($S, $b, M=$M3)  # M = Diagonal(1.0, ..., 1.0)

Current version

  19.906 ms (34 allocations: 114.13 KiB)
  19.902 ms (34 allocations: 114.13 KiB)
  21.170 ms (1059 allocations: 8.24 MiB)
  20.397 ms (54 allocations: 122.78 KiB)

and if I replace with mul!(v, M, r2)

  20.147 ms (34 allocations: 114.13 KiB)
  20.006 ms (34 allocations: 114.13 KiB)
  20.654 ms (35 allocations: 122.25 KiB)
  20.344 ms (54 allocations: 122.78 KiB)

Same behavior with cg.

EDIT: making the above modification to minres causes some tests to fail

[amontoison]

The culprit is here : https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/master/src/variants.jl#L3-L13
I only wrap matrices in a linear operator. I can add typeof(M) <: UniformScaling and typeof(N) <: UniformScaling to solve this problem when M.λ ≠ one(T) or N.λ ≠ one(T). Otherwise if λ = one(T) I can replace M and / or N by opEye().

[mtanneau]

I understand that there's a lot of legacy decisions in play, but... from what I gather, parts of the code rely on the convention that products of the form y = A * x will not allocate.

Currently this is effectively enforced for matrix input, because they are automatically wrapped in a PreallocatedLinearOperator.
The example above shows the limitations of this approach: either you ensure that non-AbstractMatrix objects are wrapped too, or the user must follow the convention that y = A * x will not allocate.
The former puts a lot of maintenance on your side, and the latter is (I think) no longer realistic given the in-place 5-arg mul!.

IMO, it would be more natural to have Krylov methods use mul!, and the user ensure that calls to mul!(y, A, x) will be efficient.

[amontoison]

We should add in the documentation that the user must follow the convention that y = A * x is not allocating for the moment.
For the mul! function, I must update all Krylov methods to support it. I will do it when 3-arg and 5-arg mul! will work as we want in LinearOperators.