performance improvement for Symmetric/Hermitian of sparse matrix time…

…s dense vector (JuliaLang#30018)

stdlib/SparseArrays/src/linalg.jl

@@ -677,6 +677,72 @@ end
 
 ## end of triangular
 
+## symmetric/Hermitian wrappers of sparse matrices
+"""
+ SparseMatrixCSCSymmHerm
+
+`Symmetric` or `Hermitian` of a `SparseMatrixCSC` or `SparseMatrixCSCView`.
+"""
+const SparseMatrixCSCSymmHerm{Tv,Ti} = Union{Symmetric{Tv,<:SparseMatrixCSCUnion{Tv,Ti}},
+ Hermitian{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}}
+
+# y .= A * x
+mul!(y::StridedVecOrMat, A::SparseMatrixCSCSymmHerm, x::StridedVecOrMat) = mul!(y,A,x,1,0)
+
+# C .= α * A * B + β * C
+function mul!(C::StridedVecOrMat{T}, sA::SparseMatrixCSCSymmHerm, B::StridedVecOrMat,
+ α::Number, β::Number) where T
+
+ fuplo = sA.uplo == 'U' ? nzrangeup : nzrangelo
+ _mul!(fuplo, C, sA, B, T(α), T(β))
+end
+
+function _mul!(nzrang::Function, C::StridedVecOrMat{T}, sA, B, α, β) where T
+ A = sA.data
+ n = size(A, 2)
+ m = size(B, 2)
+ n == size(B, 1) == size(C, 1) && m == size(C, 2) || throw(DimensionMismatch())
+ rv = rowvals(A)
+ nzv = nonzeros(A)
+ let z = T(0), sumcol=z, αxj=z, aarc=z, α = α
+ if β != 1
+ β != 0 ? rmul!(C, β) : fill!(C, z)
+ end
+ @inbounds for k = 1:m
+ for col = 1:n
+ αxj = B[col,k] * α
+ sumcol = z
+ for j = nzrang(A, col)
+ row = rv[j]
+ aarc = nzv[j]
+ if row == col
+ sumcol += (sA isa Hermitian ? real : identity)(aarc) * B[row,k]
+ else
+ C[row,k] += aarc * αxj
+ sumcol += (sA isa Hermitian ? adjoint : transpose)(aarc) * B[row,k]
+ end
+ end
+ C[col,k] += α * sumcol
+ end
+ end
+ end
+ C
+end
+
+# row range up to and including diagonal
+function nzrangeup(A, i)
+ r = nzrange(A, i); r1 = r.start; r2 = r.stop
+ rv = rowvals(A)
+ @inbounds r2 < r1 || rv[r2] <= i ? r : r1:searchsortedlast(rv, i, r1, r2, Forward)
+end
+# row range from diagonal (included) to end
+function nzrangelo(A, i)
+ r = nzrange(A, i); r1 = r.start; r2 = r.stop
+ rv = rowvals(A)
+ @inbounds r2 < r1 || rv[r1] >= i ? r : searchsortedfirst(rv, i, r1, r2, Forward):r2
+end
+## end of symmetric/Hermitian
+
 \(A::Transpose{<:Real,<:Hermitian{<:Real,<:SparseMatrixCSC}}, B::Vector) = A.parent \ B
 \(A::Transpose{<:Complex,<:Hermitian{<:Complex,<:SparseMatrixCSC}}, B::Vector) = copy(A) \ B
 \(A::Transpose{<:Number,<:Symmetric{<:Number,<:SparseMatrixCSC}}, B::Vector) = A.parent \ B

stdlib/SparseArrays/src/sparsematrix.jl

@@ -101,6 +101,7 @@ julia> nonzeros(A)
 ```
 """
 nonzeros(S::SparseMatrixCSC) = S.nzval
+nonzeros(S::SparseMatrixCSCView) = nonzeros(S.parent)
 
 """
  rowvals(A::SparseMatrixCSC)
@@ -126,6 +127,7 @@ julia> rowvals(A)
 ```
 """
 rowvals(S::SparseMatrixCSC) = S.rowval
+rowvals(S::SparseMatrixCSCView) = rowvals(S.parent)
 
 """
  nzrange(A::SparseMatrixCSC, col::Integer)
@@ -147,6 +149,7 @@ column. In conjunction with [`nonzeros`](@ref) and
  end
 """
 nzrange(S::SparseMatrixCSC, col::Integer) = S.colptr[col]:(S.colptr[col+1]-1)
+nzrange(S::SparseMatrixCSCView, col::Integer) = nzrange(S.parent, S.indices[2][col])
 
 function Base.show(io::IO, ::MIME"text/plain", S::SparseMatrixCSC)
  xnnz = nnz(S)

stdlib/SparseArrays/test/sparse.jl

@@ -2412,6 +2412,40 @@ end
  @test m2.module == SparseArrays
 end
 
+@testset "Symmetric of sparse matrix mul! dense vector" begin
+ rng = Random.MersenneTwister(1)
+ n = 1000
+ p = 0.02
+ q = 1 - sqrt(1-p)
+ Areal = sprandn(rng, n, n, p)
+ Breal = randn(rng, n)
+ Acomplex = sprandn(rng, n, n, q) + sprandn(rng, n, n, q) * im
+ Bcomplex = Breal + randn(rng, n) * im
+ @testset "symmetric/Hermitian sparse multiply with $S($U)" for S in (Symmetric, Hermitian), U in (:U, :L), (A, B) in ((Areal,Breal), (Acomplex,Bcomplex))
+ Asym = S(A, U)
+ As = sparse(Asym) # takes most time
+ @test which(mul!, (typeof(B), typeof(Asym), typeof(B))).module == SparseArrays
+ @test norm(Asym * B - As * B, Inf) <= eps() * n * p * 10
+ end
+end
+
+@testset "Symmetric of view of sparse matrix mul! dense vector" begin
+ rng = Random.MersenneTwister(1)
+ n = 1000
+ p = 0.02
+ q = 1 - sqrt(1-p)
+ Areal = view(sprandn(rng, n, n+10, p), :, 6:n+5)
+ Breal = randn(rng, n)
+ Acomplex = view(sprandn(rng, n, n+10, q) + sprandn(rng, n, n+10, q) * im, :, 6:n+5)
+ Bcomplex = Breal + randn(rng, n) * im
+ @testset "symmetric/Hermitian sparseview multiply with $S($U)" for S in (Symmetric, Hermitian), U in (:U, :L), (A, B) in ((Areal,Breal), (Acomplex,Bcomplex))
+ Asym = S(A, U)
+ As = sparse(Asym) # takes most time
+ @test which(mul!, (typeof(B), typeof(Asym), typeof(B))).module == SparseArrays
+ @test norm(Asym * B - As * B, Inf) <= eps() * n * p * 10
+ end
+end
+
 @testset "sprand" begin
  p=0.3; m=1000; n=2000;
  for s in 1:10