Boundserror with multiple right hand sides

[lijas]

I get an error when not allocating x when I have multiple RHS.

A = Symmetric(rand(10,10))
b = rand(10,2)
cg(A,b) #Gives error because x is created as an vector
cg!(zeros(size(b)),A,b) #this works

I guess cg(A,b) should also work?

[lijas]

I talked a bit with @fredrikekre . Iterative solvers might not be suited for multiple RHS since it is unclear which residual to check?

[haampie]

A trivial version would be to just run the exact equivalent of

[cg(A, b[:, 1]), cg(A, b[:, 2]), ..., cg(A, b[:, n])]

with an optimization where matrix multiplications act on a block vector, so that you have a higher flop/memop ratio.

But you need to keep in mind vectors converge at different speeds, so the matrix multiplication should only happen on the block that's not converged (you can split it like [not converged | converged] and re-shuffle at the end).

[lijas]

Do you think there would be an advantage of doing

cg!(x,A,b)

where size(x) = size(b) = (1000, 3) (ie 3 right hand sides)

over just using one rhs at the time:

cg!(x1, A, b[:,1])
cg!(x2, A, b[:,2])
cg!(x2, A, b[:,3])

Lets assume the vectors will converge with roughly the same speed. The first version works atm, but which residual will it use to check for convergence?