Backports to release 1.10

src/linalg.jl

@@ -98,6 +98,11 @@ end
     (T = promote_op(matprod, eltype(A), eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B))
 
 Base.@constprop :aggressive function LinearAlgebra.generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::AbstractSparseMatrixCSC, _add::MulAddMul)
+    if !(tA in ('N', 'T', 'C') && tB in ('N', 'T', 'C'))
+        # redirect to most generic matmatmul code
+        LinearAlgebra._generic_matmatmul!(C, 'N', 'N', LinearAlgebra.wrap(A, tA), LinearAlgebra.wrap(B, tB), _add)
+        return C
+    end
     transA = tA == 'N' ? identity : tA == 'T' ? transpose : adjoint
     if tB == 'N'
         _spmul!(C, transA(A), B, _add.alpha, _add.beta)
@@ -401,7 +406,7 @@ const WrapperMatrixTypes{T,MT} = Union{
     Hermitian{T,MT},
 }
 
-function dot(A::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,Union{DenseMatrixUnion,AbstractSparseMatrix}}}, B::AbstractSparseMatrixCSC)
+function dot(A::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,<:Union{DenseMatrixUnion,AbstractSparseMatrix}}}, B::AbstractSparseMatrixCSC)
     T = promote_type(eltype(A), eltype(B))
     (m, n) = size(A)
     if (m, n) != size(B)
@@ -423,7 +428,7 @@ function dot(A::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,Union{DenseMatri
     return s
 end
 
-function dot(A::AbstractSparseMatrixCSC, B::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,Union{DenseMatrixUnion,AbstractSparseMatrix}}})
+function dot(A::AbstractSparseMatrixCSC, B::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,<:Union{DenseMatrixUnion,AbstractSparseMatrix}}})
     return conj(dot(B, A))
 end
 
@@ -913,7 +918,7 @@ function nzrangelo(A, i, excl=false)
     @inbounds r2 < r1 || rv[r1] >= i + excl ? r : searchsortedfirst(rv, i + excl, r1, r2, Forward):r2
 end
 
-dot(x::AbstractVector, A::RealHermSymComplexHerm{<:Any,<:AbstractSparseMatrixCSC}, y::AbstractVector) =
+dot(x::AbstractVector, A::RealHermSymComplexHerm{<:Real,<:AbstractSparseMatrixCSC}, y::AbstractVector) =
     _dot(x, parent(A), y, A.uplo == 'U' ? nzrangeup : nzrangelo, A isa Symmetric ? identity : real, A isa Symmetric ? transpose : adjoint)
 function _dot(x::AbstractVector, A::AbstractSparseMatrixCSC, y::AbstractVector, rangefun::Function, diagop::Function, odiagop::Function)
     require_one_based_indexing(x, y)
@@ -944,7 +949,7 @@ function _dot(x::AbstractVector, A::AbstractSparseMatrixCSC, y::AbstractVector,
     end
     return r
 end
-dot(x::SparseVector, A::RealHermSymComplexHerm{<:Any,<:AbstractSparseMatrixCSC}, y::SparseVector) =
+dot(x::SparseVector, A::RealHermSymComplexHerm{<:Real,<:AbstractSparseMatrixCSC}, y::SparseVector) =
     _dot(x, parent(A), y, A.uplo == 'U' ? nzrangeup : nzrangelo, A isa Symmetric ? identity : real)
 function _dot(x::SparseVector, A::AbstractSparseMatrixCSC, y::SparseVector, rangefun::Function, diagop::Function)
     m, n = size(A)

src/solvers/cholmod.jl

@@ -36,7 +36,7 @@ import SparseArrays: AbstractSparseMatrix, SparseMatrixCSC, indtype, sparse, spz
 import ..increment, ..increment!
 
 using ..LibSuiteSparse
-import ..LibSuiteSparse: TRUE, FALSE, CHOLMOD_INT, CHOLMOD_INTLONG, 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
@@ -221,16 +221,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)

src/solvers/umfpack.jl

@@ -261,9 +261,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)),
@@ -295,29 +297,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)
 
 
 """
@@ -584,7 +569,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)
@@ -599,7 +584,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,
@@ -613,7 +598,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)
@@ -629,7 +614,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)
@@ -658,7 +643,7 @@ for itype in UmfpackIndexTypes
             if workspace_W_size(lu) > length(workspace.W)
                 throw(ArguementError("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())
 
@@ -683,7 +668,7 @@ for itype in UmfpackIndexTypes
             if workspace_W_size(lu) > length(workspace.W)
                 throw(ArguementError("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,
@@ -693,14 +678,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)
@@ -1016,7 +1001,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
@@ -1026,7 +1011,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

src/sparsematrix.jl

@@ -3984,9 +3984,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)
@@ -4022,6 +4022,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

test/linalg.jl

@@ -224,6 +224,13 @@ end
     end
 end
 
+@testset "Dense times symmetric/Hermitian sparse matrix multiplication" begin
+    A = [1 3; 2 4]
+    As = sparse(A)
+    B = [1 1; 1 1]
+    @test mul!(copy(B), B, Hermitian(A), true, true) == mul!(copy(B), B, Hermitian(As), true, true)
+end
+
 @testset "UniformScaling" begin
     local A = sprandn(10, 10, 0.5)
     @test A + I == Array(A) + I
@@ -395,6 +402,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
@@ -824,11 +852,12 @@ end
         @test dot(x, A, y) ≈ dot(x, Av, y)
     end
 
-    for (T, trans) in ((Float64, Symmetric), (ComplexF64, Hermitian)), uplo in (:U, :L)
+    for (T, trans) in ((Float64, Symmetric), (ComplexF64, Symmetric), (ComplexF64, Hermitian)), uplo in (:U, :L)
         B = sprandn(T, 10, 10, 0.2)
         x = sprandn(T, 10, 0.4)
         S = trans(B'B, uplo)
-        @test dot(x, S, x) ≈ dot(Vector(x), S, Vector(x)) ≈ dot(Vector(x), Matrix(S), Vector(x))
+        Sd = trans(Matrix(B'B), uplo)
+        @test dot(x, S, x) ≈ dot(x, Sd, x) ≈ dot(Vector(x), S, Vector(x)) ≈ dot(Vector(x), Sd, Vector(x))
     end
 end
 

test/umfpack.jl

@@ -121,8 +121,6 @@ end
         umfpack_report(Af)
         Af1 = lu!(copy(Af))
         umfpack_report(Af1)
-        @test trylock(Af)
-        @test trylock(Af1)
     end
 
     @testset "test similar" begin