make \ type stable (#684)

Closes #111 

Since UMFPACK is not going to add support for `Float32` anytime soon I
think the only solution is to manually convert back. Also, since the
code for dealing with the output of `lu` is not generic, there is no
point in converting back the decomposition itself, so I'm just
converting the result of `lu(A) \ b` to the appropriate type.

I also did the analogous change to achieve type stability with
`Float16`. There's also nothing that can be done about CHOLMOD, but for
SPQR the situation is better: the code around it is generic, so I just
converted back `qr(A)` itself.

Author: araujoms

src/linalg.jl

@@ -2252,7 +2252,7 @@ function \(A::AbstractSparseMatrixCSC, B::AbstractVecOrMat)
         if ishermitian(A)
             return \(Hermitian(A), B)
         end
-        return \(lu(A), B)
+        return convert(AbstractArray{typeof(one(eltype(A)) \ one(eltype(B)))}, \(lu(A), B))
     else
         return \(qr(A), B)
     end

src/solvers/cholmod.jl

@@ -1899,14 +1899,15 @@ end
 const RealHermSymComplexHermSSL{Ti, Tr} = Union{
     Symmetric{Tr, SparseMatrixCSC{Tr, Ti}},
     Hermitian{Tr, SparseMatrixCSC{Tr, Ti}},
-    Hermitian{Complex{Tr}, SparseMatrixCSC{Complex{Tr}, Ti}}} where {Ti<:ITypes, Tr<:VRealTypes}
+    Hermitian{Complex{Tr}, SparseMatrixCSC{Complex{Tr}, Ti}}} where {Ti<:ITypes, Tr<:Union{Float64, Float32, Float16}}
 
 function \(A::RealHermSymComplexHermSSL{Ti}, B::StridedVecOrMatInclAdjAndTrans) where {Ti}
+    T = typeof(one(eltype(A)) \ one(eltype(B)))
     F = cholesky(A; check = false)
     if issuccess(F)
-        return \(F, B)
+        return convert(AbstractArray{T}, \(F, B))
     else
-        return \(lu(SparseMatrixCSC{eltype(A), Ti}(A)), B)
+        return convert(AbstractArray{T}, \(lu(SparseMatrixCSC{eltype(A), Ti}(A)), B))
     end
 end
 

src/solvers/spqr.jl

