QR: handle xtype/dtype returned from LibSuiteSparse that don't match …

…matrix element type (#586) * QR: handle xtype/dtype other than input Tv * Set types by SparseMatrixCSC{Tv, Ti} when possible * One less forced copy --------- Co-authored-by: Viral B. Shah <[email protected]>
JuliaSparse · Nov 26, 2024 · 268d390 · 268d390
1 parent 9731aef
commit 268d390
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 4 deletions.
diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl
@@ -901,6 +901,15 @@ Dense(A::StridedVecOrMatInclAdjAndTrans{T}) where
 
 Dense(A::Sparse) = sparse_to_dense(A)
 
+function Dense(ptr::Ptr{cholmod_dense})
+    if ptr == C_NULL
+        throw(ArgumentError("dense matrix construction failed for " *
+            "unknown reasons. Please submit a bug report."))
+    end
+    s = unsafe_load(ptr)
+    return Dense{jlxtype(s.xtype, s.dtype)}(ptr)
+end
+
 function Base.convert(::Type{Dense{Tnew}}, A::Dense{T}) where {Tnew, T}
     GC.@preserve A begin
         Ap = unsafe_load(pointer(A))

diff --git a/src/solvers/spqr.jl b/src/solvers/spqr.jl
@@ -204,10 +204,10 @@ function LinearAlgebra.qr(A::SparseMatrixCSC{Tv, Ti}; tol=_default_tol(A), order
         C_NULL, C_NULL, C_NULL, C_NULL,
         R, E, H, HPinv, HTau)
 
-    R_ = SparseMatrixCSC(Sparse(R[]))
-    return QRSparse(SparseMatrixCSC(Sparse(H[])),
-                    vec(Array(CHOLMOD.Dense{Tv}(HTau[]))),
-                    SparseMatrixCSC(min(size(A)...),
+    R_ = SparseMatrixCSC{Tv, Ti}(Sparse(R[]))
+    return QRSparse(SparseMatrixCSC{Tv, Ti}(Sparse(H[])),
+                    vec(Array{Tv}(CHOLMOD.Dense(HTau[]))),
+                    SparseMatrixCSC{Tv, Ti}(min(size(A)...),
                                     size(R_, 2),
                                     getcolptr(R_),
                                     rowvals(R_),

diff --git a/test/spqr.jl b/test/spqr.jl
@@ -110,6 +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)
+    A = sparse(eltyA[1 0; 0 1])
+    F = qr(A)
+    @test eltype(F.Q) == eltype(F.R) == eltyA
+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)