From 268d3906e21d69ffd67429f5c052732749e7a2b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=C2=A0Tarasiuk?= <158573010+pawel-tarasiuk-quantumz@users.noreply.github.com> Date: Tue, 26 Nov 2024 20:20:26 +0100 Subject: [PATCH] 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 --- src/solvers/cholmod.jl | 9 +++++++++ src/solvers/spqr.jl | 8 ++++---- test/spqr.jl | 6 ++++++ 3 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 51688ce5..2770bf48 100644 --- 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 index 5f409463..cd1bdf6a 100644 --- 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 index 122741ce..285f3dc9 100644 --- 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)