From b8a13ef02219f41969b11d56733155dd985f89b9 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 7 Aug 2024 10:46:39 +0200 Subject: [PATCH 1/2] implement in-place `ldiv!` for Cholesky factorization (#547) * implement solve2 without workspace arguments * implement in-place ldiv! for Cholesky factorization * ldiv! with Dense not required * import ldiv! for tests * safer * fix Dense case and add more tests * rearrange wrap_dense_and_ptr * more tests --- src/solvers/cholmod.jl | 88 +++++++++++++++++++++++++++++++++++++++++- test/cholmod.jl | 44 ++++++++++++++++++--- 2 files changed, 126 insertions(+), 6 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index e155e72f..51688ce5 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -18,7 +18,7 @@ using LinearAlgebra using LinearAlgebra: RealHermSymComplexHerm, AdjOrTrans import LinearAlgebra: (\), AdjointFactorization, cholesky, cholesky!, det, diag, ishermitian, isposdef, - issuccess, issymmetric, ldlt, ldlt!, logdet, + issuccess, issymmetric, ldiv!, ldlt, ldlt!, logdet, lowrankdowndate, lowrankdowndate!, lowrankupdate, lowrankupdate! using SparseArrays @@ -913,6 +913,32 @@ function Base.convert(::Type{Dense{Tnew}}, A::Dense{T}) where {Tnew, T} end Base.convert(::Type{Dense{T}}, A::Dense{T}) where T = A +# Just calling Dense(x) or Dense(b) will allocate new +# `cholmod_dense_struct`s in CHOLMOD. Instead, we want to reuse +# the existing memory. We can do this by creating new +# `cholmod_dense_struct`s and filling them manually. +function wrap_dense_and_ptr(x::StridedVecOrMat{T}) where {T <: VTypes} + dense_x = cholmod_dense_struct() + dense_x.nrow = size(x, 1) + dense_x.ncol = size(x, 2) + dense_x.nzmax = length(x) + dense_x.d = stride(x, 2) + dense_x.x = pointer(x) + dense_x.z = C_NULL + dense_x.xtype = xtyp(eltype(x)) + dense_x.dtype = dtyp(eltype(x)) + return dense_x, pointer_from_objref(dense_x) +end +# We need to use a special handling for the case of `Dense` +# input arrays since the `pointer` refers to the pointer to the +# `cholmod_dense`, not to the array values themselves as for +# standard arrays. +function wrap_dense_and_ptr(x::Dense{T}) where {T <: VTypes} + dense_x_ptr = x.ptr + dense_x = unsafe_load(dense_x_ptr) + return dense_x, pointer_from_objref(dense_x) +end + # This constructor assumes zero based colptr and rowval function Sparse(m::Integer, n::Integer, colptr0::Vector{Ti}, rowval0::Vector{Ti}, @@ -1913,6 +1939,66 @@ const AbstractSparseVecOrMatInclAdjAndTrans = Union{AbstractSparseVecOrMat, AdjO throw(ArgumentError("self-adjoint sparse system solve not implemented for sparse rhs B," * " consider to convert B to a dense array")) +# in-place ldiv! +for TI in IndexTypes + @eval function ldiv!(x::StridedVecOrMat{T}, + L::Factor{T, $TI}, + b::StridedVecOrMat{T}) where {T<:VTypes} + if x === b + throw(ArgumentError("output array must not be aliased with input array")) + end + if size(L, 1) != size(b, 1) + throw(DimensionMismatch("Factorization and RHS should have the same number of rows. " * + "Factorization has $(size(L, 2)) rows, but RHS has $(size(b, 1)) rows.")) + end + if size(L, 2) != size(x, 1) + throw(DimensionMismatch("Factorization and solution should match sizes. " * + "Factorization has $(size(L, 1)) columns, but solution has $(size(x, 1)) rows.")) + end + if size(x, 2) != size(b, 2) + throw(DimensionMismatch("Solution and RHS should have the same number of columns. " * + "Solution has $(size(x, 2)) columns, but RHS has $(size(b, 2)) columns.")) + end + if !issuccess(L) + s = unsafe_load(pointer(L)) + if s.is_ll == 1 + throw(LinearAlgebra.PosDefException(s.minor)) + else + throw(LinearAlgebra.ZeroPivotException(s.minor)) + end + end + + # Just calling Dense(x) or Dense(b) will allocate new + # `cholmod_dense_struct`s in CHOLMOD. Instead, we want to reuse + # the existing memory. We can do this by creating new + # `cholmod_dense_struct`s and filling them manually. + dense_x, dense_x_ptr = wrap_dense_and_ptr(x) + dense_b, dense_b_ptr = wrap_dense_and_ptr(b) + + X_Handle = Ptr{cholmod_dense_struct}(dense_x_ptr) + Y_Handle = Ptr{cholmod_dense_struct}(C_NULL) + E_Handle = Ptr{cholmod_dense_struct}(C_NULL) + status = GC.@preserve x dense_x b dense_b begin + $(cholname(:solve2, TI))( + CHOLMOD_A, L, + Ref(dense_b), C_NULL, + Ref(X_Handle), C_NULL, + Ref(Y_Handle), + Ref(E_Handle), + getcommon($TI)) + end + if Y_Handle != C_NULL + free!(Y_Handle) + end + if E_Handle != C_NULL + free!(E_Handle) + end + @assert !iszero(status) + + return x + end +end + ## Other convenience methods function diag(F::Factor{Tv, Ti}) where {Tv, Ti} f = unsafe_load(typedpointer(F)) diff --git a/test/cholmod.jl b/test/cholmod.jl index 108b0ff8..cbcc0ea9 100644 --- a/test/cholmod.jl +++ b/test/cholmod.jl @@ -13,7 +13,7 @@ using Random using Serialization using LinearAlgebra: I, cholesky, cholesky!, det, diag, eigmax, ishermitian, isposdef, issuccess, - issymmetric, ldlt, ldlt!, logdet, norm, opnorm, Diagonal, Hermitian, Symmetric, + issymmetric, ldiv!, ldlt, ldlt!, logdet, norm, opnorm, Diagonal, Hermitian, Symmetric, PosDefException, ZeroPivotException, RowMaximum using SparseArrays using SparseArrays: getcolptr @@ -138,6 +138,9 @@ Random.seed!(123) @test CHOLMOD.isvalid(chma) @test unsafe_load(pointer(chma)).is_ll == 1 # check that it is in fact an LLt @test chma\b ≈ x + x2 = zero(x) + @inferred ldiv!(x2, chma, b) + @test x2 ≈ x @test nnz(chma) == 489 @test nnz(cholesky(A, perm=1:size(A,1))) > nnz(chma) @test size(chma) == size(A) @@ -281,6 +284,37 @@ end end end +@testset "ldiv! $Tv $Ti" begin + local A, x, x2, b, X, X2, B + A = sprand(10, 10, 0.1) + A = I + A * A' + A = convert(SparseMatrixCSC{Tv,Ti}, A) + factor = cholesky(A) + + x = fill(Tv(1), 10) + b = A * x + x2 = zero(x) + @inferred ldiv!(x2, factor, b) + @test x2 ≈ x + + X = fill(Tv(1), 10, 5) + B = A * X + X2 = zero(X) + @inferred ldiv!(X2, factor, B) + @test X2 ≈ X + + c = fill(Tv(1), size(x, 1) + 1) + C = fill(Tv(1), size(X, 1) + 1, size(X, 2)) + y = fill(Tv(1), size(x, 1) + 1) + Y = fill(Tv(1), size(X, 1) + 1, size(X, 2)) + @test_throws DimensionMismatch ldiv!(y, factor, b) + @test_throws DimensionMismatch ldiv!(Y, factor, B) + @test_throws DimensionMismatch ldiv!(x2, factor, c) + @test_throws DimensionMismatch ldiv!(X2, factor, C) + @test_throws DimensionMismatch ldiv!(X2, factor, b) + @test_throws DimensionMismatch ldiv!(x2, factor, B) +end + end #end for Ti ∈ itypes for Tv ∈ (Float32, Float64) @@ -365,9 +399,9 @@ end @test isa(CHOLMOD.eye(3), CHOLMOD.Dense{Float64}) end -@testset "Core functionality ($elty, $elty2)" for - elty in (Tv, Complex{Tv}), - Tv2 in (Float32, Float64), +@testset "Core functionality ($elty, $elty2)" for + elty in (Tv, Complex{Tv}), + Tv2 in (Float32, Float64), elty2 in (Tv2, Complex{Tv2}), Ti ∈ itypes A1 = sparse(Ti[1:5; 1], Ti[1:5; 2], elty <: Real ? randn(Tv, 6) : complex.(randn(Tv, 6), randn(Tv, 6))) @@ -972,7 +1006,7 @@ end f = ones(size(K, 1)) u = K \ f residual = norm(f - K * u) / norm(f) - @test residual < 1e-6 + @test residual < 1e-6 end @testset "wrapped sparse matrices" begin From 95fd7ff4383cd953b88ec5f7700bfe49f8e7a13f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 7 Aug 2024 11:10:28 +0200 Subject: [PATCH 2/2] Missing space in error message (#554) --- src/solvers/umfpack.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/umfpack.jl b/src/solvers/umfpack.jl index 8baab770..629639e6 100644 --- a/src/solvers/umfpack.jl +++ b/src/solvers/umfpack.jl @@ -397,7 +397,7 @@ lu(A::AbstractSparseMatrixCSC{<:Union{ComplexF16,ComplexF32},Ti}; lu(convert(SparseMatrixCSC{ComplexF64,Ti}, A); check = check) lu(A::Union{AbstractSparseMatrixCSC{T},AbstractSparseMatrixCSC{Complex{T}}}; check::Bool = true) where {T<:AbstractFloat} = - throw(ArgumentError(string("matrix type ", typeof(A), "not supported. ", + throw(ArgumentError(string("matrix type ", typeof(A), " not supported. ", "Try lu(convert(SparseMatrixCSC{Float64/ComplexF64,Int}, A)) for ", "sparse floating point LU using UMFPACK or lu(Array(A)) for generic ", "dense LU.")))