Diagonal-sandwiched triple product for SparseMatrixCSC

Author: jishnub

src/linalg.jl

@@ -205,6 +205,9 @@ const SparseOrTri{Tv,Ti} = Union{SparseMatrixCSCUnion{Tv,Ti},SparseTriangular{Tv
 *(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::SparseOrTri) = spmatmul(copy(A), B)
 *(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(copy(A), copy(B))
 
+(*)(Da::Diagonal, A::Union{SparseMatrixCSCUnion, AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}}, Db::Diagonal) =
+    Da * (A * Db)
+
 # Gustavson's matrix multiplication algorithm revisited.
 # The result rowval vector is already sorted by construction.
 # The auxiliary Vector{Ti} xb is replaced by a Vector{Bool} of same length.

test/linalg.jl

@@ -912,4 +912,20 @@ end
     @test sparse(3I, 4, 5) == sparse(1:4, 1:4, 3, 4, 5)
     @test sparse(3I, 5, 4) == sparse(1:4, 1:4, 3, 5, 4)
 end
+
+@testset "diagonal-sandwiched triple multiplication" begin
+    D = Diagonal(1:4)
+    S = sprand(Int, 4, 4, 0.2)
+    A = Array(S)
+    C = D * S * D
+    @test C isa SparseMatrixCSC
+    @test C ≈ D * A * D
+    C = D * S' * D
+    @test C isa SparseMatrixCSC
+    @test C ≈ D * A' * D
+    C = D * view(S, :, :) * D
+    @test C isa SparseMatrixCSC
+    @test C ≈ D * A * D
+end
+
 end