Skip to content

Commit c3e31aa

Browse files
dkarraschdkarrasch
authored
Backports to release 1.10 (#704)
- [x] #606 - [x] #617 - [x] #620 - [x] #625 - [x] #641 - [x] #642
2 parents 1882d61 + bc4fc20 commit c3e31aa

6 files changed

Lines changed: 70 additions & 47 deletions

File tree

‎src/linalg.jl‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -98,6 +98,11 @@ end
9898
(T = promote_op(matprod, eltype(A), eltype(B)); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B))
9999

100100
Base.@constprop :aggressive function LinearAlgebra.generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::AbstractSparseMatrixCSC, _add::MulAddMul)
101+
if !(tA in ('N', 'T', 'C') && tB in ('N', 'T', 'C'))
102+
# redirect to most generic matmatmul code
103+
LinearAlgebra._generic_matmatmul!(C, 'N', 'N', LinearAlgebra.wrap(A, tA), LinearAlgebra.wrap(B, tB), _add)
104+
return C
105+
end
101106
transA = tA == 'N' ? identity : tA == 'T' ? transpose : adjoint
102107
if tB == 'N'
103108
_spmul!(C, transA(A), B, _add.alpha, _add.beta)
@@ -401,7 +406,7 @@ const WrapperMatrixTypes{T,MT} = Union{
401406
Hermitian{T,MT},
402407
}
403408

404-
function dot(A::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,Union{DenseMatrixUnion,AbstractSparseMatrix}}}, B::AbstractSparseMatrixCSC)
409+
function dot(A::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,<:Union{DenseMatrixUnion,AbstractSparseMatrix}}}, B::AbstractSparseMatrixCSC)
405410
T = promote_type(eltype(A), eltype(B))
406411
(m, n) = size(A)
407412
if (m, n) != size(B)
@@ -423,7 +428,7 @@ function dot(A::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,Union{DenseMatri
423428
return s
424429
end
425430

426-
function dot(A::AbstractSparseMatrixCSC, B::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,Union{DenseMatrixUnion,AbstractSparseMatrix}}})
431+
function dot(A::AbstractSparseMatrixCSC, B::Union{DenseMatrixUnion,WrapperMatrixTypes{<:Any,<:Union{DenseMatrixUnion,AbstractSparseMatrix}}})
427432
return conj(dot(B, A))
428433
end
429434

@@ -913,7 +918,7 @@ function nzrangelo(A, i, excl=false)
913918
@inbounds r2 < r1 || rv[r1] >= i + excl ? r : searchsortedfirst(rv, i + excl, r1, r2, Forward):r2
914919
end
915920

