Skip to content

Backports for v1.12 #624

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
May 11, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
Expand All @@ -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"))
Expand Down
20 changes: 15 additions & 5 deletions src/solvers/cholmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -231,16 +231,16 @@

# 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),

Check warning on line 236 in src/solvers/cholmod.jl

View check run for this annotation

Codecov / codecov/patch

src/solvers/cholmod.jl#L236

Added line #L236 was not covered by tests
Cvoid, (Ptr{Cvoid},), cglobal(:jl_calloc, Ptr{Cvoid}))
ccall((:SuiteSparse_config_realloc_func_set, :libsuitesparseconfig),
ccall((:SuiteSparse_config_realloc_func_set, libsuitesparseconfig),

Check warning on line 238 in src/solvers/cholmod.jl

View check run for this annotation

Codecov / codecov/patch

src/solvers/cholmod.jl#L238

Added line #L238 was not covered by tests
Cvoid, (Ptr{Cvoid},), cglobal(:jl_realloc, Ptr{Cvoid}))
ccall((:SuiteSparse_config_free_func_set, :libsuitesparseconfig),
ccall((:SuiteSparse_config_free_func_set, libsuitesparseconfig),

Check warning on line 240 in src/solvers/cholmod.jl

View check run for this annotation

Codecov / codecov/patch

src/solvers/cholmod.jl#L240

Added line #L240 was not covered by tests
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})

Check warning on line 243 in src/solvers/cholmod.jl

View check run for this annotation

Codecov / codecov/patch

src/solvers/cholmod.jl#L243

Added line #L243 was not covered by tests
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)
Expand Down Expand Up @@ -1563,6 +1563,9 @@
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]`.
Expand Down Expand Up @@ -1728,6 +1731,10 @@
`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.
Expand All @@ -1737,6 +1744,9 @@
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.
Expand Down
16 changes: 11 additions & 5 deletions src/sparsematrix.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
10 changes: 5 additions & 5 deletions src/sparsevector.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
24 changes: 24 additions & 0 deletions test/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
Loading