Closed gdalle closed 3 months ago
It probably slipped through the cracks because block_gmres
is not tested, see #842
@gdalle We never defined linear operator * matrix products in LinearOperators.jl. I will have a look 🙂
Thanks! Indeed I had traced it back to https://github.com/JuliaSmoothOptimizers/LinearOperators.jl/issues/322 but I didn't know if this was a design decision or an oversight
Hey @amontoison, gentle ping cause I was curious to see if you had been able to investigate?
I thought it would just require removing some type annotations in the 3-arg and 5-arg mul!
of LinearOperators.jl, but there are some more involved changes with the storage_type
that I am not confident enough to suggest
Hi @gdalle! Sorry, I didn't find time to investigate.
Is it possible to implement your own operator temporary?
You just need to define mul!
, size
and eltype
.
I have an example here: https://github.com/JuliaSmoothOptimizers/KrylovPreconditioners.jl/blob/main/ext/CUDA/operators.jl#L1
@gdalle We never defined linear operator * matrix products in LinearOperators.jl. I will have a look 🙂
We did define *
here: https://github.com/JuliaSmoothOptimizers/LinearOperators.jl/blob/main/src/operations.jl#L132. By convention, we (or I) decided long ago that the result should be an operator. However, nothing prevents us from generalizing
mul!(res::AbstractVector, op::AbstractLinearOperator{T}, v::AbstractVector, α, β)
to matrices.
In Krylov.jl, we only rely on mul!
so we need a constructor that takes a input a function that represent:
mul!(res::AbstractMatrix, op::AbstractLinearOperator{T}, v::AbstractMatrix, α, β)
An issue is that, unlike vectors, matrices can have a different number of columns, and the operators can be dedicated to only certain sizes of matrices. For example, in KrylovPreconditioners.jl, the operators that perform triangular solves (on GPU) only work for a specific number of columns because we preallocate the internal CUSPARSE / rocSPARSE buffer for the associated number of right-hand sides. But a check of the size of the matrices can be done directly in the function provided by the user.
Guillaume, can you give more details on your linear operator?
Sorry for not answering, and thank you for taking an interest ☺️
My linear operator is the JVP/VJP for a vector-to-vector function. It is naturally defined on a single vector, but with AD packages like ForwardDiff and Enzyme, we can also apply "batched" versions which amount to multiplying the Jacobian with several vectors at once, aka a matrix. I can definitely roll out my own mul!
, but won't I lose other things by not using an official LinearOperator?
More generally, how useful is block_gmres if it only applies to actual non-lazy matrices?
I opened a PR in LinearOperators.jl to start working on a fix, would love to get some help @amontoison @dpo!
Now that the LinearOperators PR has been merged, this seems to work:
julia> using Krylov, LinearOperators
julia> A, B = rand(2, 2), rand(2, 2);
julia> block_gmres(A, B)
([0.1341553080088699 0.48826584325663686; 0.9880393854740158 0.3756128570051308], SimpleStats
niter: 1
solved: true
inconsistent: false
residuals: []
Aresiduals: []
κ₂(A): []
timer: 12.10ms
status: solution good enough given atol and rtol
)
julia> block_gmres(LinearOperator(A), B)
([0.1341553080088699 0.48826584325663686; 0.9880393854740158 0.3756128570051308], SimpleStats
niter: 1
solved: true
inconsistent: false
residuals: []
Aresiduals: []
κ₂(A): []
timer: 498.99ms
status: solution good enough given atol and rtol
)
Closing in favor of more specific issues in the future
Did you tested with a Krylov method that requires A' like lsmr
?
I thought block_gmres
was the only one of its kind?
https://jso.dev/Krylov.jl/stable/block_krylov/
Oh, you're right 🥴
I forgot that I didn't implemented a block method that requires A'
yet...
well, get cracking!
I think this is a LinearOperators.jl issue, right?