Conformity with LinearAlgebra's muladd order

[jishnub]

This package computes α * A * B + β * C whereas LinearAlgebra computes A * B * α + C * β. These are equivalent for real/complex numbers (e.g. BLAS calls), but these would differ for non-commutative types. Perhaps we should switch to the same convention?

Currently,

julia> using Quaternions, LinearAlgebra, BandedMatrices

julia> B1 = BandedMatrix(0 => [quat(rand(4)...) for i in 1:3]);

julia> B2 = BandedMatrix(0 => [quat(rand(4)...) for i in 1:3]);

julia> α, β = quat(0,1,1,0), quat(1,0,0,1);

julia> mul!(copy(B1), B1, B2, α, β) ≈ mul!(Array(B1), Array(B1), Array(B2), α, β)
false

[dlfivefifty]

I have no concerns with making this change but I'm a bit confused by why LinearAlgebra uses that order...