-
Notifications
You must be signed in to change notification settings - Fork 66
Description
From the standpoint of differential geometry, in a problem like Ax = b, x and b are fundamentally different objects: x is a vector but A'*b is a covector. (Specifically, in the associated quadratic form, you have x' * A' * A * x as the contraction of the Hessian--A'*A is a rank-2 co-tensor--and A'*b playing the role of the gradient, which contracts with x but is not "like" x.) Thus in
Krylov.jl/src/krylov_workspaces.jl
Lines 35 to 40 in 7d1b6af
| struct KrylovConstructor{S} | |
| vm::S | |
| vn::S | |
| vm_empty::S | |
| vn_empty::S | |
| end |
vm and vn should be the same type of object.
I ask because I have a highly structured problem for which I've gone to the trouble to create my own AbstractVector and AbstractMatrix types (and yes, they speed up multiplication many-fold compared to using the generic SparseMatrixCSC). Specifically, b = TopBottomVector(btop, bbottom) (conceptually these are concatenated) and btop and bbottom have different types, while x has the same type as bbottom. As I'm sure you would imagine, A also has top and bottom parts. When writing out A as an operator, I split b back into its component parts in order to efficiently implement the multiplications. For solvers like MINRES this works out fine, because the underlying problem is square and the right-hand and left-hand "vectors" can indeed have the same type. But now I'm also implementing support for CGLS where the lack of distinction in CglsWorkspace (which ultimately inherits its type structure from KrylovConstructor) is doom for the strategy I'd prefer.
I'm happy to do the work to fix this and submit a PR, but I want to raise the issue to gauge the enthusiasm (or lack thereof) for what is sure to be a certain amount of work.