From 7ec2889dc855a5d6c201552052a2ff9a7b3e1e16 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 1 Apr 2025 07:49:47 +0530 Subject: [PATCH 01/10] Fix `issymmetric` for matrices with empty columns (#606) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The `issymmetric` check tracks an `offset` that it uses to go from (row, col) to (col, row). However, currently this doesn't account for the fact that if a column is empty, entries in `colptr` will be identical. E.g., in #605, we have ```julia julia> S = sparse([2, 3, 1], [1, 1, 3], [1, 1, 1], 3, 3) 3×3 SparseMatrixCSC{Int64, Int64} with 3 stored entries: ⋅ ⋅ 1 1 ⋅ ⋅ 1 ⋅ ⋅ julia> SparseArrays.getcolptr(S) 4-element Vector{Int64}: 1 3 3 4 ``` The offset `3` corresponds to rows in the third column, as the second column is empty. This PR checks for empty columns, in which case we may exit the call immediately. Fixes #605 --- src/sparsematrix.jl | 10 ++++++++-- test/linalg.jl | 21 +++++++++++++++++++++ 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index c42c7777..7225944d 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -4074,9 +4074,9 @@ function _blockdiag(::Type{Tv}, ::Type{Ti}, X::AbstractSparseMatrixCSC...) where end ## Structure query functions -issymmetric(A::AbstractSparseMatrixCSC) = is_hermsym(A, identity) +issymmetric(A::AbstractSparseMatrixCSC) = is_hermsym(A, transpose) -ishermitian(A::AbstractSparseMatrixCSC) = is_hermsym(A, conj) +ishermitian(A::AbstractSparseMatrixCSC) = is_hermsym(A, adjoint) function is_hermsym(A::AbstractSparseMatrixCSC, check::Function) m, n = size(A) @@ -4112,6 +4112,12 @@ function is_hermsym(A::AbstractSparseMatrixCSC, check::Function) return false end else + # if nzrange(A, row) is empty, then A[:, row] is all zeros. + # Specifically, A[col, row] is zero. + # However, we know at this point that A[row, col] is not zero + # This means that the matrix is not symmetric + isempty(nzrange(A, row)) && return false + offset = tracker[row] # If the matrix is unsymmetric, there might not exist diff --git a/test/linalg.jl b/test/linalg.jl index b6113c50..afbcb2ad 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -425,6 +425,27 @@ end @test issymmetric(sparse([1 0; 1 0])) == false @test issymmetric(sparse([0 1; 1 0])) == true @test issymmetric(sparse([1 1; 1 0])) == true + + # test some non-trivial cases + local S + @testset "random matrices" begin + for sparsity in (0.1, 0.01, 0.0) + S = sparse(Symmetric(sprand(20, 20, sparsity))) + @test issymmetric(S) + @test ishermitian(S) + S = sparse(Symmetric(sprand(ComplexF64, 20, 20, sparsity))) + @test issymmetric(S) + @test !ishermitian(S) || isreal(S) + S = sparse(Hermitian(sprand(ComplexF64, 20, 20, sparsity))) + @test ishermitian(S) + @test !issymmetric(S) || isreal(S) + end + end + + @testset "issue #605" begin + S = sparse([2, 3, 1], [1, 1, 3], [1, 1, 1], 3, 3) + @test !issymmetric(S) + end end @testset "rotations" begin From 08a15ca4eade031b619ef0b09db7440bcbbfcc2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20H=C3=A9not?= <38465572+OlivierHnt@users.noreply.github.com> Date: Wed, 26 Mar 2025 13:40:50 +0100 Subject: [PATCH 02/10] Replace `v == zero(v)` with `_iszero` (#610) The purpose of this PR is to not call `==` directly, and use `_iszero` instead. This is to help integration with [IntervalArithmetic.jl](https://github.com/JuliaIntervals/IntervalArithmetic.jl) since `==` for our `Interval` type does not always return a Boolean. Closes https://github.com/JuliaSparse/SparseArrays.jl/issues/609. I did not change two checks using `===` since this would break the behaviour for `-0.0`. --- src/sparsematrix.jl | 6 +++--- src/sparsevector.jl | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 7225944d..38a6748b 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -2283,7 +2283,7 @@ function (+)(A::Array, B::SparseMatrixCSCUnion) for j in axes(B,2) for i in nzrange(B, j) rowidx = rowinds[i] - C[rowidx,j] = A[rowidx,j] + nzvals[i] + C[rowidx,j] = A[rowidx,j] + nzvals[i] end end return C @@ -2452,7 +2452,7 @@ function Base._mapreduce(f, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.I else v = f(zero(T))^(nzeros) # Bail out early if initial reduction value is zero or if there are no stored elements - (v == zero(T) || nnzA == 0) ? v : v*Base._mapreduce(f, op, nzvalview(A)) + (_iszero(v) || nnzA == 0) ? v : v*Base._mapreduce(f, op, nzvalview(A)) end end @@ -4611,6 +4611,6 @@ function copytrito!(M::AbstractMatrix, S::AbstractSparseMatrixCSC, uplo::Char) (uplo == 'U' && row <= col) || (uplo == 'L' && row >= col) || continue M[row, col] = nz[i] end - end + end return M end diff --git a/src/sparsevector.jl b/src/sparsevector.jl index e5efe4ac..6734b986 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1932,7 +1932,7 @@ function _spmul!(y::AbstractVector, A::AbstractMatrix, x::AbstractSparseVector, "Matrix A has $m rows, but vector y has a length $(length(y))")) m == 0 && return y β != one(β) && LinearAlgebra._rmul_or_fill!(y, β) - α == zero(α) && return y + _iszero(α) && return y xnzind = nonzeroinds(x) xnzval = nonzeros(x) @@ -1960,7 +1960,7 @@ function _At_or_Ac_mul_B!(tfun::Function, "Matrix A has $m columns, but vector y has a length $(length(y))")) m == 0 && return y β != one(β) && LinearAlgebra._rmul_or_fill!(y, β) - α == zero(α) && return y + _iszero(α) && return y xnzind = nonzeroinds(x) xnzval = nonzeros(x) @@ -2055,7 +2055,7 @@ function _spmul!(y::AbstractVector, A::AbstractSparseMatrixCSC, x::AbstractSpars "Matrix A has $m rows, but vector y has a length $(length(y))")) m == 0 && return y β != one(β) && LinearAlgebra._rmul_or_fill!(y, β) - α == zero(α) && return y + _iszero(α) && return y xnzind = nonzeroinds(x) xnzval = nonzeros(x) @@ -2087,7 +2087,7 @@ function _At_or_Ac_mul_B!(tfun::Function, "Matrix A has $m rows, but vector y has a length $(length(y))")) n == 0 && return y β != one(β) && LinearAlgebra._rmul_or_fill!(y, β) - α == zero(α) && return y + _iszero(α) && return y xnzind = nonzeroinds(x) xnzval = nonzeros(x) @@ -2421,7 +2421,7 @@ import Base.fill! function fill!(A::Union{AbstractCompressedVector, AbstractSparseMatrixCSC}, x) T = eltype(A) xT = convert(T, x) - if xT == zero(T) + if _iszero(xT) fill!(nonzeros(A), xT) else _fillnonzero!(A, xT) From 2ae8768d07837b956882c51bfb20f5ebe22d3464 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 7 Apr 2025 08:28:40 +0530 Subject: [PATCH 03/10] `ldiv!` for `Diagonal` and a sparse vector (#613) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This improves performance: ```julia julia> using LinearAlgebra, SparseArrays julia> D = Diagonal(rand(3000)); julia> S = sprand(size(D,1), 0.01); julia> @btime ldiv!($D, $S); 30.053 μs (0 allocations: 0 bytes) # master 1.585 μs (0 allocations: 0 bytes) # PR ``` --- src/linalg.jl | 2 +- test/linalg.jl | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index a7dbf350..33b1f1f2 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1312,7 +1312,7 @@ function rdiv!(A::AbstractSparseMatrixCSC{T}, D::Diagonal{T}) where T A end -function ldiv!(D::Diagonal{T}, A::AbstractSparseMatrixCSC{T}) where {T} +function ldiv!(D::Diagonal{T}, A::Union{AbstractSparseMatrixCSC{T}, AbstractSparseVector{T}}) where {T} # require_one_based_indexing(A) if size(A, 1) != length(D.diag) throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(size(A, 1)) rows")) diff --git a/test/linalg.jl b/test/linalg.jl index afbcb2ad..d32adcc4 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -338,6 +338,9 @@ end @test lmul!(transpose(copy(D)), copy(b)) ≈ transpose(MD)*bd @test lmul!(adjoint(copy(D)), copy(b)) ≈ MD'*bd end + + v = sprand(eltype(D), size(D,1), 0.1) + @test ldiv!(D, copy(v)) == D \ Array(v) end end From 66d65cc64c160984034a9ef68166e9205c25b00b Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 21 Apr 2025 23:50:30 +0530 Subject: [PATCH 04/10] Relax `eltype` in `Diagonal` `ldiv!`/`rdiv!` (#616) These methods do not require the `eltype`s to match exactly, and should work as long as the values may be stored in the destination. --- src/linalg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 33b1f1f2..04095dfc 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1294,7 +1294,7 @@ end \(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:AbstractSparseMatrixCSC}}, B::Vector) = copy(A) \ B -function rdiv!(A::AbstractSparseMatrixCSC{T}, D::Diagonal{T}) where T +function rdiv!(A::AbstractSparseMatrixCSC, D::Diagonal) dd = D.diag if (k = length(dd)) ≠ size(A, 2) throw(DimensionMismatch("size(A, 2)=$(size(A, 2)) should be size(D, 1)=$k")) @@ -1312,7 +1312,7 @@ function rdiv!(A::AbstractSparseMatrixCSC{T}, D::Diagonal{T}) where T A end -function ldiv!(D::Diagonal{T}, A::Union{AbstractSparseMatrixCSC{T}, AbstractSparseVector{T}}) where {T} +function ldiv!(D::Diagonal, A::Union{AbstractSparseMatrixCSC, AbstractSparseVector}) # require_one_based_indexing(A) if size(A, 1) != length(D.diag) throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(size(A, 1)) rows")) From b38f273a7dfea9a9e87e7598c59d99831f9fc0ab Mon Sep 17 00:00:00 2001 From: Elliot Saba Date: Fri, 9 May 2025 04:50:46 -0700 Subject: [PATCH 05/10] Use `libsuitesparseconfig` from JLL (#620) This prevents us from being influenced by things like `LD_LIBRARY_PATH`, and ensures that we always load the correct `libsuitesparseconfig` that came with our JLL. --- src/solvers/cholmod.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 7cca3775..116ca2e5 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -231,16 +231,16 @@ function __init__() # Register gc tracked allocator if CHOLMOD is new enough if current_version >= v"4.0.3" - ccall((:SuiteSparse_config_malloc_func_set, :libsuitesparseconfig), + ccall((:SuiteSparse_config_malloc_func_set, libsuitesparseconfig), Cvoid, (Ptr{Cvoid},), cglobal(:jl_malloc, Ptr{Cvoid})) - ccall((:SuiteSparse_config_calloc_func_set, :libsuitesparseconfig), + ccall((:SuiteSparse_config_calloc_func_set, libsuitesparseconfig), Cvoid, (Ptr{Cvoid},), cglobal(:jl_calloc, Ptr{Cvoid})) - ccall((:SuiteSparse_config_realloc_func_set, :libsuitesparseconfig), + ccall((:SuiteSparse_config_realloc_func_set, libsuitesparseconfig), Cvoid, (Ptr{Cvoid},), cglobal(:jl_realloc, Ptr{Cvoid})) - ccall((:SuiteSparse_config_free_func_set, :libsuitesparseconfig), + ccall((:SuiteSparse_config_free_func_set, libsuitesparseconfig), Cvoid, (Ptr{Cvoid},), cglobal(:jl_free, Ptr{Cvoid})) elseif current_version >= v"3.0.0" - cnfg = cglobal((:SuiteSparse_config, :libsuitesparseconfig), Ptr{Cvoid}) + cnfg = cglobal((:SuiteSparse_config, libsuitesparseconfig), Ptr{Cvoid}) unsafe_store!(cnfg, cglobal(:jl_malloc, Ptr{Cvoid}), 1) unsafe_store!(cnfg, cglobal(:jl_calloc, Ptr{Cvoid}), 2) unsafe_store!(cnfg, cglobal(:jl_realloc, Ptr{Cvoid}), 3) From e1817e8f4f31120abde39d8006ae84bbb0c0eafe Mon Sep 17 00:00:00 2001 From: Daniel Wennberg Date: Thu, 8 May 2025 05:59:27 -0700 Subject: [PATCH 06/10] Clarify pros, cons and limitations of Cholesky and LDLt (#621) I find it worth pointing out explicitly in the docs that LDLt, which mathematically looks like a drop-in replacement for Cholesky that does away with the positive definiteness requirement, comes with the following caveats: * It fails for a lot of matrices (for example, `ldlt(Symmetric(sprandn(1000, 1000, p)))` basically never succeeds for any relevant sparsity `p`) due to the requirement that all leading principal minors be well-conditioned * In CHOLMOD, `ldlt` is significantly slower than `cholesky` as it does not have a supernodal implementation So I made some docstring edits to clarify the relationship and tradeoffs between `cholesky` and `ldlt`. Citation for these claims: pages 106-107 in the CHOLMOD user guide at https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/v7.10.3/CHOLMOD/Doc/CHOLMOD_UserGuide.pdf --- src/solvers/cholmod.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 116ca2e5..557412cb 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -1563,6 +1563,9 @@ Setting the optional `shift` keyword argument computes the factorization of it should be a permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's default AMD ordering). +See also [`ldlt`](@ref) for a similar factorization that does not require +positive definiteness, but can be significantly slower than `cholesky`. + # Examples In the following example, the fill-reducing permutation used is `[3, 2, 1]`. @@ -1728,6 +1731,10 @@ To include the effects of permutation, it is typically preferable to extract `P'*L`) and `LtP = F.UP` (the equivalent of `L'*P`). The complete list of supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`. +Unlike the related Cholesky factorization, the ``LDL'`` factorization does not +require `A` to be positive definite. However, it still requires all leading +principal minors to be well-conditioned and will fail if this is not satisfied. + 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. @@ -1737,6 +1744,9 @@ Setting the optional `shift` keyword argument computes the factorization of it should be a permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's default AMD ordering). +See also [`cholesky`](@ref) for a factorization that can be significantly +faster than `ldlt`, but requires `A` to be positive definite. + !!! note This method uses the CHOLMOD[^ACM887][^DavisHager2009] library from [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse). CHOLMOD only supports real or complex types in single or double precision. From fa9736cb31229d5104b73147a1ed78fed4b62a15 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 11 Jun 2025 13:00:43 +0200 Subject: [PATCH 07/10] [release-1.12] Run CI with Julia 1.12 --- .github/workflows/ci.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 77689df4..2ce5d4ec 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -5,7 +5,7 @@ on: branches: - 'main' - 'release-*' - tags: '*' + tags: ['*'] concurrency: # Skip intermediate builds: all builds except for builds on the `main` or `release-*` branches # Cancel intermediate builds: only pull request builds @@ -19,7 +19,7 @@ jobs: fail-fast: false matrix: julia-version: - - 'nightly' + - '~1.12.0-0' os: - ubuntu-latest - windows-latest @@ -29,10 +29,10 @@ jobs: include: - os: macOS-latest julia-arch: aarch64 - julia-version: 'nightly' + julia-version: '~1.12.0-0' - os: macOS-13 julia-arch: x64 - julia-version: 'nightly' + julia-version: '~1.12.0-0' steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 @@ -56,7 +56,7 @@ jobs: strategy: matrix: julia-version: - - 'nightly' + - '~1.12.0-0' os: - ubuntu-latest julia-arch: @@ -78,7 +78,7 @@ jobs: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 with: - version: 'nightly' + version: '~1.12.0-0' - name: Generate docs run: | julia --project --color=yes -e 'using Pkg; Pkg.activate("docs"); Pkg.develop(PackageSpec(path = pwd()))' From f51dace2f566a8748bb4dd19042e104df0e8274f Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 11 Jun 2025 13:11:01 +0200 Subject: [PATCH 08/10] Fix compat for SuiteSparse_jll --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 21cba784..f53533f5 100644 --- a/Project.toml +++ b/Project.toml @@ -19,7 +19,7 @@ Pkg = "<0.0.1, 1" Printf = "<0.0.1, 1" Random = "<0.0.1, 1" Serialization = "<0.0.1, 1" -SuiteSparse_jll = "7.10.1" +SuiteSparse_jll = "7.8.3" Test = "<0.0.1, 1" julia = "1.11" From d9213082a55334b3cd737aae3f913805319c4ac0 Mon Sep 17 00:00:00 2001 From: Ian Butterworth Date: Fri, 16 May 2025 23:33:23 -0400 Subject: [PATCH 09/10] fix libsuitesparseconfig bug (cherry picked from commit 75eaa18eaf8f568975cba48e0cc3cbd5b54dab1d, PR #625) --- src/solvers/cholmod.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 557412cb..dfcbd5bb 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -37,7 +37,7 @@ import SparseArrays: AbstractSparseMatrix, SparseMatrixCSC, indtype, sparse, spz import ..increment, ..increment! using ..LibSuiteSparse -import ..LibSuiteSparse: TRUE, FALSE, CHOLMOD_INT, CHOLMOD_LONG +import ..LibSuiteSparse: TRUE, FALSE, CHOLMOD_INT, CHOLMOD_LONG, libsuitesparseconfig # # itype defines the types of integer used: # CHOLMOD_INT, # all integer arrays are int From cdbad55530fba0c7aa27d4bcc64dde2204ff133f Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 11 Jun 2025 14:15:49 +0200 Subject: [PATCH 10/10] Delete lock interface implementation from UmfpackLU (#617) (cherry picked from commit 70921fbf6b3f67c957260f0e16b8657a87417e74, PR #617) --- src/solvers/umfpack.jl | 47 ++++++++++++++---------------------------- test/umfpack.jl | 2 -- 2 files changed, 16 insertions(+), 33 deletions(-) diff --git a/src/solvers/umfpack.jl b/src/solvers/umfpack.jl index 6e35ef95..b267f0b6 100644 --- a/src/solvers/umfpack.jl +++ b/src/solvers/umfpack.jl @@ -263,9 +263,11 @@ has_refinement(F::UmfpackLU) = has_refinement(F.control) has_refinement(control::AbstractVector) = control[JL_UMFPACK_IRSTEP] > 0 # auto magick resize, should this only expand and not shrink? -getworkspace(F::UmfpackLU) = @lock F.lock begin - return resize!(F.workspace, F, has_refinement(F); expand_only=true) +function getworkspace(F::UmfpackLU) + @lock F.lock begin + return resize!(F.workspace, F, has_refinement(F); expand_only = true) end +end UmfpackWS(F::UmfpackLU{Tv, Ti}, refinement::Bool=has_refinement(F)) where {Tv, Ti} = UmfpackWS( Vector{Ti}(undef, size(F, 2)), @@ -297,29 +299,12 @@ Base.copy(F::T, ws=UmfpackWS(F)) where {T <: ATLU} = Base.transpose(F::UmfpackLU) = TransposeFactorization(F) -function Base.lock(f::Function, F::UmfpackLU) - lock(F) - try - f() - finally - unlock(F) - end -end -Base.lock(F::UmfpackLU) = if !trylock(F.lock) - @info """waiting for UmfpackLU's lock, it's safe to ignore this message. - see the documentation for Umfpack""" maxlog = 1 - lock(F.lock) -end - -@inline Base.trylock(F::UmfpackLU) = trylock(F.lock) -@inline Base.unlock(F::UmfpackLU) = unlock(F.lock) - show_umf_ctrl(F::UmfpackLU, level::Real=2.0) = - @lock F show_umf_ctrl(F.control, level) + @lock F.lock show_umf_ctrl(F.control, level) show_umf_info(F::UmfpackLU, level::Real=2.0) = - @lock F show_umf_info(F.control, F.info, level) + @lock F.lock show_umf_info(F.control, F.info, level) """ @@ -598,7 +583,7 @@ for itype in UmfpackIndexTypes @eval begin function umfpack_symbolic!(U::UmfpackLU{Float64,$itype}, q::Union{Nothing, StridedVector{$itype}}) _isnotnull(U.symbolic) && return U - @lock U begin + @lock U.lock begin tmp = Ref{Ptr{Cvoid}}(C_NULL) if q === nothing @isok $sym_r(U.m, U.n, U.colptr, U.rowval, U.nzval, tmp, U.control, U.info) @@ -613,7 +598,7 @@ for itype in UmfpackIndexTypes end function umfpack_symbolic!(U::UmfpackLU{ComplexF64,$itype}, q::Union{Nothing, StridedVector{$itype}}) _isnotnull(U.symbolic) && return U - @lock U begin + @lock U.lock begin tmp = Ref{Ptr{Cvoid}}(C_NULL) if q === nothing @isok $sym_c(U.m, U.n, U.colptr, U.rowval, real(U.nzval), imag(U.nzval), tmp, @@ -627,7 +612,7 @@ for itype in UmfpackIndexTypes return U end function umfpack_numeric!(U::UmfpackLU{Float64,$itype}; reuse_numeric=true, q=nothing) - @lock U begin + @lock U.lock begin (reuse_numeric && _isnotnull(U.numeric)) && return U if _isnull(U.symbolic) umfpack_symbolic!(U, q) @@ -643,7 +628,7 @@ for itype in UmfpackIndexTypes return U end function umfpack_numeric!(U::UmfpackLU{ComplexF64,$itype}; reuse_numeric=true, q=nothing) - @lock U begin + @lock U.lock begin (reuse_numeric && _isnotnull(U.numeric)) && return U _isnull(U.symbolic) && umfpack_symbolic!(U, q) tmp = Ref{Ptr{Cvoid}}(C_NULL) @@ -672,7 +657,7 @@ for itype in UmfpackIndexTypes if workspace_W_size(lu) > length(workspace.W) throw(ArgumentError("W should be larger than `workspace_W_size(Af)`")) end - @lock lu begin + @lock lu.lock begin umfpack_numeric!(lu) (size(b, 1) == lu.m) && (size(b) == size(x)) || throw(DimensionMismatch()) @@ -697,7 +682,7 @@ for itype in UmfpackIndexTypes if workspace_W_size(lu) > length(workspace.W) throw(ArgumentError("W should be larger than `workspace_W_size(Af)`")) end - @lock lu begin + @lock lu.lock begin umfpack_numeric!(lu) (size(b, 1) == lu.m) && (size(b) == size(x)) || throw(DimensionMismatch()) @isok $wsol_c(typ, lu.colptr, lu.rowval, lu.nzval, C_NULL, x, C_NULL, b, @@ -707,14 +692,14 @@ for itype in UmfpackIndexTypes end function det(lu::UmfpackLU{Float64,$itype}) mx = Ref{Float64}(zero(Float64)) - @lock lu @isok($det_r(mx, C_NULL, lu.numeric, lu.info)) + @lock lu.lock @isok($det_r(mx, C_NULL, lu.numeric, lu.info)) mx[] end function det(lu::UmfpackLU{ComplexF64,$itype}) mx = Ref{Float64}(zero(Float64)) mz = Ref{Float64}(zero(Float64)) - @lock lu @isok($det_z(mx, mz, C_NULL, lu.numeric, lu.info)) + @lock lu.lock @isok($det_z(mx, mz, C_NULL, lu.numeric, lu.info)) complex(mx[], mz[]) end function logabsdet(F::UmfpackLU{T, $itype}) where {T<:Union{Float64,ComplexF64}} # return log(abs(det)) and sign(det) @@ -1030,7 +1015,7 @@ for Tv in (:Float64, :ComplexF64), Ti in UmfpackIndexTypes _report_symbolic = Symbol(umf_nm("report_symbolic", Tv, Ti)) @eval umfpack_report_symbolic(lu::UmfpackLU{$Tv,$Ti}, level::Real=4; q=nothing) = - @lock lu begin + @lock lu.lock begin umfpack_symbolic!(lu, q) old_prl = lu.control[JL_UMFPACK_PRL] lu.control[JL_UMFPACK_PRL] = level @@ -1040,7 +1025,7 @@ for Tv in (:Float64, :ComplexF64), Ti in UmfpackIndexTypes end _report_numeric = Symbol(umf_nm("report_numeric", Tv, Ti)) @eval umfpack_report_numeric(lu::UmfpackLU{$Tv,$Ti}, level::Real=4; q=nothing) = - @lock lu begin + @lock lu.lock begin umfpack_numeric!(lu; q) old_prl = lu.control[JL_UMFPACK_PRL] lu.control[JL_UMFPACK_PRL] = level diff --git a/test/umfpack.jl b/test/umfpack.jl index a00c4aca..59dc9041 100644 --- a/test/umfpack.jl +++ b/test/umfpack.jl @@ -124,8 +124,6 @@ end umfpack_report(Af) Af1 = lu!(copy(Af)) umfpack_report(Af1) - @test trylock(Af) - @test trylock(Af1) end @testset "test similar" begin