@@ -108,7 +108,7 @@ function _qr!(ordering::Integer, tol::Real, econ::Integer, getCTX::Integer,
     return rnk, _E, _HPinv
 end
 
-struct QRSparseQ{Tv<:CHOLMOD.VTypes,Ti<:Integer} <: AbstractQ{Tv}
+struct QRSparseQ{Tv,Ti<:Integer} <: AbstractQ{Tv}
     factors::SparseMatrixCSC{Tv,Ti}
     τ::Vector{Tv}
     n::Int # Number of columns in original matrix
@@ -134,6 +134,14 @@ struct QRSparse{Tv,Ti} <: LinearAlgebra.Factorization{Tv}
     _ldiv_workspace::Vector{Tv}   # backing storage for work buffer (resizable)
 end
 
+function QRSparse{Tv}(F::QRSparse{<:Number, Ti}) where {Tv, Ti}
+    newfactors = convert(SparseMatrixCSC{Tv}, F.factors)
+    newτ = convert(Vector{Tv}, F.τ)
+    newR = convert(SparseMatrixCSC{Tv}, F.R)
+    newQ = QRSparseQ{Tv,Ti}(newfactors, newτ, size(newR, 2))
+    return QRSparse{Tv,Ti}(newfactors, newτ, newR, newQ, F.cpiv, F.rpivinv, ReentrantLock(), Tv[])
+end
+
 Base.size(F::QRSparse) = (size(F.factors, 1), size(F.R, 2))
 function Base.size(F::QRSparse, i::Integer)
     if i == 1
@@ -167,9 +175,9 @@ solve least squares or underdetermined problems with [`\\`](@ref). The function
 
 !!! note
     `qr(A::SparseMatrixCSC)` uses the SPQR library that is part of [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse).
-    As this library only supports sparse matrices with [`Float64`](@ref) or
-    `ComplexF64` elements, as of Julia v1.4 `qr` converts `A` into a copy that is
-    of type `SparseMatrixCSC{Float64}` or `SparseMatrixCSC{ComplexF64}` as appropriate.
+    As this library only supports sparse matrices with [`Float64`](@ref), `ComplexF64`, `Float32`, or
+    `ComplexF32` elements, calling `qr` on a matrix with a different element type will either convert it to a supported type or
+    raise an error.
 
 # Examples
 ```jldoctest
@@ -230,9 +238,9 @@ function LinearAlgebra.qr(A::SparseMatrixCSC{Tv, Ti}; tol=_default_tol(A), order
                     Tv[])              # _ldiv_workspace (lazily sized on first solve)
 end
 LinearAlgebra.qr(A::SparseMatrixCSC{Float16}; tol=_default_tol(A)) =
-    qr(convert(SparseMatrixCSC{Float32}, A); tol=tol)
+    QRSparse{Float16}(qr(convert(SparseMatrixCSC{Float32}, A); tol=tol))
 LinearAlgebra.qr(A::SparseMatrixCSC{ComplexF16}; tol=_default_tol(A)) =
-    qr(convert(SparseMatrixCSC{ComplexF32}, A); tol=tol)
+    QRSparse{ComplexF16}(qr(convert(SparseMatrixCSC{ComplexF32}, A); tol=tol))
 LinearAlgebra.qr(A::Union{SparseMatrixCSC{T},SparseMatrixCSC{Complex{T}}};
    tol=_default_tol(A)) where {T<:AbstractFloat} =
     throw(ArgumentError(string("matrix type ", typeof(A), "not supported. ",

test/cholmod.jl

@@ -853,9 +853,9 @@ end
               Symmetric(Apre + 10I), Hermitian(Apre + 10I),
               Hermitian(complex(Apre)), Hermitian(complex(Apre) + 10I))
         local A, x, b
-        x = fill(1., 10)
+        x = fill(1, 10)
         b = A*x
-        @test x ≈ A\b
+        @test @inferred A\b ≈ x
         @test transpose(A)\b ≈ A'\b
     end
 end

test/linalg.jl

@@ -1040,4 +1040,14 @@ end
     @test_throws DimensionMismatch D1 * S * D1
 end
 
+@testset "type stability of linear solve" begin
+    for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty})
+        A = sprand(elty, 2, 2, 1.0)
+        B = randn(elty, 2, 2)
+        b = randn(elty, 2)
+        @inferred A \ b
+        @inferred A \ B
+    end
+end
+
 end

test/spqr.jl

@@ -110,18 +110,12 @@ end
     @test (F.Q*F.R)::SparseMatrixCSC == A[F.prow,F.pcol]
 end
 
-@testset "Issue #585 for element type: $eltyA" for eltyA in (Float64, Float32, ComplexF64, ComplexF32)
+@testset "Issue #585 for element type: $eltyA" for eltyA in (Float64, Float32, Float16, ComplexF64, ComplexF32, ComplexF16)
     A = sparse(eltyA[1 0; 0 1])
     F = qr(A)
     @test eltype(F.Q) == eltype(F.R) == eltyA
 end
 
-@testset "Complementing issue #585 for element type: $eltyA" for (eltyA, eltyB) in [(Float16, Float32), (ComplexF16, ComplexF32)]
-    A = sparse(eltyA[1 0; 0 1])
-    F = qr(A)
-    @test eltype(F.Q) == eltype(F.R) == eltyB
-end
-
 @testset "select ordering overdetermined" begin
      A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], randn(nn), m, n)
      b = randn(m)