Give sparse \ and qr(A) \ b the eltype of dense \

src/linalg.jl

@@ -554,7 +554,7 @@ function _dot(x::AbstractSparseVector, A::SparseMatrixCSCOrColumnSubset, y::Abst
 end
 ## end of symmetric/Hermitian
 
-\(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:AbstractSparseMatrixCSC}}, B::Vector) = copy(A) \ B
+\(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:AbstractSparseMatrixCSC}}, B::VecOrMat) = copy(A) \ B
 
 ## diagonal scaling
 
@@ -1436,10 +1436,11 @@ end
 _hermitian_solve(A::AbstractSparseMatrixCSC) =
     !(eltype(A) <: Union{Integer, Complex{<:Integer}}) && ishermitian(A)
 
-# `\` and `factorize` choose the same method from the structure of `A`; `f` (or `flu`, for
-# the LU factorization) is applied in each branch so that the result stays inferrable. A
-# Hermitian matrix is only wrapped, because `\` and `factorize` finish it differently.
-function _sparse_factorize(f, A::AbstractSparseMatrixCSC, flu = f)
+# `\` and `factorize` choose the same method from the structure of `A`; `f` (or `fsolver`,
+# for the LU, QR and LQ factorizations) is applied in each branch so that the result stays
+# inferrable, and the QR and LQ factorizations are of `fqr(A)`. A Hermitian matrix is only
+# wrapped, because `\` and `factorize` finish it differently.
+function _sparse_factorize(f, A::AbstractSparseMatrixCSC, fsolver = f, fqr = identity)
     m, n = size(A)
     if m == n
         if istril(A)
@@ -1449,28 +1450,37 @@ function _sparse_factorize(f, A::AbstractSparseMatrixCSC, flu = f)
         elseif _hermitian_solve(A)
             return f(Hermitian(A))
         end
-        return flu(lu(A))
+        return fsolver(lu(A))
     elseif m > n
-        return f(qr(A))
+        return fsolver(qr(fqr(A)))
     else
         # A is wide, so the LQ factorization gives the minimum-norm solution
-        return f(lq(A))
+        return fsolver(lq(fqr(A)))
     end
 end
 
 # the solution is dense, so a diagonal is densified for the solve
 _solve_factor(F) = F
 _solve_factor(D::Diagonal) = Diagonal(Vector(diag(D)))
 _sparse_solve(op, F, B) = \(op(_solve_factor(F)), B)
-# the sparse LU works in double precision, so its solution is converted to the eltype that
-# dense `\` gives
-_sparse_lusolve(op, F, A, B) =
-    convert(AbstractArray{typeof(one(eltype(A)) \ one(eltype(B)))}, _sparse_solve(op, F, B))
+# The eltype that dense `\` gives. It is `nothing` for an abstract eltype, which has no `one`
+# and whose solution is returned as the solver gives it.
+_ldiv_eltype(A, B) = isconcretetype(eltype(A)) && isconcretetype(eltype(B)) ?
+    typeof(one(eltype(A)) \ one(eltype(B))) : nothing
+# The sparse LU works in double precision and the QR and LQ in the precision of the matrix,
+# so their solution is converted to the eltype that dense `\` gives, and for the QR and LQ a
+# matrix of a narrower eltype than that is converted before it is factorized.
+_sparse_solversolve(op, F, A, B) = _convert_solution(_ldiv_eltype(A, B), _sparse_solve(op, F, B))
+_convert_solution(T, X) = convert(AbstractArray{T}, X)
+_convert_solution(::Nothing, X) = X
+_qr_matrix(A, B) = _qr_matrix(A, _ldiv_eltype(A, B))
+_qr_matrix(A, T::Type) = convert(AbstractMatrix{promote_type(float(eltype(A)), real(T))}, A)
+_qr_matrix(A, ::Nothing) = A
 
 function \(A::AbstractSparseMatrixCSC, B::AbstractVecOrMat)
     require_one_based_indexing(A, B)
     return _sparse_factorize(F -> _sparse_solve(identity, F, B), A,
-                             F -> _sparse_lusolve(identity, F, A, B))
+                             F -> _sparse_solversolve(identity, F, A, B), A -> _qr_matrix(A, B))
 end
 for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose))
     @eval function \(xformA::($xformtype){<:Any,<:AbstractSparseMatrixCSC}, B::AbstractVecOrMat)
@@ -1479,12 +1489,12 @@ for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose))
         m, n = size(A)
         if m < n
             # A' is tall, so the least squares solve needs a factorization of A' itself
