Coming from Matlab, I was expecting the same behavior too.
As far as I can tell, Julia’s \ accepts as right-hand side b only vectors like rand(4) or single-column matrices like rand(4,1). The result is of the same type as b.
The vector form seems to be preferred, not all linear solvers accept a single-column matrix.
using LinearAlgebra, SparseArrays
function solve_multiRHS!(A, b)
lu!(A)
for i = 1:length(b)
b[i] = A\b[i]
end
return b
end
A = rand(4,4)
b = Vector{Float64}[]
for i = 1:3
push!(b, rand(4))
end
solve_multiRHS!(A, b)
The sparse operator, when needed, is then applied to A and the individual vectors of b.