Introduce lu! for UmfpackLU to make use of its symbolic factorization (#33738)

* Introduce `lu!` for `UmfpackLU` to make use of its symbolic factorization

* Remove spaces at keyword argument

Co-authored-by: Viral B. Shah <[email protected]>

Author: 2 people

NEWS.md

@@ -43,6 +43,24 @@ Standard library changes
 
 
 #### SparseArrays
+* `lu!` accepts `UmfpackLU` as an argument to make use of its symbolic factorization.
+
+#### Dates
+
+#### Statistics
+
+
+#### Sockets
+
+
+Deprecated or removed
+---------------------
+
+External dependencies
+---------------------
+
+Tooling Improvements
+---------------------
 
 
 <!--- generated by NEWS-update.jl: -->

stdlib/SuiteSparse/src/umfpack.jl

@@ -6,7 +6,7 @@ export UmfpackLU
 
 import Base: (\), getproperty, show, size
 using LinearAlgebra
-import LinearAlgebra: Factorization, det, lu, ldiv!
+import LinearAlgebra: Factorization, det, lu, lu!, ldiv!
 
 using SparseArrays
 using SparseArrays: getcolptr
@@ -176,6 +176,65 @@ lu(A::Union{SparseMatrixCSC{T},SparseMatrixCSC{Complex{T}}};
     "dense LU.")))
 lu(A::SparseMatrixCSC; check::Bool = true) = lu(float(A); check = check)
 
+"""
+    lu!(F::UmfpackLU, A::SparseMatrixCSC; check=true) -> F::UmfpackLU
+
+Compute the LU factorization of a sparse matrix `A`, reusing the symbolic
+factorization of an already existing LU factorization stored in `F`. The
+sparse matrix `A` must have an identical nonzero pattern as the matrix used
+to create the LU factorization `F`, otherwise an error is thrown.
+
+When `check = true`, an error is thrown if the decomposition fails.
+When `check = false`, responsibility for checking the decomposition's
+validity (via [`issuccess`](@ref)) lies with the user.
+
+!!! note
+    `lu!(F::UmfpackLU, A::SparseMatrixCSC)` uses the UMFPACK library that is part of
+    SuiteSparse. As this library only supports sparse matrices with [`Float64`](@ref) or
+    `ComplexF64` elements, `lu!` converts `A` into a copy that is of type
+    `SparseMatrixCSC{Float64}` or `SparseMatrixCSC{ComplexF64}` as appropriate.
+
+!!! compat "Julia 1.4"
+    `lu!` for `UmfpackLU` requires at least Julia 1.4.
+
+# Examples
+```jldoctest
+julia> A = sparse(Float64[1.0 2.0; 0.0 3.0]);
+
+julia> F = lu(A);
+
+julia> B = sparse(Float64[1.0 1.0; 0.0 1.0]);
+
+julia> lu!(F, B);
+
+julia> F \\ ones(2)
+2-element Array{Float64,1}:
+ 0.0
+ 1.0
+```
+"""
+function lu!(F::UmfpackLU, S::SparseMatrixCSC{<:UMFVTypes,<:UMFITypes}; check::Bool=true)
+    zerobased = getcolptr(S)[1] == 0
+    F.m = size(S, 1)
+    F.n = size(S, 2)
+    F.colptr = zerobased ? copy(getcolptr(S)) : decrement(getcolptr(S))
+    F.rowval = zerobased ? copy(rowvals(S)) : decrement(rowvals(S))
+    F.nzval = copy(nonzeros(S))
+
+    umfpack_numeric!(F)
+    check && (issuccess(F) || throw(LinearAlgebra.SingularException(0)))
+    return F
+end
+lu!(F::UmfpackLU, A::SparseMatrixCSC{<:Union{Float16,Float32},Ti};
+   check::Bool = true) where {Ti<:UMFITypes} =
+    lu!(F, convert(SparseMatrixCSC{Float64,Ti}, A); check = check)
+lu!(F::UmfpackLU, A::SparseMatrixCSC{<:Union{ComplexF16,ComplexF32},Ti};
+   check::Bool = true) where {Ti<:UMFITypes} =
+    lu!(F, convert(SparseMatrixCSC{ComplexF64,Ti}, A); check = check)
+lu!(F::UmfpackLU, A::Union{SparseMatrixCSC{T},SparseMatrixCSC{Complex{T}}};
+   check::Bool = true) where {T<:AbstractFloat} =
+    throw(ArgumentError(string("matrix type ", typeof(A), "not supported.")))
+lu!(F::UmfpackLU, A::SparseMatrixCSC; check::Bool = true) = lu!(F, float(A); check = check)
 
 size(F::UmfpackLU) = (F.m, F.n)
 function size(F::UmfpackLU, dim::Integer)
