Allow vm and vn to be of different types #1038
Conversation
|
We need to update a few things in |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1038 +/- ##
==========================================
+ Coverage 94.68% 97.65% +2.97%
==========================================
Files 45 49 +4
Lines 8027 9512 +1485
==========================================
+ Hits 7600 9289 +1689
+ Misses 427 223 -204 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
We need a few more modifications in I should have done that in the previous breaking-release 😢 |
I am sure you're right, but I'll need some hints about where. Note that @eval krylov_solve!(workspace :: $workspace{T,FC,S}, $(def_args...); $(def_kwargs...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: AbstractVector{FC}} = $(krylov!)(workspace, $(args...); $(kwargs...))seems unlikely to break, although perhaps they could be updated. Is it the subtypes of
Some of those appear to be "chains": for example, AdaptiveRegularization and DCISolver depend on RipQP which does indeed break. But fix RipQP and those two (I didn't check any more) might work fine. What do you think about #1037 (comment)? If that would be acceptable and it succeeds in making it non-breaking, then could this be merged? |
|
If we find a way to make it non-breaking, I am ok to merge it! |
|
Locally this seems non-breaking. |
|
I'm unclear on whether the macos failures are "my fault" (I don't think so). And, can I see the updated breakage report (i.e., like #1038 (comment)) myself? I clicked on the |
I don't think it is your fault. I will restart them now.
You need to go to |
|
Now this is a bit more serious: in particular I added some tests. I also now support 5 methods, CGLS, LSQR, TRICG, TRIMR, and GPMR. I only supported the latter 3 in the generic interface, because there you can specify both (If/when this merges, it should be squashed.) |
timholy
force-pushed
the
teh/split_types
branch
from
98c6c66 to
79b4180
Compare
|
Now I've squashed this down to a single commit and all tests pass locally. JSOSolvers using this branch also passes tests without modification. |
|
@timholy Should I do a new patch release with just these modifications: # TODO: in the next breaking release, change to KrylovWorkspace{T,FC,Sm,Sn} and delete the alias below
abstract type _KrylovWorkspace{T,FC,Sm,Sn} end
"Abstract type for using Krylov solvers in-place."
const KrylovWorkspace{T,FC,S} = _KrylovWorkspace{T,FC,S,S}Such that we can update |
|
If I understand what you're asking (though perhaps I don't), that's already part of this PR. Here's the squashed commit message: The breakage test seems to report no new breakages. |
|
Oh ok, you did something even better than what I expected. I thought that we needed to add It is a very neat solution 🙂 |
|
@gdalle Does it also solve your problem? |
|
For @gdalle's benefit: since I didn't convert diff --git a/src/gmres.jl b/src/gmres.jl
index 5841c90b5..0ecbf7541 100644
--- a/src/gmres.jl
+++ b/src/gmres.jl
@@ -114,7 +114,7 @@ optargs_gmres = (:x0,)
kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream)
@eval begin
- function gmres!(workspace :: GmresWorkspace{T,FC,S}, $(def_args_gmres...); $(def_kwargs_gmres...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: AbstractVector{FC}}
+ function gmres!(workspace :: GmresWorkspace{T,FC,Sm,Sn}, $(def_args_gmres...); $(def_kwargs_gmres...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, Sm <: AbstractVector{FC}, Sn <: AbstractVector{FC}}
# Timer
start_time = time_ns()
@@ -132,12 +132,12 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it
# Check type consistency
eltype(A) == FC || @warn "eltype(A) ≠ $FC. This could lead to errors or additional allocations in operator-vector products."
- ktypeof(b) == S || error("ktypeof(b) must be equal to $S")
+ # ktypeof(b) == Sm || error("ktypeof(b) must be equal to $Sm")
# Set up workspace.
- allocate_if(!MisI , workspace, :q , S, workspace.x) # The length of q is n
- allocate_if(!NisI , workspace, :p , S, workspace.x) # The length of p is n
- allocate_if(restart, workspace, :Δx, S, workspace.x) # The length of Δx is n
+ allocate_if(!MisI , workspace, :q , Sn, workspace.x) # The length of q is n
+ allocate_if(!NisI , workspace, :p , Sn, workspace.x) # The length of p is n
+ allocate_if(restart, workspace, :Δx, Sn, workspace.x) # The length of Δx is n
Δx, x, w, V, z = workspace.Δx, workspace.x, workspace.w, workspace.V, workspace.z
c, s, R, stats = workspace.c, workspace.s, workspace.R, workspace.stats
warm_start = workspace.warm_start
@@ -297,7 +297,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it
# Stopping conditions that do not depend on user input.
# This is to guard against tolerances that are unreasonably small.
resid_decrease_mach = (rNorm + one(T) ≤ one(T))
-
+
# Update stopping criterion.
user_requested_exit = callback(workspace) :: Bool
resid_decrease_lim = rNorm ≤ ε
diff --git a/src/krylov_workspaces.jl b/src/krylov_workspaces.jl
index 34b4be6c2..1b5f0df93 100644
--- a/src/krylov_workspaces.jl
+++ b/src/krylov_workspaces.jl
@@ -2391,15 +2391,15 @@ The following outer constructors can be used to initialize this workspace:
`memory` is set to `n` if the value given is larger than `n`.
"""
-mutable struct GmresWorkspace{T,FC,S} <: KrylovWorkspace{T,FC,S}
+mutable struct GmresWorkspace{T,FC,Sm,Sn} <: KrylovWorkspaceNext{T,FC,Sm,Sn}
m :: Int
n :: Int
- Δx :: S
- x :: S
- w :: S
- p :: S
- q :: S
- V :: Vector{S}
+ Δx :: Sn
+ x :: Sn
+ w :: Sn
+ p :: Sn
+ q :: Sn
+ V :: Vector{Sn}
c :: Vector{T}
s :: Vector{FC}
z :: Vector{FC}
@@ -2409,9 +2409,8 @@ mutable struct GmresWorkspace{T,FC,S} <: KrylovWorkspace{T,FC,S}
stats :: SimpleStats{T}
end
-function GmresWorkspace(kc::KrylovConstructor; memory::Integer = 20)
- S = typeof(kc.vm)
- FC = eltype(S)
+function GmresWorkspace(kc::KrylovConstructor{Sm,Sn}; memory::Integer = 20) where {Sm,Sn}
+ FC = eltype(Sn)
T = real(FC)
m = length(kc.vm)
n = length(kc.vn)
@@ -2421,13 +2420,13 @@ function GmresWorkspace(kc::KrylovConstructor; memory::Integer = 20)
w = similar(kc.vn)
p = similar(kc.vn_empty)
q = similar(kc.vn_empty)
- V = S[similar(kc.vn) for i = 1 : memory]
+ V = Sn[similar(kc.vn) for i = 1 : memory]
c = Vector{T}(undef, memory)
s = Vector{FC}(undef, memory)
z = Vector{FC}(undef, memory)
R = Vector{FC}(undef, div(memory * (memory+1), 2))
stats = SimpleStats(0, false, false, false, 0, T[], T[], T[], 0.0, "unknown")
- workspace = GmresWorkspace{T,FC,S}(m, n, Δx, x, w, p, q, V, c, s, z, R, false, 0, stats)
+ workspace = GmresWorkspace{T,FC,Sm,Sn}(m, n, Δx, x, w, p, q, V, c, s, z, R, false, 0, stats)
return workspace
endWith that: julia> using Krylov, LinearAlgebra, StaticArrays
Precompiling Krylov finished.
1 dependency successfully precompiled in 4 seconds. 2 already precompiled.
julia> A = rand(2, 2)
2×2 Matrix{Float64}:
0.991096 0.709605
0.894476 0.0801216
julia> b = SVector{2}(rand(2))
2-element SVector{2, Float64} with indices SOneTo(2):
0.48823909171743063
0.1536619453559741
julia> workspace = GmresWorkspace(KrylovConstructor(b, rand(2)));
julia> gmres!(workspace, A, b)
GmresWorkspace{Float64, Float64, SVector{2, Float64}, Vector{Float64}}(2, 2, [0.0, 0.0], [0.1259115797700803, 0.5121841870497231], [5.551115123125783e-17, 0.0], [0.0, 0.0], [0.0, 0.0], [[0.9538732844068873, 0.30020952232535425], [-0.30020952232535403, 0.9538732844068875]], [0.941666626182997, 1.0], [0.33654712171274304, 1.45256928199235e-16], [0.27386626227694344, 0.4507589575043843], [1.4531063190880973, 0.1864285293723651, 0.38215837219907683], false, 2, SimpleStats
niter: 2
solved: true
inconsistent: false
indefinite: false
npcCount: 0
residuals: []
Aresiduals: []
κ₂(A): []
timer: 107.95μs
status: solution good enough given atol and rtol
)
julia> Krylov.solution(workspace)
2-element Vector{Float64}:
0.1259115797700803
0.5121841870497231
julia> A \ b
2-element Vector{Float64}:
0.12591157977008033
0.512184187049723
julia> using ComponentArrays
julia> b = ComponentVector(b1=rand(), b2=rand())
ComponentVector{Float64}(b1 = 0.15879398420905388, b2 = 0.013338194763834355)
julia> workspace = GmresWorkspace(KrylovConstructor(b, rand(2)));
julia> gmres!(workspace, A, b);
julia> Krylov.solution(workspace)
2-element Vector{Float64}:
-0.0058668876892307275
0.231972118028271
julia> A \ b
2-element Vector{Float64}:
-0.005866887689230701
0.23197211802827095 |
timholy
force-pushed
the
teh/split_types
branch
from
79b4180 to
d89759d
Compare
|
@amontoison force-push was because I discovered that I needed to update Line 66 in 7d1b6af Otherwise the new workspaces were not being pretty-printed. You may want to play with this a bit before making a new release, just to make sure everything is working as expected. Or, if there is a type that is often used in "general" tests, we could convert that one over in hopes of discovering any additional breakages. |
|
I don't understand why you need If the issue is to provide a right-hand side that is a different type of the |
For problems `Ax=b`, this allows the vector space of `x` (`vn`) to be
different from the vector space of `b` (`vm`); also for some bipartite
problems
[A B; C D] [x; y] = [b; c]
it allows `x` and `y` to be of different types.
The key change is to introduce `KrylovWorkspaceNext{T,FC,Sm,Sn}`, for
which `Sm` and `Sn` are the types of the length-`m` and length-`n`
vectors, respectively. `KrylovWorkspace{T,FC,S}` is now just a type
alias for `KrylovWorkspaceNext{T,FC,S,S}`. This is done for reasons of
backwards compatibility; in the next breaking release, we should rename
`KrylovWorkspaceNext` to `KrylovWorkspace`.
This PR exploits the new flexibility for CGLS, LSQR, TRICG, TRIMR, and
GPMR. Extending to other solvers is left as an exercise for the reader.
Fixes JuliaSmoothOptimizers#1037
Yes, I agree. Do you need these features soon or can you wait the second half of December? I also want to test that I don't break |
timholy
force-pushed
the
teh/split_types
branch
from
d89759d to
d4d1780
Compare
|
Just knowing this is on the schedule for merging soon is sufficient; I can go ahead and use it locally. The package I'm developing this for still has a couple of months of work left before it goes public, and it sounds like this will be merged by then. |
@timholy The feature requested by Guillaume was not related to #766 / #769. He asked me the support of two different types for vectors of size I didn't find a related issue. Maybe we discussed about somewhere else. |
|
I don't have time to review in detail but on the whole it looks like a clever solution to our common problem. It might even be enough to convince me to revert to Krylov.jl inside ImplicitDifferentiation.jl, and to use it for upcoming GPU work! As for GMRES and square systems, I think we had a discussion elsewhere where I said that just because input and output have the same size doesn't mean they need to have the same type. The real key aspect here is whether a given solver is designed assuming that |
I will do a new release before the end of the year. |
Great!
My issue is that you will do I documented that we only need Adding an option to bypass these checks could be a solution for you Guillaume. |
|
You may not need user intervention: you could use an error hint if user code throws an error and perhaps check for In my case, the only two efficient combinations defined are |
|
I kind of agree here: why go through the trouble of throwing an error if the code might work? If it's gonna error, it will error downstream anyway, but being overly careful prevents unexpected uses of your library |
3 participants
This is a proof-of-concept that in aFor problemsAx=b, this allows the vector spaceof
x(vn) to be different from the vector space ofb(vm); also for some bipartite problemsit allows
xandyto be of different types.This PR implements the actual change
only for CGLSfor CGLS, LSQR, TRICG, TRIMR, and GPMR.Extending to other solvers is left as an exercise for the reader.
Fixes #1037