From ddcfb041ee4ad28b55a1472b1aed3ae60826d61b Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 24 Jun 2019 22:10:02 +0200 Subject: [PATCH] add specific Adjoint{StridedMatrix} * SparseVector method add complex and adjoint tests --- stdlib/SparseArrays/src/sparsevector.jl | 72 ++++++++++++++++++------ stdlib/SparseArrays/test/sparsevector.jl | 58 ++++++++++++------- 2 files changed, 93 insertions(+), 37 deletions(-) diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index 5cadbe1524ca1..28201f4f49219 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -1494,7 +1494,7 @@ function (*)(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx} require_one_based_indexing(A, x) m, n = size(A) length(x) == n || throw(DimensionMismatch()) - Ty = promote_type(Ta, Tx) + Ty = promote_op(matprod, Ta, Tx) y = Vector{Ty}(undef, m) mul!(y, A, x) end @@ -1531,23 +1531,62 @@ end function *(transA::Transpose{<:Any,<:StridedMatrix{Ta}}, x::AbstractSparseVector{Tx}) where {Ta,Tx} require_one_based_indexing(transA, x) - A = transA.parent - m, n = size(A) - length(x) == m || throw(DimensionMismatch()) - Ty = promote_type(Ta, Tx) - y = Vector{Ty}(undef, n) - mul!(y, transpose(A), x) + m, n = size(transA) + length(x) == n || throw(DimensionMismatch()) + Ty = promote_op(matprod, Ta, Tx) + y = Vector{Ty}(undef, m) + mul!(y, transA, x) end mul!(y::AbstractVector{Ty}, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = - (A = transA.parent; mul!(y, transpose(A), x, one(Tx), zero(Ty))) + mul!(y, transA, x, one(Tx), zero(Ty)) function mul!(y::AbstractVector, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector, α::Number, β::Number) + require_one_based_indexing(y, transA, x) + m, n = size(transA) + length(x) == n && length(y) == m || throw(DimensionMismatch()) + m == 0 && return y + if β != one(β) + β == zero(β) ? fill!(y, zero(eltype(y))) : rmul!(y, β) + end + α == zero(α) && return y + + xnzind = nonzeroinds(x) + xnzval = nonzeros(x) + _nnz = length(xnzind) + _nnz == 0 && return y + A = transA.parent - require_one_based_indexing(y, A, x) - m, n = size(A) - length(x) == m && length(y) == n || throw(DimensionMismatch()) - n == 0 && return y + Ty = promote_op(matprod, eltype(A), eltype(x)) + @inbounds for j = 1:m + s = zero(Ty) + for i = 1:_nnz + s += transpose(A[xnzind[i], j]) * xnzval[i] + end + y[j] += s * α + end + return y +end + +# * and mul!(C, adjoint(A), B) + +function *(adjA::Adjoint{<:Any,<:StridedMatrix{Ta}}, x::AbstractSparseVector{Tx}) where {Ta,Tx} + require_one_based_indexing(adjA, x) + m, n = size(adjA) + length(x) == n || throw(DimensionMismatch()) + Ty = promote_op(matprod, Ta, Tx) + y = Vector{Ty}(undef, m) + mul!(y, adjA, x) +end + +mul!(y::AbstractVector{Ty}, adjA::Adjoint{<:Any,<:StridedMatrix}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = + mul!(y, adjA, x, one(Tx), zero(Ty)) + +function mul!(y::AbstractVector, adjA::Adjoint{<:Any,<:StridedMatrix}, x::AbstractSparseVector, α::Number, β::Number) + require_one_based_indexing(y, adjA, x) + m, n = size(adjA) + length(x) == n && length(y) == m || throw(DimensionMismatch()) + m == 0 && return y if β != one(β) β == zero(β) ? fill!(y, zero(eltype(y))) : rmul!(y, β) end @@ -1558,11 +1597,12 @@ function mul!(y::AbstractVector, transA::Transpose{<:Any,<:StridedMatrix}, x::Ab _nnz = length(xnzind) _nnz == 0 && return y - s0 = zero(eltype(A)) * zero(eltype(x)) - @inbounds for j = 1:n - s = zero(s0) + A = adjA.parent + Ty = promote_op(matprod, eltype(A), eltype(x)) + @inbounds for j = 1:m + s = zero(Ty) for i = 1:_nnz - s += A[xnzind[i], j] * xnzval[i] + s += adjoint(A[xnzind[i], j]) * xnzval[i] end y[j] += s * α end diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index 4adfc3398edd8..1715d5a49ca56 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -841,30 +841,46 @@ end @testset "BLAS Level-2" begin @testset "dense A * sparse x -> dense y" begin - let A = randn(9, 16), x = sprand(16, 0.7) - xf = Array(x) - for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0] - y = rand(9) - rr = α*A*xf + β*y - @test mul!(y, A, x, α, β) === y - @test y ≈ rr + for TA in (Float64, ComplexF64), Tx in (Float64, ComplexF64) + T = Base.promote_op(LinearAlgebra.matprod, TA, Tx) + let A = randn(TA, 9, 16), x = sprand(Tx, 16, 0.7) + xf = Array(x) + for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0] + y = rand(T, 9) + rr = α*A*xf + β*y + @test mul!(y, A, x, α, β) === y + @test y ≈ rr + end + y = A*x + @test isa(y, Vector{T}) + @test A*x ≈ A*xf end - y = A*x - @test isa(y, Vector{Float64}) - @test A*x ≈ A*xf - end - let A = randn(16, 9), x = sprand(16, 0.7) - xf = Array(x) - for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0] - y = rand(9) - rr = α*A'xf + β*y - @test mul!(y, transpose(A), x, α, β) === y - @test y ≈ rr + let A = randn(TA, 16, 9), x = sprand(Tx, 16, 0.7) + xf = Array(x) + for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0] + y = rand(T, 9) + rr = α*transpose(A)*xf + β*y + @test mul!(y, transpose(A), x, α, β) === y + @test y ≈ rr + end + y = *(transpose(A), x) + @test isa(y, Vector{T}) + @test y ≈ *(transpose(A), xf) + end + + let A = randn(TA, 16, 9), x = sprand(Tx, 16, 0.7) + xf = Array(x) + for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0] + y = rand(T, 9) + rr = α*A'xf + β*y + @test mul!(y, adjoint(A), x, α, β) === y + @test y ≈ rr + end + y = *(adjoint(A), x) + @test isa(y, Vector{T}) + @test y ≈ *(adjoint(A), xf) end - y = *(transpose(A), x) - @test isa(y, Vector{Float64}) - @test y ≈ *(transpose(A), xf) end end @testset "sparse A * sparse x -> dense y" begin