Introduce lu! for UmfpackLU to make use of its symbolic factoriza…

…tion (JuliaLang#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]>

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