KrylovJL solvers returning incorrect answer when u0 is a SubArray

[vpuri3]

julia> using LinearSolve, IterativeSolvers

julia> A = rand(4, 4); b = rand(4); u0 = zeros(4);

julia> lp = LinearProblem(A, b; u0 = view(u0, :));

julia> solve(lp, KrylovJL_GMRES()) # broken
u: 4-element view(::Vector{Float64}, :) with eltype Float64:
 0.0
 0.0
 0.0
 0.0

julia> solve(lp, IterativeSolversJL_GMRES())
u: 4-element view(::Vector{Float64}, :) with eltype Float64:
  1.4228823954515826
 -1.3845427042881076
  3.1810595009865623
 -1.512734087737276

julia> solve(lp, LUFactorization())
u: 4-element view(::Vector{Float64}, :) with eltype Float64:
  1.4228823954515823
 -1.3845427042881064
  3.1810595009865623
 -1.512734087737277

[fredrikekre]

This is caused by

LinearSolve.jl/src/iterative_wrappers.jl

Line 220 in 0886d05

solver.x = u

. Because the RHS is ::Vector then the solver struct enforces the same type also for x, and thus u0 will be converted to Vector in the line above. A potential fix is:

diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl
index b471065..6de56a9 100644
--- a/src/iterative_wrappers.jl
+++ b/src/iterative_wrappers.jl
@@ -282,6 +282,11 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...)
         ReturnCode.Success
     end

+    # Copy the solution to the allocated output vector
+    if cache.u !== cache.cacheval.x
+        cache.u .= cache.cacheval.x
+    end
+
     return SciMLBase.build_linear_solution(alg, cache.u, resid, cache;
         iters = stats.niter)
 end