[SparseArrays] Respect order of mul in (l)mul!(::Diagonal,::Sparse) (#…

…30163)

* order of mul in (l)mul!(::Diagonal,::Sparse)

add tests for non-commutative mul

* Create Quaternions.jl

remove quaternions from generic tests

stdlib/LinearAlgebra/test/generic.jl

@@ -3,33 +3,10 @@
 module TestGeneric
 
 using Test, LinearAlgebra, Random
-import Base: -, *, /, \
-
-# A custom Quaternion type with minimal defined interface and methods.
-# Used to test mul and mul! methods to show non-commutativity.
-struct Quaternion{T<:Real} <: Number
-    s::T
-    v1::T
-    v2::T
-    v3::T
-end
-Quaternion(s::Real, v1::Real, v2::Real, v3::Real) = Quaternion(promote(s, v1, v2, v3)...)
-Base.abs2(q::Quaternion) = q.s*q.s + q.v1*q.v1 + q.v2*q.v2 + q.v3*q.v3
-Base.abs(q::Quaternion) = sqrt(abs2(q))
-Base.real(::Type{Quaternion{T}}) where {T} = T
-Base.conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3)
-Base.isfinite(q::Quaternion) = isfinite(q.s) & isfinite(q.v1) & isfinite(q.v2) & isfinite(q.v3)
-
-(-)(ql::Quaternion, qr::Quaternion) =
-    Quaternion(ql.s - qr.s, ql.v1 - qr.v1, ql.v2 - qr.v2, ql.v3 - qr.v3)
-(*)(q::Quaternion, w::Quaternion) = Quaternion(q.s*w.s - q.v1*w.v1 - q.v2*w.v2 - q.v3*w.v3,
-                                               q.s*w.v1 + q.v1*w.s + q.v2*w.v3 - q.v3*w.v2,
-                                               q.s*w.v2 - q.v1*w.v3 + q.v2*w.s + q.v3*w.v1,
-                                               q.s*w.v3 + q.v1*w.v2 - q.v2*w.v1 + q.v3*w.s)
-(*)(q::Quaternion, r::Real) = Quaternion(q.s*r, q.v1*r, q.v2*r, q.v3*r)
-(*)(q::Quaternion, b::Bool) = b * q # remove method ambiguity
-(/)(q::Quaternion, w::Quaternion) = q * conj(w) * (1.0 / abs2(w))
-(\)(q::Quaternion, w::Quaternion) = conj(q) * w * (1.0 / abs2(q))
+
+const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+isdefined(Main, :Quaternions) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Quaternions.jl"))
+using .Main.Quaternions
 
 Random.seed!(123)
 

stdlib/SparseArrays/src/linalg.jl

@@ -1202,7 +1202,7 @@ function mul!(C::SparseMatrixCSC, D::Diagonal{T, <:Vector}, A::SparseMatrixCSC)
     Arowval = A.rowval
     resize!(Cnzval, length(Anzval))
     for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
-        @inbounds Cnzval[p] = Anzval[p] * b[Arowval[p]]
+        @inbounds Cnzval[p] = b[Arowval[p]] * Anzval[p]
     end
     C
 end
@@ -1238,7 +1238,7 @@ function rmul!(A::SparseMatrixCSC, D::Diagonal)
     (n == size(D, 1)) || throw(DimensionMismatch())
     Anzval = A.nzval
     @inbounds for col = 1:n, p = A.colptr[col]:(A.colptr[col + 1] - 1)
-         Anzval[p] *= D.diag[col]
+         Anzval[p] = Anzval[p] * D.diag[col]
     end
     return A
 end
@@ -1249,7 +1249,7 @@ function lmul!(D::Diagonal, A::SparseMatrixCSC)
     Anzval = A.nzval
     Arowval = A.rowval
     @inbounds for col = 1:n, p = A.colptr[col]:(A.colptr[col + 1] - 1)
-        Anzval[p] *= D.diag[Arowval[p]]
+        Anzval[p] = D.diag[Arowval[p]] * Anzval[p]
     end
     return A
 end

stdlib/SparseArrays/test/sparse.jl

@@ -365,6 +365,10 @@ end
     @test_throws DimensionMismatch dot(sprand(5,5,0.2),sprand(5,6,0.2))
 end
 
