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()))' 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" diff --git a/src/linalg.jl b/src/linalg.jl index a7dbf350..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::AbstractSparseMatrixCSC{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")) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 7cca3775..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 @@ -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) @@ -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. 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/src/sparsematrix.jl b/src/sparsematrix.jl index c42c7777..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 @@ -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 @@ -4605,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) diff --git a/test/linalg.jl b/test/linalg.jl index b6113c50..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 @@ -425,6 +428,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 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