diff --git a/src/linalg.jl b/src/linalg.jl index c76f26a6..7436c68c 100644 --- a/src/linalg.jl +++ b/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) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 97ca426e..3a359293 100644 --- a/src/solvers/cholmod.jl +++ b/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) diff --git a/src/solvers/umfpack.jl b/src/solvers/umfpack.jl index 26d6c8d6..a85efb34 100644 --- a/src/solvers/umfpack.jl +++ b/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 diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index ecadccfd..5dcf5402 100644 --- a/src/sparsematrix.jl +++ b/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 diff --git a/test/linalg.jl b/test/linalg.jl index 7b8aba54..e7533255 100644 --- a/test/linalg.jl +++ b/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 diff --git a/test/umfpack.jl b/test/umfpack.jl index b0d53e4e..92bc3986 100644 --- a/test/umfpack.jl +++ b/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