+const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+isdefined(Main, :Quaternions) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Quaternions.jl"))
+using .Main.Quaternions
+
 sA = sprandn(3, 7, 0.5)
 sC = similar(sA)
 dA = Array(sA)
@@ -403,6 +407,26 @@ dA = Array(sA)
         @test_throws DimensionMismatch rdiv!(copy(sAt), Diagonal(fill(1., length(b)+1)))
         @test_throws LinearAlgebra.SingularException rdiv!(copy(sAt), Diagonal(zeros(length(b))))
     end
+
+    @testset "non-commutative multiplication" begin
+        # non-commutative multiplication
+        Avals = Quaternion.(randn(10), randn(10), randn(10), randn(10))
+        sA = sparse(rand(1:3, 10), rand(1:7, 10), Avals, 3, 7)
+        sC = copy(sA)
+        dA = Array(sA)
+
+        b = Quaternion.(randn(7), randn(7), randn(7), randn(7))
+        D = Diagonal(b)
+        @test Array(sA * D) ≈ dA * D
+        @test rmul!(copy(sA), D) ≈ dA * D
+        @test mul!(sC, copy(sA), D) ≈ dA * D
+
+        b = Quaternion.(randn(3), randn(3), randn(3), randn(3))
+        D = Diagonal(b)
+        @test Array(D * sA) ≈ D * dA
+        @test lmul!(D, copy(sA)) ≈ D * dA
+        @test mul!(sC, D, copy(sA)) ≈ D * dA
+    end
 end
 
 @testset "copyto!" begin

test/testhelpers/Quaternions.jl

@@ -0,0 +1,34 @@
+module Quaternions
+
+export Quaternion
+
+# A custom Quaternion type with minimal defined interface and methods.
+# Used to test mul and mul! methods to show non-commutativity.
+struct Quaternion{T<:Real} <: Number
+    s::T
+    v1::T
+    v2::T
+    v3::T
+end
+Quaternion(s::Real, v1::Real, v2::Real, v3::Real) = Quaternion(promote(s, v1, v2, v3)...)
+Base.abs2(q::Quaternion) = q.s*q.s + q.v1*q.v1 + q.v2*q.v2 + q.v3*q.v3
+Base.abs(q::Quaternion) = sqrt(abs2(q))
+Base.real(::Type{Quaternion{T}}) where {T} = T
+Base.conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3)
+Base.isfinite(q::Quaternion) = isfinite(q.s) & isfinite(q.v1) & isfinite(q.v2) & isfinite(q.v3)
+Base.zero(::Type{Quaternion{T}}) where T = Quaternion{T}(zero(T), zero(T), zero(T), zero(T))
+
+Base.:(+)(ql::Quaternion, qr::Quaternion) =
+ Quaternion(ql.s + qr.s, ql.v1 + qr.v1, ql.v2 + qr.v2, ql.v3 + qr.v3)
+Base.:(-)(ql::Quaternion, qr::Quaternion) =
+ Quaternion(ql.s - qr.s, ql.v1 - qr.v1, ql.v2 - qr.v2, ql.v3 - qr.v3)
+Base.:(*)(q::Quaternion, w::Quaternion) = Quaternion(q.s*w.s - q.v1*w.v1 - q.v2*w.v2 - q.v3*w.v3,
+                                            q.s*w.v1 + q.v1*w.s + q.v2*w.v3 - q.v3*w.v2,
+                                            q.s*w.v2 - q.v1*w.v3 + q.v2*w.s + q.v3*w.v1,
+                                            q.s*w.v3 + q.v1*w.v2 - q.v2*w.v1 + q.v3*w.s)
+Base.:(*)(q::Quaternion, r::Real) = Quaternion(q.s*r, q.v1*r, q.v2*r, q.v3*r)
+Base.:(*)(q::Quaternion, b::Bool) = b * q # remove method ambiguity
+Base.:(/)(q::Quaternion, w::Quaternion) = q * conj(w) * (1.0 / abs2(w))
+Base.:(\)(q::Quaternion, w::Quaternion) = conj(q) * w * (1.0 / abs2(q))
+
+end