Fix mat * vec with non-commutative numbers (#413)

* Fix mat * vec with non-commutative numbers

* Update banded * diagonal

* Change some lmul to rmul

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "BandedMatrices"
 uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "1.3"
+version = "1.3.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -24,6 +24,7 @@ GenericLinearAlgebra = "0.3"
 InfiniteArrays = "0.12, 0.13"
 LinearAlgebra = "1.6"
 PrecompileTools = "1"
+Quaternions = "0.7"
 Random = "1.6"
 SparseArrays = "1.6"
 Test = "1.6"
@@ -34,9 +35,10 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
+Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Documenter", "GenericLinearAlgebra", "InfiniteArrays", "Random", "SparseArrays", "Test"]
+test = ["Aqua", "Documenter", "GenericLinearAlgebra", "InfiniteArrays", "Random", "SparseArrays", "Test", "Quaternions"]

src/BandedMatrices.jl

@@ -29,7 +29,7 @@ import ArrayLayouts: MemoryLayout, transposelayout, triangulardata,
                     triangularlayout, MatLdivVec, hermitianlayout, hermitiandata,
                     materialize, materialize!, BlasMatMulMatAdd, BlasMatMulVecAdd, BlasMatLmulVec, BlasMatLdivVec,
                     colsupport, rowsupport, symmetricuplo, MatMulMatAdd, MatMulVecAdd,
-                    sublayout, sub_materialize, _fill_lmul!, _copy_oftype,
+                    sublayout, sub_materialize, _copy_oftype, zero!,
                     reflector!, reflectorApply!, _copyto!, checkdimensions,
                     _qr!, _qr, _lu!, _lu, _factorize, AbstractTridiagonalLayout, TridiagonalLayout,
                     BidiagonalLayout, bidiagonaluplo, diagonaldata, supdiagonaldata, subdiagonaldata, copymutable_oftype_layout, dualadjoint

src/generic/matmul.jl

@@ -30,7 +30,7 @@ banded_gbmv!(tA, α, A, x, β, y) =
     =#
     length(y) == 0 && return y
     if length(x) == 0
-        _fill_lmul!(β, y)
+        _fill_rmul!(y, β)
     else
         xc = Base.unalias(y, x)
         banded_gbmv!(tA, α, A, xc, β, y)
@@ -42,15 +42,15 @@ function _banded_muladd!(α, A, x::AbstractVector, β, y)
     m, n = size(A)
     l, u = bandwidths(A)
     if -l > u # no bands
-        _fill_lmul!(β, y)
+        _fill_rmul!(y, β)
     elseif l < 0 # with u >= -l > 0, that is, all bands lie above the diagonal
         # E.g. (l,u) = (-1,2)
         # set lview = 0 and uview = u + l >= 0
         _banded_gbmv!('N', α, view(A, :, 1-l:n), view(x, 1-l:n), β, y)
     elseif u < 0 # with -l <= u < 0, that is, all bands lie below the diagnoal.
         # E.g. (l,u) = (2,-1)
         # set lview = l + u >= 0 and uview = 0
-        _fill_lmul!(β, @view(y[1:-u]))
+        _fill_rmul!(@view(y[1:-u]), β)
         _banded_gbmv!('N', α, view(A, 1-u:m, :), x, β, view(y, 1-u:m))
         y
     else
@@ -67,11 +67,11 @@ function _banded_muladd_row!(tA, α, At, x, β, y)
     n, m = size(At)
     u, l = bandwidths(At)
     if -l > u # no bands
-        _fill_lmul!(β, y)
+        _fill_rmul!(y, β)
     elseif l < 0
         _banded_gbmv!(tA, α, view(At, 1-l:n, :,), view(x, 1-l:n), β, y)
     elseif u < 0
-        _fill_lmul!(β, @view(y[1:-u]))
+        _fill_rmul!(@view(y[1:-u]), β)
         _banded_gbmv!(tA, α, view(At, :, 1-u:m), x, β, view(y, 1-u:m))
         y
     else
@@ -100,10 +100,10 @@ end
 @inline function materialize!(M::MatMulVecAdd{<:AbstractBandedLayout})
     checkdimensions(M)
     α,A,B,β,C = M.α,M.A,M.B,M.β,M.C
-    _fill_lmul!(β, C)
+    _fill_rmul!(C, β)
     @inbounds for j = intersect(rowsupport(A), colsupport(B))
         for k = colrange(A,j)
-            C[k] += α*inbands_getindex(A,k,j)*B[j]
+            C[k] += inbands_getindex(A,k,j) * B[j] * α
         end
     end
     C
@@ -113,11 +113,11 @@ end
     checkdimensions(M)
     α,At,B,β,C = M.α,M.A,M.B,M.β,M.C
     A = transpose(At)
-    _fill_lmul!(β, C)
+    _fill_rmul!(C, β)
 
     @inbounds for j = rowsupport(A)
         for k = intersect(colrange(A,j), colsupport(B))
-            C[j] +=  α*transpose(inbands_getindex(A,k,j))*B[k]
+            C[j] +=  transpose(inbands_getindex(A,k,j)) * B[k] * α
         end
     end
     C
@@ -127,10 +127,10 @@ end
     checkdimensions(M)
     α,Ac,B,β,C = M.α,M.A,M.B,M.β,M.C
     A = Ac'
-    _fill_lmul!(β, C)
+    _fill_rmul!(C, β)
     @inbounds for j = rowsupport(A)
         for k = intersect(colrange(A,j), colsupport(B))
-            C[j] += α*inbands_getindex(A,k,j)'*B[k]
+            C[j] += inbands_getindex(A,k,j)' * B[k] * α
         end
     end
     C
@@ -227,13 +227,13 @@ end
 
 function materialize!(M::MatMulMatAdd{<:DiagonalLayout{<:AbstractFillLayout},<:AbstractBandedLayout})
     checkdimensions(M)
-    M.C .= (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C
+    M.C .= getindex_value(M.A.diag) .* M.B .* M.α .+ M.C .* M.β
     M.C
 end
 
 function materialize!(M::MatMulMatAdd{<:AbstractBandedLayout,<:DiagonalLayout{<:AbstractFillLayout}})
     checkdimensions(M)
-    M.C .= (M.α * getindex_value(M.B.diag)) .* M.A .+ M.β .* M.C
+    M.C .= M.A .* getindex_value(M.B.diag) .* M.α .+ M.C .* M.β
     M.C
 end
 
@@ -244,7 +244,7 @@ function materialize!(M::MatMulMatAdd{<:BandedColumns, <:AbstractStridedLayout,
     α, β, A, B, C = M.α, M.β, M.A, M.B, M.C
 
     if iszero(α)
-        lmul!(β, C)
+        rmul!(C, β)
     else
         for (colC, colB) in zip(eachcol(C), eachcol(B))
             mul!(colC, A, colB, α, β)
@@ -259,7 +259,7 @@ function materialize!(M::MatMulMatAdd{<:AbstractStridedLayout, <:BandedColumns,
     α, β, A, B, C = M.α, M.β, M.A, M.B, M.C
 
     if iszero(α)
-        lmul!(β, C)
+        rmul!(C, β)
     else
         for (rowC, rowA) in zip(eachrow(C), eachrow(A))
             mul!(rowC, transpose(B), rowA, α, β)

src/generic/utils.jl

@@ -24,3 +24,7 @@ prodbandwidths(A...) = broadcast(+, bandwidths.(A)...)
 function sumbandwidths(A::AbstractMatrix, B::AbstractMatrix)
     max(bandwidth(A, 1), bandwidth(B, 1)), max(bandwidth(A, 2), bandwidth(B, 2))
 end
+
+
+_fill_lmul!(β, A::AbstractArray{T}) where T = iszero(β) ? zero!(A) : lmul!(β, A)
+_fill_rmul!(A::AbstractArray{T}, β) where T = iszero(β) ? zero!(A) : rmul!(A, β)

test/test_linalg.jl

@@ -1,4 +1,9 @@
-using ArrayLayouts, BandedMatrices, FillArrays, LinearAlgebra, Test
+using ArrayLayouts
+using BandedMatrices
+using FillArrays
+using LinearAlgebra
+using Quaternions
+using Test
 
 import Base.Broadcast: materialize, broadcasted
 import BandedMatrices: BandedColumns, _BandedMatrix
@@ -109,6 +114,16 @@ ArrayLayouts.colsupport(::UnknownLayout, A::MyOneElement{<:Any,1}, _) =
             @test B*M ≈ Matrix(B)*M
             @test M*B ≈ M*Matrix(B)
         end
+        @testset "non-commutative" begin
+            B1 = BandedMatrix(0 => [quat(rand(4)...) for i in 1:3])
+            v = [quat(rand(4)...) for i in 1:3]
+            α, β = quat(0,1,1,0), quat(1,0,0,1)
+            @test mul!(zero(v), B1, v, α, β) ≈ mul!(zero(v), Array(B1), v, α, β)
+
+            D = Diagonal(Fill(quat(rand(4)...), size(B1,2)))
+            @test mul!(D * B1, D, B1, α, β) ≈ mul!(D * Array(B1), D, Array(B1), α, β)
+            @test mul!(B1 * D, B1, D, α, β) ≈ mul!(Array(B1) * D, Array(B1), D, α, β)
+        end
     end
 
     @testset "BandedMatrix * sparse" begin