add specific Adjoint{StridedMatrix} * SparseVector method (#32403)

add complex and adjoint tests

Author: dkarrasch

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

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