Fix matrix multiplication with non-commutative elements (#321)

* Fix matrix multiplication with array elements

* fix fillmatrix * stridedvec

* Fix for adjoint

* Tests with non-zero beta

* Bump version to v1.9.3

* Fix 5-term mul order and test against Quaternions

* Move imports to one line

---------

Co-authored-by: Sheehan Olver <[email protected]>

Author: jishnub

Project.toml

@@ -26,6 +26,7 @@ Documenter = "1"
 Infinities = "0.1"
 LinearAlgebra = "1.6"
 PDMats = "0.11.17"
+Quaternions = "0.7"
 Random = "1.6"
 ReverseDiff = "1"
 SparseArrays = "1.6"
@@ -40,11 +41,12 @@ Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test", "Base64", "Infinities", "PDMats", "ReverseDiff", "SparseArrays", "StaticArrays", "Statistics", "Documenter"]
+test = ["Aqua", "Test", "Base64", "Infinities", "PDMats", "ReverseDiff", "SparseArrays", "StaticArrays", "Statistics", "Quaternions", "Documenter"]

src/fillalgebra.jl

@@ -150,60 +150,56 @@ end
 function mul!(y::AbstractVector, A::AbstractFillMatrix, b::AbstractFillVector, alpha::Number, beta::Number)
     check_matmul_sizes(y, A, b)
 
-    αAb = alpha * getindex_value(A) * getindex_value(b) * length(b)
+    Abα = Ref(getindex_value(A) * getindex_value(b) * alpha * length(b))
 
     if iszero(beta)
-        y .= αAb
+        y .= Abα
     else
-        y .= αAb .+ beta .* y
+        y .= Abα .+ y .* beta
     end
     y
 end
 
 function mul!(y::StridedVector, A::StridedMatrix, b::AbstractFillVector, alpha::Number, beta::Number)
     check_matmul_sizes(y, A, b)
 
-    αb = alpha * getindex_value(b)
+    bα = Ref(getindex_value(b) * alpha)
 
     if iszero(beta)
-        y .= zero(eltype(y))
-        for col in eachcol(A)
-            y .+= αb .* col
-        end
+        y .= Ref(zero(eltype(y)))
     else
-        lmul!(beta, y)
-        for col in eachcol(A)
-            y .+= αb .* col
-        end
+        rmul!(y, beta)
+    end
+    for Acol in eachcol(A)
+        @. y += Acol * bα
     end
     y
 end
 
 function mul!(y::StridedVector, A::AbstractFillMatrix, b::StridedVector, alpha::Number, beta::Number)
     check_matmul_sizes(y, A, b)
 
-    αA = alpha * getindex_value(A)
+    Abα = Ref(getindex_value(A) * sum(b) * alpha)
 
     if iszero(beta)
-        y .= αA .* sum(b)
+        y .= Abα
     else
-        y .= αA .* sum(b) .+ beta .* y
+        y .= Abα .+ y .* beta
     end
     y
 end
 
-function _mul_adjtrans!(y::AbstractVector, A::AbstractMatrix, b::AbstractVector, alpha, beta, f)
-    α = alpha * getindex_value(b)
-
+function _mul_adjtrans!(y::AbstractVector, A::AbstractMatrix, b::AbstractFillVector, alpha, beta, f)
+    bα = getindex_value(b) * alpha
     At = f(A)
 
     if iszero(beta)
-        for (ind, col) in zip(eachindex(y), eachcol(At))
-            y[ind] = α .* f(sum(col))
+        for (ind, Atcol) in zip(eachindex(y), eachcol(At))
+            y[ind] = f(sum(Atcol)) * bα
         end
     else
-        for (ind, col) in zip(eachindex(y), eachcol(At))
-            y[ind] = α .* f(sum(col)) .+ beta .* y[ind]
+        for (ind, Atcol) in zip(eachindex(y), eachcol(At))
+            y[ind] = f(sum(Atcol)) * bα .+ y[ind] .* beta
         end
     end
     y
@@ -218,11 +214,11 @@ end
 
 function mul!(C::AbstractMatrix, A::AbstractFillMatrix, B::AbstractFillMatrix, alpha::Number, beta::Number)
     check_matmul_sizes(C, A, B)
-    αAB = alpha * getindex_value(A) * getindex_value(B) * size(B,1)
+    ABα = getindex_value(A) * getindex_value(B) * alpha * size(B,1)
     if iszero(beta)
-        C .= αAB
+        C .= ABα
     else
-        C .= αAB .+ beta .* C
+        C .= ABα .+ C .* beta
     end
     C
 end
@@ -248,7 +244,7 @@ _firstcol(C::Union{Adjoint, Transpose}) = view(parent(C), 1, :)
 function _mulfill!(C, A, B::AbstractFillMatrix, alpha, beta)
     check_matmul_sizes(C, A, B)
     if iszero(size(B,2))
-        return lmul!(beta, C)
+        return rmul!(C, beta)
     end
     mul!(_firstcol(C), A, view(B, :, 1), alpha, beta)
     copyfirstcol!(C)

test/runtests.jl

@@ -1,4 +1,4 @@
-using FillArrays, LinearAlgebra, PDMats, SparseArrays, StaticArrays, ReverseDiff, Random, Base64, Test, Statistics
+using FillArrays, LinearAlgebra, PDMats, SparseArrays, StaticArrays, ReverseDiff, Random, Base64, Test, Statistics, Quaternions
 import FillArrays: AbstractFill, RectDiagonal, SquareEye
 
 using Documenter
@@ -1558,11 +1558,25 @@ end
         Z = Zeros(SMatrix{2,3,Float64,6}, 2)
         @test Z' * D == Array(Z)' * D
 
+        S = SMatrix{2,3}(1:6)
+        A = reshape([S,2S,3S,4S],2,2)
+        F = Fill(S',2,2)
+        @test A * F == A * fill(S',size(F))
+        @test mul!(A * F, A, F, 2, 1) == 3 * A * fill(S',size(F))
+        @test F * A == fill(S',size(F)) * A
+        @test mul!(F * A, F, A, 2, 1) == 3 * fill(S',size(F)) * A
+
         # doubly nested
         A = [[[1,2]]]'
         Z = Zeros(SMatrix{1,1,SMatrix{2,2,Int,4},1},1)
         Z2 = zeros(SMatrix{1,1,SMatrix{2,2,Int,4},1},1)
         @test A * Z == A * Z2
+
+        x = [1 2 3; 4 5 6]
+        A = reshape([x,2x,3x,4x],2,2)
+        F = Fill(x,2,2)
+        @test A' * F == A' * fill(x,size(F))
+        @test mul!(A' * F, A', F, 2, 1) == 3 * A' * fill(x,size(F))
     end
 
     for W in (zeros(3,4), @MMatrix zeros(3,4))
@@ -1697,6 +1711,40 @@ end
         @test adjoint(A)*fillmat ≈ adjoint(A)*Array(fillmat)
     end
 
+    @testset "non-commutative" begin
+        A = Fill(quat(rand(4)...), 2, 2)
+        M = Array(A)
+        α, β = quat(0,1,1,0), quat(1,0,0,1)
+        @testset "matvec" begin
+            f = Fill(quat(rand(4)...), size(A,2))
+            v = Array(f)
+            D = copy(v)
+            exp_res = M * v * α + D * β
+            @test mul!(copy(D), A, f, α, β) ≈ mul!(copy(D), M, v, α, β) ≈ exp_res
+            @test mul!(copy(D), M, f, α, β) ≈ mul!(copy(D), M, v, α, β) ≈ exp_res
+            @test mul!(copy(D), A, v, α, β) ≈ mul!(copy(D), M, v, α, β) ≈ exp_res
+
+            @test mul!(copy(D), M', f, α, β) ≈ mul!(copy(D), M', v, α, β) ≈ M' * v * α + D * β
+            @test mul!(copy(D), transpose(M), f, α, β) ≈ mul!(copy(D), transpose(M), v, α, β) ≈ transpose(M) * v * α + D * β
+        end
+
+        @testset "matmat" begin
+            B = Fill(quat(rand(4)...), 2, 2)
+            N = Array(B)
+            D = copy(N)
+            exp_res = M * N * α + D * β
+            @test mul!(copy(D), A, B, α, β) ≈ mul!(copy(D), M, N, α, β) ≈ exp_res
+            @test mul!(copy(D), M, B, α, β) ≈ mul!(copy(D), M, N, α, β) ≈ exp_res
+            @test mul!(copy(D), A, N, α, β) ≈ mul!(copy(D), M, N, α, β) ≈ exp_res
+
+            @test mul!(copy(D), M', B, α, β) ≈ mul!(copy(D), M', N, α, β) ≈ M' * N * α + D * β
+            @test mul!(copy(D), transpose(M), B, α, β) ≈ mul!(copy(D), transpose(M), N, α, β) ≈ transpose(M) * N * α + D * β
+
+            @test mul!(copy(D), A, N', α, β) ≈ mul!(copy(D), M, N', α, β) ≈ M * N' * α + D * β
+            @test mul!(copy(D), A, transpose(N), α, β) ≈ mul!(copy(D), M, transpose(N), α, β) ≈ M * transpose(N) * α + D * β
+        end
+    end
+
     @testset "ambiguities" begin
         UT33 = UpperTriangular(ones(3,3))
         UT11 = UpperTriangular(ones(1,1))