diff --git a/Project.toml b/Project.toml
index a29a3f75..b7dbe159 100644
--- a/Project.toml
+++ b/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"]
diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl
index 5ce905d2..9beaad7b 100644
--- a/src/BandedMatrices.jl
+++ b/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
diff --git a/src/generic/matmul.jl b/src/generic/matmul.jl
index 93581c75..34704500 100644
--- a/src/generic/matmul.jl
+++ b/src/generic/matmul.jl
@@ -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 .= getindex_value(M.B.diag) .* M.A .* M.α .+ M.C .* M.β
     M.C
 end
 
diff --git a/src/generic/utils.jl b/src/generic/utils.jl
index a329c270..5f6978aa 100644
--- a/src/generic/utils.jl
+++ b/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, β)
diff --git a/test/test_linalg.jl b/test/test_linalg.jl
index fd102f04..8b398620 100644
--- a/test/test_linalg.jl
+++ b/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,12 @@ 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, α, β)
+        end
     end
 
     @testset "BandedMatrix * sparse" begin