916-
dot(x::AbstractVector, A::RealHermSymComplexHerm{<:Any,<:AbstractSparseMatrixCSC}, y::AbstractVector) =
921+
dot(x::AbstractVector, A::RealHermSymComplexHerm{<:Real,<:AbstractSparseMatrixCSC}, y::AbstractVector) =
917922
_dot(x, parent(A), y, A.uplo == 'U' ? nzrangeup : nzrangelo, A isa Symmetric ? identity : real, A isa Symmetric ? transpose : adjoint)
918923
function _dot(x::AbstractVector, A::AbstractSparseMatrixCSC, y::AbstractVector, rangefun::Function, diagop::Function, odiagop::Function)
919924
require_one_based_indexing(x, y)
@@ -944,7 +949,7 @@ function _dot(x::AbstractVector, A::AbstractSparseMatrixCSC, y::AbstractVector,
944949
end
945950
return r
946951
end
947-
dot(x::SparseVector, A::RealHermSymComplexHerm{<:Any,<:AbstractSparseMatrixCSC}, y::SparseVector) =
952+
dot(x::SparseVector, A::RealHermSymComplexHerm{<:Real,<:AbstractSparseMatrixCSC}, y::SparseVector) =
948953
_dot(x, parent(A), y, A.uplo == 'U' ? nzrangeup : nzrangelo, A isa Symmetric ? identity : real)
949954
function _dot(x::SparseVector, A::AbstractSparseMatrixCSC, y::SparseVector, rangefun::Function, diagop::Function)
950955
m, n = size(A)

‎src/solvers/cholmod.jl‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ import SparseArrays: AbstractSparseMatrix, SparseMatrixCSC, indtype, sparse, spz
3636
import ..increment, ..increment!
3737

3838
using ..LibSuiteSparse
39-
import ..LibSuiteSparse: TRUE, FALSE, CHOLMOD_INT, CHOLMOD_INTLONG, CHOLMOD_LONG
39+
import ..LibSuiteSparse: TRUE, FALSE, CHOLMOD_INT, CHOLMOD_LONG, libsuitesparseconfig
4040

4141
# # itype defines the types of integer used:
4242
# CHOLMOD_INT, # all integer arrays are int
@@ -221,16 +221,16 @@ function __init__()
221221

222222
# Register gc tracked allocator if CHOLMOD is new enough
223223
if current_version >= v"4.0.3"
224-
ccall((:SuiteSparse_config_malloc_func_set, :libsuitesparseconfig),
224+
ccall((:SuiteSparse_config_malloc_func_set, libsuitesparseconfig),
225225
Cvoid, (Ptr{Cvoid},), cglobal(:jl_malloc, Ptr{Cvoid}))
226-
ccall((:SuiteSparse_config_calloc_func_set, :libsuitesparseconfig),
226+
ccall((:SuiteSparse_config_calloc_func_set, libsuitesparseconfig),
227227
Cvoid, (Ptr{Cvoid},), cglobal(:jl_calloc, Ptr{Cvoid}))
228-
ccall((:SuiteSparse_config_realloc_func_set, :libsuitesparseconfig),
228+
ccall((:SuiteSparse_config_realloc_func_set, libsuitesparseconfig),
229229
Cvoid, (Ptr{Cvoid},), cglobal(:jl_realloc, Ptr{Cvoid}))
230-
ccall((:SuiteSparse_config_free_func_set, :libsuitesparseconfig),
230+
ccall((:SuiteSparse_config_free_func_set, libsuitesparseconfig),
231231
Cvoid, (Ptr{Cvoid},), cglobal(:jl_free, Ptr{Cvoid}))
232232
elseif current_version >= v"3.0.0"
233-
cnfg = cglobal((:SuiteSparse_config, :libsuitesparseconfig), Ptr{Cvoid})
233+
cnfg = cglobal((:SuiteSparse_config, libsuitesparseconfig), Ptr{Cvoid})
234234
unsafe_store!(cnfg, cglobal(:jl_malloc, Ptr{Cvoid}), 1)
235235
unsafe_store!(cnfg, cglobal(:jl_calloc, Ptr{Cvoid}), 2)
236236
unsafe_store!(cnfg, cglobal(:jl_realloc, Ptr{Cvoid}), 3)

‎src/solvers/umfpack.jl‎

Lines changed: 16 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -261,9 +261,11 @@ has_refinement(F::UmfpackLU) = has_refinement(F.control)
261261
has_refinement(control::AbstractVector) = control[JL_UMFPACK_IRSTEP] > 0
262262

263263
# auto magick resize, should this only expand and not shrink?
264-
getworkspace(F::UmfpackLU) = @lock F.lock begin
265-
return resize!(F.workspace, F, has_refinement(F); expand_only=true)
264+
function getworkspace(F::UmfpackLU)
265+
@lock F.lock begin
266+
return resize!(F.workspace, F, has_refinement(F); expand_only = true)
266267
end
268+
end
267269

268270
UmfpackWS(F::UmfpackLU{Tv, Ti}, refinement::Bool=has_refinement(F)) where {Tv, Ti} = UmfpackWS(
269271
Vector{Ti}(undef, size(F, 2)),
@@ -295,29 +297,12 @@ Base.copy(F::T, ws=UmfpackWS(F)) where {T <: ATLU} =
295297

296298
Base.transpose(F::UmfpackLU) = TransposeFactorization(F)
297299

298-
function Base.lock(f::Function, F::UmfpackLU)
299-
lock(F)
300-
try
301-
f()
302-
finally
303-
unlock(F)
304-
end
305-
end
306-
Base.lock(F::UmfpackLU) = if !trylock(F.lock)
307-
@info """waiting for UmfpackLU's lock, it's safe to ignore this message.
308-
see the documentation for Umfpack""" maxlog = 1
309-
lock(F.lock)
310-
end
311-
312-
@inline Base.trylock(F::UmfpackLU) = trylock(F.lock)
313-
@inline Base.unlock(F::UmfpackLU) = unlock(F.lock)
314-
315300
show_umf_ctrl(F::UmfpackLU, level::Real=2.0) =
316-
@lock F show_umf_ctrl(F.control, level)
301+
@lock F.lock show_umf_ctrl(F.control, level)
317302

318303

319304
show_umf_info(F::UmfpackLU, level::Real=2.0) =
320-
@lock F show_umf_info(F.control, F.info, level)
305+
@lock F.lock show_umf_info(F.control, F.info, level)
321306

322307

323308
"""
@@ -584,7 +569,7 @@ for itype in UmfpackIndexTypes
584569
@eval begin
585570
function umfpack_symbolic!(U::UmfpackLU{Float64,$itype}, q::Union{Nothing, StridedVector{$itype}})
586571
_isnotnull(U.symbolic) && return U
587-
@lock U begin
572+
@lock U.lock begin
588573
tmp = Ref{Ptr{Cvoid}}(C_NULL)
589574
if q === nothing
590575
@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
599584
end
600585
function umfpack_symbolic!(U::UmfpackLU{ComplexF64,$itype}, q::Union{Nothing, StridedVector{$itype}})
601586
_isnotnull(U.symbolic) && return U
602-
@lock U begin
587+
@lock U.lock begin
603588
tmp = Ref{Ptr{Cvoid}}(C_NULL)
604589
if q === nothing
605590
@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
613598
return U
614599
end
615600
function umfpack_numeric!(U::UmfpackLU{Float64,$itype}; reuse_numeric=true, q=nothing)
616-
@lock U begin
601+
@lock U.lock begin
617602
(reuse_numeric && _isnotnull(U.numeric)) && return U
618603
if _isnull(U.symbolic)
619604
umfpack_symbolic!(U, q)
@@ -629,7 +614,7 @@ for itype in UmfpackIndexTypes
629614
return U
630615
end
631616
function umfpack_numeric!(U::UmfpackLU{ComplexF64,$itype}; reuse_numeric=true, q=nothing)
632-
@lock U begin
617+
@lock U.lock begin
633618
(reuse_numeric && _isnotnull(U.numeric)) && return U
634619
_isnull(U.symbolic) && umfpack_symbolic!(U, q)
635620
tmp = Ref{Ptr{Cvoid}}(C_NULL)
@@ -658,7 +643,7 @@ for itype in UmfpackIndexTypes
658643
if workspace_W_size(lu) > length(workspace.W)
659644
throw(ArguementError("W should be larger than `workspace_W_size(Af)`"))
660645
end
661-
@lock lu begin
646+
@lock lu.lock begin
662647
umfpack_numeric!(lu)
663648
(size(b, 1) == lu.m) && (size(b) == size(x)) || throw(DimensionMismatch())
664649

@@ -683,7 +668,7 @@ for itype in UmfpackIndexTypes
683668
if workspace_W_size(lu) > length(workspace.W)
684669
throw(ArguementError("W should be larger than `workspace_W_size(Af)`"))
685670
end
686-
@lock lu begin
671+
@lock lu.lock begin
687672
umfpack_numeric!(lu)
688673
(size(b, 1) == lu.m) && (size(b) == size(x)) || throw(DimensionMismatch())
689674
@isok $wsol_c(typ, lu.colptr, lu.rowval, lu.nzval, C_NULL, x, C_NULL, b,
@@ -693,14 +678,14 @@ for itype in UmfpackIndexTypes
693678
end
694679
function det(lu::UmfpackLU{Float64,$itype})
695680
mx = Ref{Float64}(zero(Float64))
696-
@lock lu @isok($det_r(mx, C_NULL, lu.numeric, lu.info))
681+
@lock lu.lock @isok($det_r(mx, C_NULL, lu.numeric, lu.info))
697682
mx[]
698683
end
699684

700685
function det(lu::UmfpackLU{ComplexF64,$itype})
701686
mx = Ref{Float64}(zero(Float64))
702687
mz = Ref{Float64}(zero(Float64))
703-
@lock lu @isok($det_z(mx, mz, C_NULL, lu.numeric, lu.info))
688+
@lock lu.lock @isok($det_z(mx, mz, C_NULL, lu.numeric, lu.info))
704689
complex(mx[], mz[])
705690
end
706691
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
10161001

10171002
_report_symbolic = Symbol(umf_nm("report_symbolic", Tv, Ti))
10181003
@eval umfpack_report_symbolic(lu::UmfpackLU{$Tv,$Ti}, level::Real=4; q=nothing) =
1019-
@lock lu begin
1004+
@lock lu.lock begin
10201005
umfpack_symbolic!(lu, q)
10211006
old_prl = lu.control[JL_UMFPACK_PRL]
10221007
lu.control[JL_UMFPACK_PRL] = level
@@ -1026,7 +1011,7 @@ for Tv in (:Float64, :ComplexF64), Ti in UmfpackIndexTypes
10261011
end
10271012
_report_numeric = Symbol(umf_nm("report_numeric", Tv, Ti))
10281013
@eval umfpack_report_numeric(lu::UmfpackLU{$Tv,$Ti}, level::Real=4; q=nothing) =
1029-
@lock lu begin
1014+
@lock lu.lock begin
10301015
umfpack_numeric!(lu; q)
10311016
old_prl = lu.control[JL_UMFPACK_PRL]
10321017
lu.control[JL_UMFPACK_PRL] = level

‎src/sparsematrix.jl‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3984,9 +3984,9 @@ function _blockdiag(::Type{Tv}, ::Type{Ti}, X::AbstractSparseMatrixCSC...) where
39843984
end
39853985

39863986
## Structure query functions
3987-
issymmetric(A::AbstractSparseMatrixCSC) = is_hermsym(A, identity)
3987+
issymmetric(A::AbstractSparseMatrixCSC) = is_hermsym(A, transpose)
39883988

3989-
ishermitian(A::AbstractSparseMatrixCSC) = is_hermsym(A, conj)
3989+
ishermitian(A::AbstractSparseMatrixCSC) = is_hermsym(A, adjoint)
39903990

39913991
function is_hermsym(A::AbstractSparseMatrixCSC, check::Function)
39923992
m, n = size(A)
@@ -4022,6 +4022,12 @@ function is_hermsym(A::AbstractSparseMatrixCSC, check::Function)
40224022
return false
40234023
end
40244024
else
4025+
# if nzrange(A, row) is empty, then A[:, row] is all zeros.
4026+
# Specifically, A[col, row] is zero.
4027+
# However, we know at this point that A[row, col] is not zero
4028+
# This means that the matrix is not symmetric
4029+
isempty(nzrange(A, row)) && return false
4030+
40254031
offset = tracker[row]
40264032

40274033
# If the matrix is unsymmetric, there might not exist

‎test/linalg.jl‎

Lines changed: 31 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,13 @@ end
224224
end
225225
end
226226

227+
@testset "Dense times symmetric/Hermitian sparse matrix multiplication" begin
228+
A = [1 3; 2 4]
229+
As = sparse(A)
230+
B = [1 1; 1 1]
231+
@test mul!(copy(B), B, Hermitian(A), true, true) == mul!(copy(B), B, Hermitian(As), true, true)
232+
end
233+
227234
@testset "UniformScaling" begin
228235
local A = sprandn(10, 10, 0.5)
229236
@test A + I == Array(A) + I
@@ -395,6 +402,27 @@ end
395402
@test issymmetric(sparse([1 0; 1 0])) == false
396403
@test issymmetric(sparse([0 1; 1 0])) == true
397404
@test issymmetric(sparse([1 1; 1 0])) == true
405+
406+
# test some non-trivial cases
407+
local S
408+
@testset "random matrices" begin
409+
for sparsity in (0.1, 0.01, 0.0)
410+
S = sparse(Symmetric(sprand(20, 20, sparsity)))
411+
@test issymmetric(S)
412+
@test ishermitian(S)
413+
S = sparse(Symmetric(sprand(ComplexF64, 20, 20, sparsity)))
414+
@test issymmetric(S)
415+
@test !ishermitian(S) || isreal(S)
416+
S = sparse(Hermitian(sprand(ComplexF64, 20, 20, sparsity)))
417+
@test ishermitian(S)
418+
@test !issymmetric(S) || isreal(S)
419+
end
420+
end
421+
422+
@testset "issue #605" begin
423+
S = sparse([2, 3, 1], [1, 1, 3], [1, 1, 1], 3, 3)
424+
@test !issymmetric(S)
425+
end
398426
end
399427

400428
@testset "rotations" begin
@@ -824,11 +852,12 @@ end
824852
@test dot(x, A, y) ≈ dot(x, Av, y)
825853
end
826854

827-
for (T, trans) in ((Float64, Symmetric), (ComplexF64, Hermitian)), uplo in (:U, :L)
855+
for (T, trans) in ((Float64, Symmetric), (ComplexF64, Symmetric), (ComplexF64, Hermitian)), uplo in (:U, :L)
828856
B = sprandn(T, 10, 10, 0.2)
829857
x = sprandn(T, 10, 0.4)
830858
S = trans(B'B, uplo)
831-
@test dot(x, S, x) ≈ dot(Vector(x), S, Vector(x)) ≈ dot(Vector(x), Matrix(S), Vector(x))
859+
Sd = trans(Matrix(B'B), uplo)
860+
@test dot(x, S, x) ≈ dot(x, Sd, x) ≈ dot(Vector(x), S, Vector(x)) ≈ dot(Vector(x), Sd, Vector(x))
832861
end
833862
end
834863

‎test/umfpack.jl‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -121,8 +121,6 @@ end
121121
umfpack_report(Af)
122122
Af1 = lu!(copy(Af))
123123
umfpack_report(Af1)
124-
@test trylock(Af)
125-
@test trylock(Af1)
126124
end
127125

128126
@testset "test similar" begin

0 commit comments

Comments
 (0)