From 260e88d62e2910b297a717bba10867826a16d6ef Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Sat, 31 Oct 2026 10:58:03 -0400 Subject: [PATCH] Require Julia 1.12 and fix memory-safety and failure-path problems Bump the version to 0.7.0, as several of these change behaviour. Julia 1.12 no longer sets up SuiteSparse when SparseArrays loads, so call init_suitesparse from _common, as its own solvers do on entry, rather than from __init__. Memory safety: - A factorization shares vectors with the caller, so dropzeros!(A) or resize! could leave K.nzval too short, and libklu read past its end. Operations that read these vectors now check the lengths first. Sharing stays, to avoid using more memory, and the docs now say what the caller must not do. - The pattern check in klu! shifted colptr and rowval in place, which other factorizations sharing them could see while inside libklu. - klu_analyze! validated the caller's permutations and then used them, so another task could change them in between. It now validates and uses copies. - Assigning the current libklu pointer freed it and left it dangling, and replacing the symbolic pointer left numeric factors that depended on it. - Fields that don't need to change are const, which leaves _symbolic, _numeric, _perms and nzval. Memory leak: - deserialize does not call the constructor, so a deserialized factorization had no finalizer, and never freed the libklu objects it made on demand. Register the finalizer when the first of them is attached instead. Failure paths: - When klu! stopped partway it left a mix of old and new factors, and a failed klu_factor! kept the earlier factors, so K \ b returned finite wrong answers. Both now release them, so later use throws until a factorization succeeds. - rcond and condest keep the status, as they could make issuccess return true. - An earlier failure no longer stops klu_factor! from retrying. - solve! with check=false doesn't throw when it factors on demand. - Check the result of klu_sort. - klu rejects empty matrices with a clear error, and KLU_INVALID has a clearer message. - SingularException reports a column with no pivot. Locking: - issuccess, nnz and show read state under the lock, and show prints factors from one factorization. Also: - klu_refactor! accepts keywords - add logabsdet, and so det and logdet, including for K' and transpose(K) - ldiv! accepts a single row with any column stride - K.symbolic and K.numeric, which hold raw pointers, are no longer listed as public - the docs for the scale setting describe rows rather than columns - remove the AdjointFact and TransposeFact aliases, and the unused `@isok` macro Co-Authored-By: Claude Opus 5.5 --- Project.toml | 12 +- README.md | 11 +- docs/Project.toml | 2 +- docs/src/index.md | 12 +- src/KLU.jl | 292 ++++++++++++++++++++++++++++---------------- test/runtests.jl | 301 +++++++++++++++++++++++++++++++++++++++++++++- 6 files changed, 515 insertions(+), 115 deletions(-) diff --git a/Project.toml b/Project.toml index 9fe71c5..c2103e1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "KLU" uuid = "ef3ab10e-7fda-4108-b977-705223b18434" authors = ["Raye Kimmerer and contributors"] -version = "0.6.0" +version = "0.7.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -9,14 +9,16 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SuiteSparse_jll = "bea87d4a-7f5b-5778-9afe-8cc45184846c" [compat] -LinearAlgebra = "1.10" +LinearAlgebra = "1.12" +Serialization = "1.11" SuiteSparse_jll = "7.2" Test = "1.10" -SparseArrays = "1.10" -julia = "1.10" +SparseArrays = "1.12" +julia = "1.12" [extras] +Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Serialization", "Test"] diff --git a/README.md b/README.md index a088a57..5e735d7 100644 --- a/README.md +++ b/README.md @@ -48,6 +48,14 @@ julia> klu!(factor, B); # B has the same pattern as A julia> klu!(factor, nonzeros(B)); # or just its values ``` +If `klu!` throws a `SingularException`, `factor` has no numeric factorization until a later +`klu!` or `klu_factor!` succeeds, and solving with it throws as well. + +To avoid using more memory, `factor` shares the vectors it is given: after `klu(A)` or +`klu!(factor, B)` with `Float64` or `ComplexF64` values, `nonzeros(factor)` is the vector of `A` +or `B`. It is the caller's responsibility not to change them in a way that breaks `factor`, as +described in the docstring of `klu`. + `klu(A; full_factor = false)` only analyzes the pattern of `A`, and `klu!` then computes the first numeric factorization. Use `copy(factor)` to get an independent factorization, for example to refactor it without affecting the original. This costs as much as a new @@ -56,7 +64,8 @@ factorization. ## Checking a factorization `rcond` is a cheap estimate of the reciprocal condition number, `condest` is a more accurate -estimate of the 1-norm condition number, and `rgrowth` is the reciprocal pivot growth. +estimate of the 1-norm condition number, and `rgrowth` is the reciprocal pivot growth. `det`, +`logabsdet`, and `logdet` give the determinant of the factored matrix. ## Thread safety diff --git a/docs/Project.toml b/docs/Project.toml index 14da2a4..bf35cbd 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,4 +3,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" [compat] Documenter = "1" -julia = "1.10" +julia = "1.12" diff --git a/docs/src/index.md b/docs/src/index.md index 402e627..c13bb77 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -49,6 +49,14 @@ julia> klu!(factor, B); # B has the same pattern as A julia> klu!(factor, nonzeros(B)); # or just its values ``` +If `klu!` throws a `SingularException`, `factor` has no numeric factorization until a later +`klu!` or `klu_factor!` succeeds, and solving with it throws as well. + +To avoid using more memory, `factor` shares the vectors it is given: after `klu(A)` or +`klu!(factor, B)` with `Float64` or `ComplexF64` values, `nonzeros(factor)` is the vector of `A` +or `B`. It is the caller's responsibility not to change them in a way that breaks `factor`, as +described in the docstring of [`klu`](@ref). + `klu(A; full_factor = false)` only analyzes the pattern of `A`, and `klu!` then computes the first numeric factorization. Use `copy(factor)` to get an independent factorization, for example to refactor it without affecting the original. This costs as much as a new @@ -57,7 +65,8 @@ factorization. ## Checking a factorization `rcond` is a cheap estimate of the reciprocal condition number, `condest` is a more accurate -estimate of the 1-norm condition number, and `rgrowth` is the reciprocal pivot growth. +estimate of the 1-norm condition number, and `rgrowth` is the reciprocal pivot growth. `det`, +`logabsdet`, and `logdet` give the determinant of the factored matrix. ## Thread safety @@ -91,6 +100,7 @@ solve! rcond condest rgrowth +KLU.LinearAlgebra.logabsdet(::KLU.KLUFactorization) KLU.KLUFactorization Base.copy(::KLU.KLUFactorization) ``` diff --git a/src/KLU.jl b/src/KLU.jl index 449177e..844c06d 100644 --- a/src/KLU.jl +++ b/src/KLU.jl @@ -29,9 +29,7 @@ end increment(A::AbstractArray{<:Integer}) = increment!(copy(A)) using LinearAlgebra - -const AdjointFact = isdefined(LinearAlgebra, :AdjointFactorization) ? LinearAlgebra.AdjointFactorization : Adjoint -const TransposeFact = isdefined(LinearAlgebra, :TransposeFactorization) ? LinearAlgebra.TransposeFactorization : Transpose +using LinearAlgebra: AdjointFactorization, TransposeFactorization const KLUTypes = Union{Float64, ComplexF64} const KLUValueTypes = (:Float64, :ComplexF64) @@ -52,15 +50,19 @@ function kluerror(status::Integer) elseif status == KLU_OUT_OF_MEMORY throw(OutOfMemoryError()) elseif status == KLU_INVALID - throw(ArgumentError("Invalid Status")) + throw(ArgumentError("libklu reported invalid input (KLU_INVALID), such as a malformed sparsity pattern or invalid settings in `common`.")) elseif status == KLU_TOO_LARGE - throw(OverflowError("Integer overflow has occured")) + throw(OverflowError("libklu reported an integer overflow (KLU_TOO_LARGE). Try Int64 indices.")) else throw(ErrorException("Unknown KLU error code: $status")) end end -kluerror(common::Union{klu_l_common, klu_common}) = kluerror(common.status) +function kluerror(common::Union{klu_l_common, klu_common}) + # `singular_col` is the zero-based index of a column for which libklu found no pivot. + common.status == KLU_SINGULAR && throw(LinearAlgebra.SingularException(common.singular_col + 1)) + return kluerror(common.status) +end """ Data structure for parameters of and statistics generated by KLU functions. @@ -70,8 +72,8 @@ Data structure for parameters of and statistics generated by KLU functions. - `btf::Int64`: If `btf != 0` use BTF pre-ordering - `ordering::Int64`: If `ordering == 0` use AMD to permute, if `ordering == 1` use COLAMD, \ if `ordering == 2` use the permutations given to [`klu_analyze!`](@ref), if `ordering == 3` use the user provided ordering function. -- `scale::Int64`: If `scale == 1` then `A[:,i] ./= sum(abs.(A[:,i]))`, if `scale == 2` then \ -`A[:,i] ./= maximum(abs.(A[:,i]))`. If `scale == 0` no scaling is done, and the input is \ +- `scale::Int64`: If `scale == 1` then `A[i,:] ./= sum(abs.(A[i,:]))`, if `scale == 2` then \ +`A[i,:] ./= maximum(abs.(A[i,:]))`. If `scale == 0` no scaling is done, and the input is \ checked for errors if `scale >= 0`. See the [KLU User Guide](https://github.com/DrTimothyAldenDavis/SuiteSparse/raw/dev/KLU/Doc/KLU_UserGuide.pdf) @@ -86,15 +88,15 @@ This is the `Int32` version of [`klu_l_common`](@ref). """ klu_common -macro isok(A) - :(kluerror($(esc(A)))) -end - function _klu_name(name, Tv, Ti) outname = "klu_" * (Tv === :Float64 ? "" : "z") * (Ti === :Int64 ? "l_" : "_") * name return Symbol(replace(outname, "__"=>"_")) end function _common(T) + # SparseArrays requires this before any use of libklu. Otherwise its solvers switch the + # allocator on first use, and memory from libc's malloc is then freed through `jl_free`. + # Each KLUFactorization gets its `common` here, before libklu allocates anything for it. + SparseArrays.LibSuiteSparse.init_suitesparse() if T == Int64 common = klu_l_common() ok = klu_l_defaults(Ref(common)) @@ -119,16 +121,19 @@ abstract type AbstractKLUFactorization{Tv, Ti} <: LinearAlgebra.Factorization{Tv const _EXITING = Threads.Atomic{Bool}(false) function __init__() atexit(() -> (_EXITING[] = true; nothing)) - # SparseArrays requires this before any use of libklu. Otherwise its solvers switch the - # allocator on first use, and memory from libc's malloc is then freed through `jl_free`. - # Julia before 1.12 does this when SparseArrays loads. - @static if isdefined(SparseArrays.LibSuiteSparse, :init_suitesparse) - SparseArrays.LibSuiteSparse.init_suitesparse() - end return end +# libklu only sees pointers, so it would read past the end of vectors that are too short. +# It also rejects n <= 0 itself, but only with a generic KLU_INVALID. +function _check_lengths(n, colptr, rowval, nzval) + n >= 1 || throw(ArgumentError("KLU cannot factor an empty matrix: n must be positive, but it is $n.")) + length(colptr) == n + 1 && colptr[end] <= min(length(rowval), length(nzval)) || + throw(ArgumentError("colptr must have length n + 1, and rowval and nzval at least colptr[end] entries.")) + return nothing +end + """ KLUFactorization <: Factorization @@ -154,37 +159,29 @@ corrupt each other's results. Distinct `KLUFactorization` objects share no state used from different threads simultaneously without contention. """ mutable struct KLUFactorization{Tv<:KLUTypes, Ti<:KLUITypes} <: AbstractKLUFactorization{Tv, Ti} - common::Union{klu_l_common, klu_common} + const common::Union{klu_l_common, klu_common} + # Only these, `nzval` and `_perms` stay mutable, as they change when `K` is analyzed or factored. _symbolic::Ptr{Cvoid} _numeric::Ptr{Cvoid} - n::Int - colptr::Vector{Ti} - rowval::Vector{Ti} + const n::Int + const colptr::Vector{Ti} + const rowval::Vector{Ti} nzval::Vector{Tv} # Serializes all libklu calls on this object. KLU keeps per-call scratch space inside # the Numeric object, so concurrent solves/refactors on one factorization are a data race. - lock::ReentrantLock + const lock::ReentrantLock # The permutations given to `klu_analyze!`, if any, so that a copy can analyze with them. _perms::Union{Nothing, Tuple{Vector{Ti}, Vector{Ti}}} - function KLUFactorization(n, colptr, rowval, nzval) - # libklu only sees pointers, so it would read past the end of vectors that are too short. - n >= 0 && length(colptr) == n + 1 && colptr[end] <= min(length(rowval), length(nzval)) || - throw(ArgumentError("colptr must have length n + 1, and rowval and nzval at least colptr[end] entries.")) + function KLUFactorization(n, colptr, rowval, nzval, common = _common(eltype(colptr))) + _check_lengths(n, colptr, rowval, nzval) Ti = eltype(colptr) - common = _common(Ti) - obj = new{eltype(nzval), Ti}(common, C_NULL, C_NULL, n, colptr, rowval, nzval, ReentrantLock(), nothing) - # During normal operation an object being finalized is unreachable, so no other task - # can be holding its lock or be inside a libklu call on it. (Any in-flight call keeps - # the object rooted.) That does not hold at process exit, see `_EXITING`. - function f(klu) - _EXITING[] && return - _free_symbolic(klu) - _free_numeric(klu) - end - return finalizer(f, obj) + return new{eltype(nzval), Ti}(common, C_NULL, C_NULL, n, colptr, rowval, nzval, ReentrantLock(), nothing) end end +# `K` shares its vectors with the caller, who could have resized them since they were checked. +_check_lengths(K::KLUFactorization, nzval = K.nzval) = _check_lengths(K.n, K.colptr, K.rowval, nzval) + """ lock(K::KLUFactorization) unlock(K::KLUFactorization) @@ -243,31 +240,46 @@ function size(K::AbstractKLUFactorization, dim::Integer) end end -nnz(K::AbstractKLUFactorization) = K.lnz + K.unz + K.nzoff -"""The nonzeros of the factorization `K`. +nnz(K::AbstractKLUFactorization) = Base.@lock K K.lnz + K.unz + K.nzoff +""" + nonzeros(K::KLUFactorization) -> Vector + +The values of the matrix factored by `K`. + !!! warning "nonzeros(K) may or may not point to nonzeros(A)" - After `K = klu(A)`, it may or may not be the case that `nonzeros(K) === nonzeros(A)`. \ - If `A` is of type `SparseMatrixCSC{Float32, Int64}`, then `nonzeros(K)` will be \ - a `Vector{Float64}`, while `nonzeros(A)` is a `Vector{Float32}`, so they're definitely distinct. \ - On the other hand, if `A` has value type `Float64`, then changing the values in `nonzeros(A)` \ - will change `nonzeros(K)` as well. + After `K = klu(A)`, `nonzeros(K) === nonzeros(A)` if `A` has `Float64` or `ComplexF64` values, + so changing one changes the other. Other value types are converted to a new vector. + See [`klu`](@ref) for what changing these values means for `K`. """ nonzeros(K::AbstractKLUFactorization) = K.nzval -if !isdefined(LinearAlgebra, :AdjointFactorization) # VERSION < v"1.10-" - Base.adjoint(K::AbstractKLUFactorization) = Adjoint(K) +Base.transpose(K::AbstractKLUFactorization) = TransposeFactorization(K) + +# During normal operation an object being finalized is unreachable, so no other task +# can be holding its lock or be inside a libklu call on it. (Any in-flight call keeps +# the object rooted.) That does not hold at process exit, see `_EXITING`. +function _release(klu::AbstractKLUFactorization) + _EXITING[] && return + _free_symbolic(klu) + _free_numeric(klu) end -Base.transpose(K::AbstractKLUFactorization) = TransposeFact(K) # Assigning a new symbolic/numeric pointer releases the one it replaces, so that # re-running `klu_factor!` on an already factored object does not leak the old factor. +# Numeric factors belong to their symbolic analysis, so they are released with it. function setproperty!(klu::AbstractKLUFactorization, s::Symbol, x) - if s === :_symbolic + new = convert(fieldtype(typeof(klu), s), x) + if s === :_symbolic && new != getfield(klu, s) + _free_numeric(klu) _free_symbolic(klu) - elseif s === :_numeric + # Registered here rather than in the constructor, which `deserialize` does not call. + # Numeric factors need a symbolic analysis, so this covers them, and analyzing again + # only adds another call which finds nothing left to free. + new != C_NULL && finalizer(_release, klu) + elseif s === :_numeric && new != getfield(klu, s) _free_numeric(klu) end - return setfield!(klu, s, convert(fieldtype(typeof(klu), s), x)) + return setfield!(klu, s, new) end # Certain sets of inputs must be non-null *together*: @@ -290,7 +302,7 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes Base.@lock klu begin klu._numeric == C_NULL && throw(ArgumentError("This KLUFactorization has not yet been factored. Try `klu_factor!`.")) status = klu.common.status - $sort(klu._symbolic, klu._numeric, Ref(klu.common)) + $sort(klu._symbolic, klu._numeric, Ref(klu.common)) == 1 || kluerror(klu.common) ok = $call if ok == 1 klu.common.status = status # libklu resets it, which would hide KLU_SINGULAR @@ -304,8 +316,9 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes end function Base.propertynames(::AbstractKLUFactorization, private::Bool=false) - publicnames = (:lnz, :unz, :nzoff, :L, :U, :F, :q, :p, :Rs, :symbolic, :numeric,) - privatenames = (:nblocks, :maxblock, :lock, :(_L), :(_U), :(_F)) + publicnames = (:lnz, :unz, :nzoff, :L, :U, :F, :q, :p, :Rs,) + # `symbolic` and `numeric` are copies of libklu's structs, with pointers into its memory. + privatenames = (:symbolic, :numeric, :nblocks, :maxblock, :lock, :(_L), :(_U), :(_F)) if private return (publicnames..., privatenames...) else @@ -316,11 +329,9 @@ end function getproperty(klu::AbstractKLUFactorization{Tv, Ti}, s::Symbol) where {Tv<:KLUTypes, Ti<:KLUITypes} # Forwards to the numeric struct: if s ∈ (:lnz, :unz, :nzoff) - klu._numeric == C_NULL && throw(ArgumentError("This KLUFactorization has not yet been factored. Try `klu_factor!`.")) return getproperty(klu.numeric, s) end if s ∈ (:nblocks, :maxblock) - klu._symbolic == C_NULL && throw(ArgumentError("This KLUFactorization has not yet been analyzed. Try `klu_analyze!`.")) return getproperty(klu.symbolic, s) end # libklu has no getters, so these copy its structs. They are the only reads of libklu's memory. @@ -442,22 +453,23 @@ function getproperty(klu::AbstractKLUFactorization{Tv, Ti}, s::Symbol) where {Tv end function LinearAlgebra.issuccess(K::AbstractKLUFactorization; allowsingular=false) - return (allowsingular ? K.common.status >= KLU_OK : K.common.status == KLU_OK) && K._numeric != C_NULL + return Base.@lock K (allowsingular ? K.common.status >= KLU_OK : K.common.status == KLU_OK) && K._numeric != C_NULL end function show(io::IO, mime::MIME{Symbol("text/plain")}, K::AbstractKLUFactorization) summary(io, K); println(io) - if K._numeric != C_NULL + # Extract all the factors first, so that they are from one factorization, and printing to + # a slow `io` does not hold up other tasks. + factors = Base.@lock K K._numeric == C_NULL ? nothing : (K.L, K.U, K.F) + if factors === nothing + println(io, "Incomplete Factorization, please try `klu_factor!(K)`.") + else + L, U, F = factors println(io, "L factor:") - show(io, mime, K.L) + show(io, mime, L) println(io, "\nU factor:") - show(io, mime, K.U) - F = K.F - if F !== nothing - println(io, "\nF factor:") - show(io, mime, K.F) - end - else - println(io, "Incomplete Factorization, please try `klu_factor!(K)`.") + show(io, mime, U) + println(io, "\nF factor:") + show(io, mime, F) end end @@ -469,6 +481,7 @@ represented by the fields of `K`.""" function klu_analyze!(K::KLUFactorization{Tv, Ti}; check=true) where {Tv, Ti<:KLUITypes} Base.@lock K begin if K._symbolic != C_NULL return K end + _check_lengths(K) if Ti == Int64 sym = klu_l_analyze(K.n, K.colptr, K.rowval, Ref(K.common)) else @@ -492,10 +505,13 @@ which are zero-based, unlike `K.p` and `K.q`. Like `klu_analyze!`, this does not has already been analyzed.""" function klu_analyze!(K::KLUFactorization{Tv, Ti}, P::Vector{Ti}, Q::Vector{Ti}; check=true) where {Tv, Ti<:KLUITypes} # libklu does not validate these, and reads out of bounds if they are not permutations. + # Validate copies, as another task could change the originals after the check. + P, Q = copy(P), copy(Q) all(x -> length(x) == K.n && isperm(x .+ one(Ti)), (P, Q)) || throw(ArgumentError("P and Q must be zero-based permutations of length $(K.n).")) Base.@lock K begin if K._symbolic != C_NULL return K end + _check_lengths(K) if Ti == Int64 sym = klu_l_analyze_given(K.n, K.colptr, K.rowval, P, Q, Ref(K.common)) else @@ -505,7 +521,7 @@ function klu_analyze!(K::KLUFactorization{Tv, Ti}, P::Vector{Ti}, Q::Vector{Ti}; kluerror(K.common) else K._symbolic = sym - K._perms = (copy(P), copy(Q)) + K._perms = (P, Q) end end return K @@ -516,14 +532,17 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes @eval begin function klu_factor!(K::KLUFactorization{$Tv, $Ti}; check=true, allowsingular=false) Base.@lock K begin - K._symbolic == C_NULL && K.common.status >= KLU_OK && klu_analyze!(K; check) - if K._symbolic != C_NULL && K.common.status >= KLU_OK + _check_lengths(K) + K._symbolic == C_NULL && klu_analyze!(K; check) + if K._symbolic != C_NULL K.common.halt_if_singular = !allowsingular && check num = $factor(K.colptr, K.rowval, K.nzval, K._symbolic, Ref(K.common)) K.common.halt_if_singular = true else num = C_NULL end + # Also on failure: the factors of the previous values must not stay in use. + K._numeric = num if num == C_NULL && check kluerror(K.common) else @@ -533,7 +552,6 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes (K.common.status == KLU_OK) || (check && kluerror(K.common)) end end - K._numeric = num end return K end @@ -548,6 +566,7 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes function rgrowth(K::KLUFactorization{$Tv, $Ti}) Base.@lock K begin K._numeric == C_NULL && klu_factor!(K) + _check_lengths(K) status = K.common.status ok = $rgrowth(K.colptr, K.rowval, K.nzval, K._symbolic, K._numeric, Ref(K.common)) if ok == 0 @@ -562,10 +581,12 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes function rcond(K::AbstractKLUFactorization{$Tv, $Ti}) Base.@lock K begin K._numeric == C_NULL && klu_factor!(K) + status = K.common.status ok = $rcond(K._symbolic, K._numeric, Ref(K.common)) if ok == 0 kluerror(K.common) else + K.common.status = status # libklu resets it, which would hide KLU_SINGULAR return K.common.rcond end end @@ -574,10 +595,13 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes function condest(K::KLUFactorization{$Tv, $Ti}) Base.@lock K begin K._numeric == C_NULL && klu_factor!(K) + _check_lengths(K) + status = K.common.status ok = $condest(K.colptr, K.nzval, K._symbolic, K._numeric, Ref(K.common)) if ok == 0 kluerror(K.common) else + K.common.status = status # libklu resets it, which would hide KLU_SINGULAR return K.common.condest end end @@ -607,6 +631,52 @@ Accurately estimate the 1-norm condition number of the factorization. It is `Inf """ condest +# The sign of a permutation, from the number of its cycles of even length. +function _permsign(p::AbstractVector{<:Integer}) + seen = falses(length(p)) + s = 1 + for i in eachindex(p) + seen[i] && continue + len = 0 + j = i + while !seen[j] + seen[j] = true + j = Int(p[j]) + len += 1 + end + iseven(len) && (s = -s) + end + return s +end + +""" + logabsdet(K::KLUFactorization) -> (logabsdet, sign) + det(K::KLUFactorization) + +The log of the absolute value of the determinant of the factored matrix, and its sign, which +has unit magnitude, or is 0 for a singular matrix. `det` is computed from these, and `logdet` is +also available. `K` is factored if it has not been yet. +""" +function LinearAlgebra.logabsdet(K::KLUFactorization{Tv}) where {Tv} + # `R \\ A[p, q] = L * U + F`, where L has a unit diagonal and `L * U + F` is block upper + # triangular, with `F` strictly above the diagonal blocks. + d, Rs, p, q = Base.@lock K begin + K._numeric == C_NULL && klu_factor!(K) + diag(K.U), K.Rs, K.p, K.q + end + s = Tv(_permsign(p) * _permsign(q)) + logdet = sum(log, Rs; init = 0.0) + for x in d + s *= sign(x) + logdet += log(abs(x)) + end + return logdet, s +end + +LinearAlgebra.logabsdet(K::AdjointFactorization{<:Any, <:KLUFactorization}) = + ((l, s) = logabsdet(parent(K)); (l, conj(s))) +LinearAlgebra.logabsdet(K::TransposeFactorization{<:Any, <:KLUFactorization}) = logabsdet(parent(K)) + """ klu_factor!(K::KLUFactorization; check=true, allowsingular=false) -> K::KLUFactorization @@ -658,6 +728,7 @@ The relation between `K` and `A` is - `\\`, `ldiv!`, and [`solve!`](@ref) to solve linear systems - [`klu!`](@ref) to refactor with new values - [`rcond`](@ref), [`condest`](@ref), and [`rgrowth`](@ref) to check the quality of the factorization +- `det`, `logabsdet`, and `logdet` - `copy` and `deepcopy` # Arguments @@ -671,8 +742,20 @@ The relation between `K` and `A` is !!! note `klu(A::SparseMatrixCSC)` uses the KLU[^ACM907] library that is part of SuiteSparse. As this library only supports sparse matrices with `Float64` or - `ComplexF64` elements, `klu` converts `A` into a copy that is of type - `SparseMatrixCSC{Float64}` or `SparseMatrixCSC{ComplexF64}` as appropriate. + `ComplexF64` elements, `klu` converts other values to one of these, in a new vector. + +!!! warning "`K` shares vectors with the caller" + To avoid using more memory, `K` does not copy what it is given. `klu(A)` keeps `nonzeros(A)`, + `klu(n, colptr, rowval, nzval)` keeps all three vectors, and `klu!(K, nzval)` keeps `nzval`, + unless the values had to be converted. It is the caller's responsibility not to change + them in a way that breaks `K`: + + - Changing values does not change the current factors, but later calls to + [`klu_factor!`](@ref), [`condest`](@ref), [`rgrowth`](@ref), and `copy` read the new values. + - Resizing a vector, for example with `dropzeros!(A)`, makes the operations that read it throw + an `ArgumentError` if it became too short. + - `colptr` and `rowval` must not be changed at all while `K` is in use, as libklu relies on + them matching its analysis. This is not checked. !!! note "Thread safety" A single `KLUFactorization` may be shared between tasks/threads: every libklu call on @@ -738,6 +821,10 @@ Recompute the KLU factorization `K` using the values of `A` or `nzval`. The patt matrix used to create `K` must match the pattern of `A`. If `K` has not been factored yet, as after `klu(A; full_factor = false)`, this computes its first numeric factorization. +If this throws because the new values are singular, `K` has no numeric factorization until a +later `klu!` or [`klu_factor!`](@ref) succeeds. Solving with `K` in the meantime factors the +singular values again, and so throws as well. + For sparse `A` with real or complex element type, the return type of `K` is `KLUFactorization{Tv, Ti}`, with `Tv` = `Float64` or `ComplexF64` respectively and `Ti` is an integer type (`Int32` or `Int64`). @@ -754,8 +841,8 @@ See also: [`klu`](@ref) !!! note `klu(A::SparseMatrixCSC)` uses the KLU[^ACM907] library that is part of SuiteSparse. As this library only supports sparse matrices with `Float64` or - `ComplexF64` elements, `klu` converts `A` into a copy that is of type - `SparseMatrixCSC{Float64}` or `SparseMatrixCSC{ComplexF64}` as appropriate. + `ComplexF64` elements, `klu!` converts other values to one of these, in a new vector. + Otherwise `K` keeps the vector it is given, see [`klu`](@ref) for what this means. [^ACM907]: Davis, Timothy A., & Palamadai Natarajan, E. (2010). Algorithm 907: KLU, A Direct Sparse Solver for Circuit Simulation Problems. ACM Trans. Math. Softw., 37(3). doi:10.1145/1824801.1824814 """ @@ -767,6 +854,7 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes function klu!(K::KLUFactorization{$Tv, $Ti}, nzval::Vector{$Tv}; check=true, allowsingular=false) length(nzval) != length(K.nzval) && throw(DimensionMismatch()) Base.@lock K begin + _check_lengths(K, nzval) K.nzval = nzval K._numeric == C_NULL && return klu_factor!(K; check, allowsingular) K.common.halt_if_singular = !allowsingular && check @@ -775,6 +863,8 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes if (ok == 1 || !check || (allowsingular && K.common.status >= KLU_OK)) return K else + # libklu stopped partway, leaving a mix of the old and new factors. + K._numeric = C_NULL kluerror(K.common) end end @@ -790,29 +880,26 @@ function klu!(K::AbstractKLUFactorization{Float64}, nzval::Vector{U}; check=true return klu!(K, convert(Vector{Float64}, nzval); check, allowsingular) end +# Compares without shifting the zero-based vector in place, as it may be shared with the +# caller and with other factorizations, which could be inside libklu at the time. +function _same_indices(zero_based::Vector, one_based::Vector) + length(zero_based) == length(one_based) || return false + return all(i -> zero_based[i] + 1 == one_based[i], eachindex(zero_based, one_based)) +end + function klu!(K::KLUFactorization{U}, S::SparseMatrixCSC{U}; check=true, allowsingular=false) where {U} size(K) == size(S) || throw(ArgumentError("Sizes of K and S must match.")) - # The pattern check below temporarily shifts K.colptr/K.rowval to 1-based indices, - # so it must not interleave with any libklu call on K. - Base.@lock K begin - increment!(K.colptr) - increment!(K.rowval) - # what should happen here when check = false? This is not really a KLU error code. - K.colptr == S.colptr && K.rowval == S.rowval || - (decrement!(K.colptr); decrement!(K.rowval); - throw(ArgumentError("The pattern of the original matrix must match the pattern of the refactor.")) - ) - decrement!(K.colptr) - decrement!(K.rowval) - return klu!(K, S.nzval; check, allowsingular) - end + # what should happen here when check = false? This is not really a KLU error code. + _same_indices(K.colptr, S.colptr) && _same_indices(K.rowval, S.rowval) || + throw(ArgumentError("The pattern of the original matrix must match the pattern of the refactor.")) + return klu!(K, S.nzval; check, allowsingular) end """ klu_refactor!(...) -> K::KLUFactorization Same as [`klu!`](@ref): alias for naming consistency reasons.""" -klu_refactor!(args...) = klu!(args...) +klu_refactor!(args...; kwargs...) = klu!(args...; kwargs...) """ copy(K::KLUFactorization) -> KLUFactorization @@ -842,10 +929,10 @@ function Base.deepcopy_internal(K::KLUFactorization{Tv, Ti}, stackdict::IdDict) C = KLUFactorization(K.n, Base.deepcopy_internal(K.colptr, stackdict), Base.deepcopy_internal(K.rowval, stackdict), - Base.deepcopy_internal(K.nzval, stackdict)) + Base.deepcopy_internal(K.nzval, stackdict), + Base.deepcopy_internal(K.common, stackdict)) stackdict[K] = C - C.common = Base.deepcopy_internal(K.common, stackdict) - C.common.status = KLU_OK # klu_factor! refuses to proceed after an earlier failure + C.common.status = KLU_OK # so that the outcome of the steps below can be told apart C.common.memusage = C.common.mempeak = 0 # these count the libklu objects of `C` only if K._symbolic != C_NULL K._perms === nothing ? klu_analyze!(C) : klu_analyze!(C, K._perms...) @@ -862,8 +949,9 @@ end # libklu's `ldim`: the distance between the starts of columns of `B`, in elements of `B`. libklu # only steps forward by it, so it must be at least the column length. Unused for a single column. function _ldim(B::StridedVecOrMat) - stride(B, 1) == 1 || throw(ArgumentError("B must have unit stride along its first dimension")) n, nrhs = size(B, 1), size(B, 2) + # libklu never steps along a column with a single entry. + n <= 1 || stride(B, 1) == 1 || throw(ArgumentError("B must have unit stride along its first dimension")) nrhs <= 1 && return n ldim = stride(B, 2) ldim >= n || throw(ArgumentError("The columns of B must be disjoint and in increasing memory order, but B has $n rows and a column stride of $ldim.")) @@ -876,7 +964,7 @@ end # with check=false it *is technically* possible to enter a seriously invalid state wherein `klu.common` # is undefined. I don't think this is possible in practice though, since we throw on invalid `common` on construction. """ - solve!(klu::Union{KLUFactorization, AdjointFact{Tv, KLUFactorization}, TransposeFact{Tv, KLUFactorization}}, B::StridedVecOrMat) -> B + solve!(klu::Union{KLUFactorization, AdjointFactorization{Tv, KLUFactorization}, TransposeFactorization{Tv, KLUFactorization}}, B::StridedVecOrMat) -> B Solve the linear system `AX = B` or `A'X = B` using the matrix factorization `klu` and right-hand side `B`. @@ -910,7 +998,7 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes ldim = _ldim(B) size(B, 1) == size(klu, 1) || throw(DimensionMismatch()) Base.@lock klu begin - klu._numeric == C_NULL && klu_factor!(klu) + klu._numeric == C_NULL && klu_factor!(klu; check) status = klu.common.status isok = $solve_name(klu._symbolic, klu._numeric, ldim, size(B, 2), B, Ref(klu.common)) isok == 0 ? (check && kluerror(klu.common)) : (klu.common.status = status) @@ -928,26 +1016,26 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes call = :($tsolve(klu._symbolic, klu._numeric, ldim, size(B, 2), B, Ref(klu.common))) end @eval begin - function solve!(klu::AdjointFact{$Tv, K}, B::StridedVecOrMat{$Tv}; check=true) where {K<:AbstractKLUFactorization{$Tv, $Ti}} + function solve!(klu::AdjointFactorization{$Tv, K}, B::StridedVecOrMat{$Tv}; check=true) where {K<:AbstractKLUFactorization{$Tv, $Ti}} conj = 1 klu = parent(klu) ldim = _ldim(B) size(B, 1) == size(klu, 1) || throw(DimensionMismatch()) Base.@lock klu begin - klu._numeric == C_NULL && klu_factor!(klu) + klu._numeric == C_NULL && klu_factor!(klu; check) status = klu.common.status isok = $call isok == 0 ? (check && kluerror(klu.common)) : (klu.common.status = status) end return B end - function solve!(klu::TransposeFact{$Tv, K}, B::StridedVecOrMat{$Tv}; check=true) where {K<: AbstractKLUFactorization{$Tv, $Ti}} + function solve!(klu::TransposeFactorization{$Tv, K}, B::StridedVecOrMat{$Tv}; check=true) where {K<: AbstractKLUFactorization{$Tv, $Ti}} conj = 0 klu = parent(klu) ldim = _ldim(B) size(B, 1) == size(klu, 1) || throw(DimensionMismatch()) Base.@lock klu begin - klu._numeric == C_NULL && klu_factor!(klu) + klu._numeric == C_NULL && klu_factor!(klu; check) status = klu.common.status isok = $call isok == 0 ? (check && kluerror(klu.common)) : (klu.common.status = status) @@ -965,7 +1053,7 @@ end # we are not adding check to `ldiv!` since it is not in the contract. LinearAlgebra.ldiv!(klu::AbstractKLUFactorization{Tv}, B::StridedVecOrMat{Tv}) where {Tv<:KLUTypes} = solve!(klu, B) -LinearAlgebra.ldiv!(klu::Union{AdjointFact{Tv, K},TransposeFact{Tv, K}}, B::StridedVecOrMat{Tv}) where {Tv, Ti, K<:AbstractKLUFactorization{Tv, Ti}} = +LinearAlgebra.ldiv!(klu::Union{AdjointFactorization{Tv, K},TransposeFactorization{Tv, K}}, B::StridedVecOrMat{Tv}) where {Tv, Ti, K<:AbstractKLUFactorization{Tv, Ti}} = solve!(klu, B) function LinearAlgebra.ldiv!(klu::AbstractKLUFactorization{<:AbstractFloat}, B::StridedVecOrMat{<:Complex}) imagX = solve(klu, imag(B)) @@ -973,7 +1061,7 @@ function LinearAlgebra.ldiv!(klu::AbstractKLUFactorization{<:AbstractFloat}, B:: map!(complex, B, realX, imagX) end -function LinearAlgebra.ldiv!(klu::Union{AdjointFact{Tv, K},TransposeFact{Tv, K}}, B::StridedVecOrMat{<:Complex}) where {Tv<:AbstractFloat, Ti, K<:AbstractKLUFactorization{Tv, Ti}} +function LinearAlgebra.ldiv!(klu::Union{AdjointFactorization{Tv, K},TransposeFactorization{Tv, K}}, B::StridedVecOrMat{<:Complex}) where {Tv<:AbstractFloat, Ti, K<:AbstractKLUFactorization{Tv, Ti}} imagX = solve(klu, imag(B)) realX = solve(klu, real(B)) map!(complex, B, realX, imagX) diff --git a/test/runtests.jl b/test/runtests.jl index 904c9ec..6cc799e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,8 @@ using KLU using KLU: increment!, KLUITypes, decrement, klu!, KLUFactorization, klu_analyze!, klu_factor!, klu_refactor! using Test -using SparseArrays: SparseMatrixCSC, sparse, sprand, nnz, nonzeros +using Serialization +using SparseArrays: SparseMatrixCSC, sparse, sprand, spzeros, dropzeros!, nnz, nonzeros using LinearAlgebra @testset "KLU Wrappers" begin @@ -95,10 +96,10 @@ using LinearAlgebra klu_factor!(K) @test K.numeric.lnz == 7 == K.numeric.unz @test K.nblocks == 3 - @test propertynames(K) == (:lnz, :unz, :nzoff, :L, :U, :F, :q, :p, :Rs, :symbolic, :numeric,) + @test propertynames(K) == (:lnz, :unz, :nzoff, :L, :U, :F, :q, :p, :Rs,) @test_throws ArgumentError size(K, -1) @test K.R[[1, end]] == [1, 6] && length(K.R) == K.nblocks + 1 - @test (:nblocks, :lock) ⊆ propertynames(K, true) + @test (:symbolic, :numeric, :nblocks, :lock) ⊆ propertynames(K, true) R = convert(SparseMatrixCSC{Tv, Ti}, sparse(ones(2, 3))) @test_throws DimensionMismatch klu(R) @test_throws ArgumentError KLUFactorization(R) @@ -113,7 +114,7 @@ using LinearAlgebra @test klu!(F, B.nzval)\b ≈ B\b ≈ Matrix(B)\b #test just supplying nzval for recompute @test klu_refactor!(F, A)\b ≈ Matrix(A)\b - # The pattern check temporarily shifts the indices of F, which must be undone. + # A rejected pattern leaves F usable. P = copy(A) P[findfirst(iszero, Matrix(A))] = one(Tv) @test_throws ArgumentError klu!(F, P) @@ -129,6 +130,26 @@ using LinearAlgebra @test klu!(G, nonzeros(C)) \ ones(Tv, 2) ≈ Matrix(C) \ ones(Tv, 2) @test !issuccess(klu!(G, singular; check = false)) @test issuccess(klu!(G, singular; allowsingular = true); allowsingular = true) + @test !issuccess(klu_refactor!(klu(C), singular; check = false)) + + # A refactorization that stops partway must not leave a mix of old and new factors. + G = klu(C) + @test_throws SingularException klu!(G, singular) + @test G._numeric == C_NULL + @test !issuccess(G) + @test_throws SingularException G \ ones(Tv, 2) + @test_throws SingularException rcond(G) + @test_throws SingularException condest(G) + @test !issuccess(G) + @test klu!(G, nonzeros(C)) \ ones(Tv, 2) ≈ Matrix(C) \ ones(Tv, 2) + @test issuccess(G) + + # A failed klu_factor! must not leave the factors of the earlier values in use. + G = klu(copy(C)) + nonzeros(G) .= one(Tv) + @test_throws SingularException klu_factor!(G) + @test !issuccess(G) + @test_throws SingularException G \ ones(Tv, 2) end @testset "Singular matrix" begin @@ -320,8 +341,10 @@ end @testset "Later solvers do not switch the allocator" begin script = """ using KLU, SparseArrays, LinearAlgebra + S = sparse([4.0 1.0; 1.0 3.0]) + F = klu(S) before = KLU.SuiteSparse_config_malloc_func_get() - lu(sparse([4.0 1.0; 1.0 3.0])) + lu(S) KLU.SuiteSparse_config_malloc_func_get() == before || exit(1) """ cmd = `$(Base.julia_cmd()) --startup-file=no --project=$(Base.active_project()) -e $script` @@ -437,6 +460,30 @@ end end end +@testset "Factorizations sharing pattern vectors" begin + n = 20_000 + A = sparse(Tridiagonal(fill(-1.0, n - 1), fill(4.0, n), fill(-1.0, n - 1))) + colptr, rowval = A.colptr .- 1, A.rowval .- 1 + K1 = klu(n, colptr, rowval, copy(nonzeros(A))) + K2 = klu(n, colptr, rowval, copy(nonzeros(A))) + @test K1.colptr === K2.colptr + # libklu would read out of bounds on K2 if refactoring K1 shifted the shared vectors. + stop = Threads.Atomic{Bool}(false) + shifted = Threads.Atomic{Bool}(false) + watcher = Threads.@spawn while !stop[] + colptr[1] == 0 || (shifted[] = true) + yield() + end + for _ in 1:50 + klu!(K1, A) + end + stop[] = true + wait(watcher) + @test !shifted[] + b = ones(n) + @test A * (K2 \ b) ≈ b +end + @testset "klu_factor! on an already factored object releases the old factor" begin A = sparse([1, 2], [1, 2], [1.0, 2.0]) F = klu(A) @@ -491,5 +538,249 @@ end rgrowth(F) # libklu resets the status here @test F.common.status == KLU.KLU_SINGULAR @test !issuccess(F) + # All of these keep a status they did not set themselves, as libklu resets it. + F = klu(sparse(Tv[1 2; 0 1])) + F.common.status = KLU.KLU_SINGULAR + for op in (rcond, condest, rgrowth) + op(F) + @test F.common.status == KLU.KLU_SINGULAR + end + end +end + +@testset "A factorization shares the vectors it is given" begin + b = ones(3) + # An explicitly stored zero, so that dropzeros! shrinks the vectors of A. + A = sparse([1, 2, 3, 1], [1, 2, 3, 3], [1.0, 2.0, 3.0, 0.0]) + Ad = Matrix(A) + + K = klu(A; full_factor = false) + @test nonzeros(K) === nonzeros(A) + nonzeros(A) .*= 2 + @test klu_factor!(K) \ b ≈ 2 * Ad \ b + nonzeros(A) ./= 2 + @test nonzeros(KLUFactorization(A)) !== nonzeros(A) + @test nonzeros(klu(SparseMatrixCSC{Float32, Int}(A))) isa Vector{Float64} + + colptr, rowval, nzval = [0, 1, 2, 4], [0, 1, 0, 2], [1.0, 2.0, 0.0, 3.0] + K = klu(3, colptr, rowval, nzval) + @test K.colptr === colptr && K.rowval === rowval && nonzeros(K) === nzval + v = 2 * nzval + @test klu!(K, v) \ b ≈ 2 * Ad \ b + @test nonzeros(K) === v + # Changing the values doesn't change the current factors. + v .*= 2 + @test K \ b ≈ 2 * Ad \ b + + # libklu would read past the end of a vector that the caller made too short. + A2 = copy(A) + K = klu(A2) + x = K \ b + full = copy(nonzeros(K)) + dropzeros!(A2) + @test length(nonzeros(K)) < K.colptr[end] + for op in (klu_factor!, rgrowth, condest, copy) + @test_throws ArgumentError op(K) + end + @test_throws DimensionMismatch klu!(K, full) + @test K \ b ≈ x # the factors are unaffected + A2 = copy(A) + U = klu(A2; full_factor = false) + dropzeros!(A2) + @test_throws ArgumentError klu_factor!(U) + @test_throws ArgumentError U \ b + colptr, rowval = [0, 1, 2, 4], [0, 1, 0, 2] + U = KLUFactorization(3, colptr, rowval, [1.0, 2.0, 0.0, 3.0]) + resize!(rowval, 2) + @test_throws ArgumentError klu_analyze!(U) + @test_throws ArgumentError klu_analyze!(U, [0, 1, 2], [0, 1, 2]) +end + +@testset "Fields that can't be reassigned" begin + A = sparse([2.0 1.0; 0.0 2.0]) + K = klu(A) + for s in (:common, :n, :colptr, :rowval, :lock) + @test_throws ErrorException setproperty!(K, s, getfield(K, s)) + end + @test K \ [1.0, 2.0] ≈ Matrix(A) \ [1.0, 2.0] +end + +@testset "Empty matrices" begin + for f in (() -> klu(spzeros(0, 0)), () -> klu(spzeros(0, 0); check = false), + () -> KLUFactorization(spzeros(0, 0)), () -> klu(0, [0], Int[], Float64[])) + @test_throws "empty matrix" f() + end +end + +@testset "Determinant, $Tv elements, $Ti indices" for Tv in (Float64, ComplexF64), Ti in Base.uniontypes(KLUITypes) + # Several blocks, nontrivial permutations, and off-diagonal entries in F. + A0 = sparse([1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 1], [1, 3, 2, 5, 3, 1, 5, 4, 2, 5, 4, 2], + [2.0, 3.0, -1.0, 4.0, 5.0, 1.0, 7.0, -2.0, 1.0, 3.0, 6.0, 9.0]) + @test klu(A0).nblocks > 1 + for B in (A0, A0[[2, 1, 3, 4, 5], :], A0[:, [5, 4, 3, 2, 1]], A0[[3, 1, 2, 5, 4], [2, 1, 4, 3, 5]]) + A = convert(SparseMatrixCSC{Tv, Ti}, B) + Tv <: Complex && (nonzeros(A) .*= cis.(1:nnz(A))) + Ad = Matrix(A) + F = klu(A) + ld, s = logabsdet(F) + ld0, s0 = logabsdet(Ad) + @test ld ≈ ld0 + @test s ≈ s0 + @test abs(s) ≈ 1 + @test det(F) ≈ det(Ad) + @test det(F') ≈ det(Ad') + @test det(transpose(F)) ≈ det(transpose(Ad)) + @test logabsdet(F')[2] ≈ conj(s0) + @test logabsdet(transpose(F))[2] ≈ s0 + # With other scalings, which change Rs but not the determinant. + for scale in (0, 1, 2) + K = KLUFactorization(A) + K.common.scale = scale + @test det(klu_factor!(K)) ≈ det(Ad) + end + # A factorization which was not yet computed is computed on demand. + @test det(klu(A; full_factor = false)) ≈ det(Ad) + @test det(klu!(klu(A), nonzeros(2 * A))) ≈ 2^5 * det(Ad) + end + S = convert(SparseMatrixCSC{Tv, Ti}, sparse([1.0 2.0; 0.0 0.0])) + F = klu(S; allowsingular = true) + @test logabsdet(F) == (-Inf, 0) + @test det(F) == 0 + @test F.common.status == KLU.KLU_SINGULAR + @test_throws SingularException det(klu(S; check = false, full_factor = false)) + # A large determinant doesn't overflow in logabsdet. + L = convert(SparseMatrixCSC{Tv, Ti}, sparse(1:400, 1:400, fill(1e10, 400))) + @test logabsdet(klu(L))[1] ≈ 400 * log(1e10) +end + +@testset "The column reported by SingularException" begin + # A status alone doesn't say which column. + @test_throws SingularException(0) KLU.kluerror(KLU.KLU_SINGULAR) + for Tv in (Float64, ComplexF64), Ti in Base.uniontypes(KLUITypes) + @test_throws SingularException(2) klu(convert(SparseMatrixCSC{Tv, Ti}, sparse([1.0 2.0; 0.0 0.0]))) + @test_throws SingularException(3) klu(convert(SparseMatrixCSC{Tv, Ti}, sparse([1, 2], [1, 2], [1.0, 1.0], 3, 3))) + C = convert(SparseMatrixCSC{Tv, Ti}, sparse([1, 2, 1, 2], [1, 1, 2, 2], [1.0, 1.0, 1.0, 2.0])) + e = try klu!(klu(C), Tv[1, 1, 1, 1]) catch e; e end + @test e isa SingularException && e.info in 1:2 end end + +@testset "Failure paths and state changes" begin + A = sparse([2.0 1.0; 0.0 2.0]) + b = [1.0, 2.0] + x = Matrix(A) \ b + + # An earlier failure doesn't stop a retry, which a temporary lack of memory should allow. + K = klu(A) + K.common.status = KLU.KLU_OUT_OF_MEMORY + @test issuccess(klu_factor!(K)) + @test K \ b ≈ x + K = klu(A; full_factor = false) + K.common.status = KLU.KLU_OUT_OF_MEMORY + @test issuccess(klu_factor!(K)) + K = KLUFactorization(A) + K.common.status = KLU.KLU_OUT_OF_MEMORY + @test issuccess(klu_factor!(K)) + # What libklu rejects is still rejected. + K = klu(2, [0, 5, 1], [0], [1.0]; check = false) + @test_throws ArgumentError klu_factor!(K) + @test K.common.status == KLU.KLU_INVALID + + # check = false also applies to a factorization that a solve computes on demand. + for op in (identity, adjoint, transpose) + K = klu(2, [0, 5, 1], [0], [1.0]; check = false) + B = ones(2) + @test solve!(op(K), B; check = false) === B + @test K.common.status == KLU.KLU_INVALID + @test !issuccess(K) + @test_throws ArgumentError solve!(op(K), ones(2)) + S = sparse([1.0 2.0; 0.0 0.0]) + K = klu(S; full_factor = false) + @test solve!(op(K), ones(2); check = false) isa Vector{Float64} + @test K.common.status == KLU.KLU_SINGULAR + K = klu(S; full_factor = false) + @test_throws SingularException solve!(op(K), ones(2)) + end + + # Assigning the current pointer must not free it, which would leave it dangling. + K = klu(A) + K._numeric = K._numeric + K._symbolic = K._symbolic + @test K \ b ≈ x + # Numeric factors must not outlive the analysis they belong to. + K._symbolic = C_NULL + @test K._numeric == C_NULL + @test !issuccess(K) + @test K \ b ≈ x + + # The permutations are the ones that were validated, not the caller's vectors. + P, Q = [1, 0], [0, 1] + K = klu_analyze!(KLUFactorization(A), P, Q) + @test K._perms[1] !== P && K._perms[2] !== Q + P .= 0 + @test K._perms == ([1, 0], [0, 1]) + @test copy(K) \ b ≈ x + + # Column stride does not matter with a single row, as libklu never steps along a column. + F = klu(sparse([2.0;;])) + M = [1.0 2.0; 3.0 4.0] + @test ldiv!(F, view(M, 1, 1:1)) == [0.5] + @test M == [0.5 2.0; 3.0 4.0] + + # Reading state takes the lock, so these wait for another task holding it. + K = klu(A) + for f in (issuccess, nnz, K -> sprint(show, MIME("text/plain"), K)) + lock(K) + t = Threads.@spawn f(K) + sleep(0.05) + @test !istaskdone(t) + unlock(K) + wait(t) + @test !istaskfailed(t) + end + + # Printing has all the factors from one factorization, without holding the lock while writing. + K = klu(A) + @test sprint(show, MIME("text/plain"), K) == + "$(summary(K))\nL factor:\n$(sprint(show, MIME("text/plain"), K.L))\nU factor:\n$(sprint(show, MIME("text/plain"), K.U))\nF factor:\n$(sprint(show, MIME("text/plain"), K.F))" + @test !islocked(K) +end + +@testset "Finalizers release the libklu objects" begin + A = sparse([2.0 1.0; 0.0 2.0]) + b = [1.0, 2.0] + x = Matrix(A) \ b + for make in (() -> klu(A), () -> klu(A; full_factor = false), () -> klu_factor!(KLUFactorization(A))) + K = make() + @test K._symbolic != C_NULL + K \ b + @test K._numeric != C_NULL + finalize(K) + @test K._symbolic == C_NULL && K._numeric == C_NULL + end + # An object with nothing to release is fine, and can still be used afterwards. + K = KLUFactorization(A) + finalize(K) + @test K \ b ≈ x + # Analyzing again adds a finalizer call, which finds nothing left to free. + K = klu(A) + K._symbolic = C_NULL + @test K \ b ≈ x + finalize(K) + @test K._symbolic == C_NULL && K._numeric == C_NULL + + # `deserialize` does not call the constructor, and libklu's objects are not serialized. + K = klu(A) + io = IOBuffer() + serialize(io, K) + seekstart(io) + D = deserialize(io) + @test D isa KLUFactorization{Float64, Int} + @test D._symbolic == C_NULL && D._numeric == C_NULL + @test D \ b ≈ x + @test D._symbolic != C_NULL && D._numeric != C_NULL + @test D._numeric != K._numeric + finalize(D) + @test D._symbolic == C_NULL && D._numeric == C_NULL + @test K \ b ≈ x +end