Interface Block-GMRES from Krylov.jl

[amontoison]

#552
You probably want to use block_gmres! instead of gmres! by default if the right-hand is an AbstractMatrix.
Note that for block-Krylov methods, the solution is stored in solver.X and not solver.x.

[ChrisRackauckas]

this is misplaced if it's a hard dep.

[ChrisRackauckas]

You probably want to use block_gmres! instead of gmres! by default if the right-hand is an AbstractMatrix.

That's not the right interpretation because it doesn't solve the vec form? So it would get the wrong answer?

[amontoison]

I expect to solve Ax = b if b is a vector or AX = B is B is matrix and not Ax = vec(B).
Otherwise, it's not directly a linear system.
From a linear algebra point of view, it seems wrong to me.

I can understand that your want to apply PDE operators on 2D or 3D matrices but the right-hand side should be already vectorized, like I did in this example:
https://jso.dev/Krylov.jl/dev/matrix_free/#Example-with-discretized-PDE

[ChrisRackauckas]

Otherwise, it's not directly a linear system.

No? In that case B is not an element of the vector space that A acts on, it's a set of b's in the vector space. It's playing fast and loose which is fine when everything is just an Array and gives some computational improvements in a special case but it breaks down when dealing with arbitrary AbstractArray types very fast. If you always interpret the right hand side as an element of the vector space on which A is solving for then you can expect things to actually work more naturally, and correctness always comes first.

[amontoison]

Solving AX = vec(B) and AX = B is completely different, do we agree on that?
The second system can be reformulated as a linear system with a vector as right-hand but the left-hand side is not anymore A but block_diagonal(A, ..., A).
But A, X and B are still linearly related.

It's not the case with AX = vec(B), which is counter intuitive for a package called LinearSolve.jl.

[ChrisRackauckas]

The second system can be reformulated as a linear system with a vector as right-hand but the left-hand side is not anymore A but block_diagonal(A, ..., A).

That's not the A you wrote down. It's really weird to solve a problem for a different linear operator than the one the user gave you.

[amontoison]

You don't solve a problem with a different operator. I just explained how a linear system with a vector as right-hand and a linear system with a matrix as right-hand side are equivalent.

[ChrisRackauckas]

They aren't equivalent, it's overloaded syntax. The latter is actually solving with block_diagonal(A, ..., A). block_diagonal(A, ..., A) != A. It's fancy overloaded syntax that lets it work for vectors and matrices but it breaks down really fast when generalizing.

[amontoison]

I don't understand why you talk about overloading.
The comment is just about linear algebra and what you want to solve.

Solving AX = B where B has p columns is equivalent to solving p systems Ax_i = b_i where b_i is the column i of B and x_i the column i of X.
All these systems can be merged into one:
block_diagonal(A, ..., A) [x_1; x2; ...; x_p] = [b_1; ...; b_p] that we can reformulate as Dy = c.

But the main goal of my messages was to explain that it's completely different than solving AX = vec(B).
A \ B doesn't return the solution of AX = vec(B), so why do you to interprete the problem like this if the right-hand side is a matrix?

Conclusion: block_gmres(A, B) returns the same solution than A \ B if B is a matrix.
gmres(A, vec(B)) doesn't return the solution of AX = B.

[ChrisRackauckas]

This is literally the same thing as taking vector transposes seriously. Yes, obviously it's homomorphic to solving with block_diagonal(A, ..., A). But that's not A, and making it work in the general case requires specifying a system for which you have homomorphisms of linear systems. Take it to a separate issue since this would be a major breaking things and there's many other things that would have to be handled in the type system. Also, it would need more general support as sparse solvers generally do not support this.

[amontoison]

Ok, I understand your issue now 👍