@@ -262,35 +321,33 @@ for itype in UmfpackIndexTypes
             return U
         end
         function umfpack_numeric!(U::UmfpackLU{Float64,$itype})
-            if U.numeric != C_NULL return U end
             if U.symbolic == C_NULL umfpack_symbolic!(U) end
             tmp = Vector{Ptr{Cvoid}}(undef, 1)
             status = ccall(($num_r, :libumfpack), $itype,
                            (Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Cvoid},
                             Ptr{Float64}, Ptr{Float64}),
                            U.colptr, U.rowval, U.nzval, U.symbolic, tmp,
                            umf_ctrl, umf_info)
+            U.status = status
             if status != UMFPACK_WARNING_singular_matrix
                 umferror(status)
             end
             U.numeric = tmp[1]
-            U.status = status
             return U
         end
         function umfpack_numeric!(U::UmfpackLU{ComplexF64,$itype})
-            if U.numeric != C_NULL return U end
             if U.symbolic == C_NULL umfpack_symbolic!(U) end
             tmp = Vector{Ptr{Cvoid}}(undef, 1)
             status = ccall(($num_c, :libumfpack), $itype,
                            (Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Cvoid},
                             Ptr{Float64}, Ptr{Float64}),
                            U.colptr, U.rowval, real(U.nzval), imag(U.nzval), U.symbolic, tmp,
                            umf_ctrl, umf_info)
+            U.status = status
             if status != UMFPACK_WARNING_singular_matrix
                 umferror(status)
             end
             U.numeric = tmp[1]
-            U.status = status
             return U
         end
         function solve!(x::StridedVector{Float64}, lu::UmfpackLU{Float64,$itype}, b::StridedVector{Float64}, typ::Integer)

stdlib/SuiteSparse/test/umfpack.jl

@@ -189,6 +189,34 @@ using LinearAlgebra: Adjoint, Transpose, SingularException
         end
     end
 
+    @testset "Reuse symbolic LU factorization" begin
+        A1 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
+                    increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
+                    [2.,1.,3.,4.,-1.,-3.,3.,9.,2.,1.,4.,2.], 5, 5)
+        for Tv in (Float64, ComplexF64, Float16, Float32, ComplexF16, ComplexF32)
+            for Ti in Base.uniontypes(SuiteSparse.UMFPACK.UMFITypes)
+                A = convert(SparseMatrixCSC{Tv,Ti}, A0)
+                B = convert(SparseMatrixCSC{Tv,Ti}, A1)
+                b = Tv[8., 45., -3., 3., 19.]
+                F = lu(A)
+                lu!(F, B)
+                @test F\b ≈ B\b ≈ Matrix(B)\b
+
+                # singular matrix
+                C = copy(B)
+                C[4, 3] = Tv(0)
+                F = lu(A)
+                @test_throws SingularException lu!(F, C)
+
+                # change of nonzero pattern
+                D = copy(B)
+                D[5, 1] = Tv(1.0)
+                F = lu(A)
+                @test_throws ArgumentError lu!(F, D)
+            end
+        end
+    end
+
 end
 
 @testset "REPL printing of UmfpackLU" begin