-            return _sparse_solve(identity, qr($xformop(A)), B)
+            return _sparse_solversolve(identity, qr($xformop(_qr_matrix(A, B))), A, B)
         end
         # Otherwise the transformed choice for A serves. For a wide A' that is the
         # factorization of A, which gives the minimum-norm solution.
         return _sparse_factorize(F -> _sparse_solve($xformop, F, B), A,
-                                 F -> _sparse_lusolve($xformop, F, A, B))
+                                 F -> _sparse_solversolve($xformop, F, A, B), A -> _qr_matrix(A, B))
     end
 end
 

src/solvers/spqr.jl

@@ -654,6 +654,15 @@ LinearAlgebra.rank(F::AdjointQRSparse) = rank(parent(F))
 _rhs_eltype(::Type{T}, ::Type{<:Complex}) where {T<:LinearAlgebra.BlasReal} = Complex{T}
 _rhs_eltype(::Type{T}, ::Type) where {T} = T
 _rhs_array(F, B::AbstractVecOrMat) = convert(Array{_rhs_eltype(eltype(F), eltype(B))}, B)
+# Rather than round the right-hand side to the precision of the factorization, a solve whose
+# solution dense `\` gives a wider eltype converts the factorization to that precision, as
+# dense `\` does. A real factorization stays real for a complex right-hand side.
+_solve_eltype(F, B) = isconcretetype(eltype(B)) ?
+    promote_type(eltype(F), real(typeof(zero(eltype(B)) / oneunit(eltype(F))))) : eltype(F)
+_promoted_factor(F::QRSparse, B) = _promoted_factor(F, _solve_eltype(F, B))
+_promoted_factor(F::QRSparse{T}, ::Type{T}) where {T} = F
+_promoted_factor(F::QRSparse, ::Type{T}) where {T<:LinearAlgebra.BlasFloat} = QRSparse{T}(F)
+_promoted_factor(F::QRSparse, ::Type) = F
 
 ## Two helper methods
 _ret_size(F::Union{QRSparse,AdjointQRSparse}, b::AbstractVector) = (size(F, 2),)
@@ -733,7 +742,7 @@ julia> qr(A)\\fill(1.0, 4)
  0.0
 ```
 """
-(\)(F::QRSparse, B::AbstractVecOrMat) = F\_rhs_array(F, B)
+(\)(F::QRSparse, B::AbstractVecOrMat) = (G = _promoted_factor(F, B); G\_rhs_array(G, B))
 
 function LinearAlgebra.ldiv!(X::StridedVecOrMat{T}, F::QRSparse{T}, B::StridedVecOrMat{T};
                              workspace::Union{Nothing,SpqrWS{T}} = nothing) where {T}
@@ -842,7 +851,8 @@ julia> A'x
  0.0
 ```
 """
-(\)(Fadj::AdjointQRSparse, B::AbstractVecOrMat) = Fadj\_rhs_array(Fadj, B)
+(\)(Fadj::AdjointQRSparse, B::AbstractVecOrMat) =
+    (G = adjoint(_promoted_factor(parent(Fadj), B)); G\_rhs_array(G, B))
 
 function LinearAlgebra.ldiv!(X::StridedVecOrMat{T}, Fadj::AdjointQRSparse{T}, B::StridedVecOrMat{T};
                              workspace::Union{Nothing,SpqrWS{T}} = nothing) where {T}

test/solvers/solvers.jl

@@ -129,7 +129,14 @@ end
     # the adjoint solve converts its result to the eltype the plain solve gives
     x = square' \ Float32.(b)
     @test x isa Vector{Float32} && x ≈ Matrix(square') \ b
-    # only the LU solve converts its result, so a least squares solve takes any right-hand side
+    # the least squares solve gives the eltype of dense `\`, and keeps the precision of `b`
+    x = @inferred Float32.(tall) \ b
+    @test x isa Vector{Float64} && x ≈ Matrix(tall) \ b
+    # so does the solve of a transposed complex Hermitian matrix with a matrix right-hand side
+    Ht = transpose(sparse(ComplexF32[4 im; -im 3]))
+    @test @inferred(Ht \ Float32[1 2; 3 4]) isa Matrix{ComplexF32}
+    # the solution for an abstract eltype is not converted, so a least squares solve takes
+    # any right-hand side
     @test tall \ Any[1.0, 2, 3] ≈ Matrix(tall) \ b
     end
     # when the Cholesky factorization of a Hermitian matrix fails, `factorize` returns

test/solvers/spqr.jl

@@ -254,6 +254,8 @@ end
     b, B = collect(1.0:m), reshape(collect(1.0:3m), 3, m)
     eltyA <: Real && @test F.Q * b ≈ F.Q * eltyA.(b)
     eltyA <: Real && @test B * F.Q' ≈ eltyA.(B) * F.Q'
+    # a solve converts F to the precision of the solution, as dense `\` does
+    @test @inferred(F \ b) isa Vector{widen(eltyA)}
     end
 end