Which Krylov package solves A*x=b where x and b can be any "vector-like" object?

IIUC, this is searching for f(X) = b, where b is not vector-like, but rather a Matrix.

Either try to define the missing orthogonalize! method for your output type,
or convert to a regular array (the collect used below) ?

julia> using KrylovKit
julia> b = rand(3, 3)
julia> res, info  = linsolve(b) do x
              return transpose(x) |> collect
       end;
julia> transpose(res) ≈ b
true