Backports for v1.12 (#624)

Author: dkarrasch

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"))

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)
@@ -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.

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

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)

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