Fix generic axpy! for vector element types without commutative multiplication,

e.g. quaternions.

Author: andreasnoack

base/linalg/generic.jl

@@ -439,20 +439,20 @@ function peakflops(n::Integer=2000; parallel::Bool=false)
     parallel ? sum(pmap(peakflops, [ n for i in 1:nworkers()])) : (2*Float64(n)^3/t)
 end
 
-# BLAS-like in-place y=alpha*x+y function (see also the version in blas.jl
+# BLAS-like in-place y = x*α+y function (see also the version in blas.jl
 #                                          for BlasFloat Arrays)
-function axpy!(alpha, x::AbstractArray, y::AbstractArray)
+function axpy!(α, x::AbstractArray, y::AbstractArray)
     n = length(x)
     if n != length(y)
         throw(DimensionMismatch("x has length $n, but y has length $(length(y))"))
     end
     for i = 1:n
-        @inbounds y[i] += alpha * x[i]
+        @inbounds y[i] += x[i]*α
     end
     y
 end
 
-function axpy!{Ti<:Integer,Tj<:Integer}(alpha, x::AbstractArray, rx::AbstractArray{Ti}, y::AbstractArray, ry::AbstractArray{Tj})
+function axpy!{Ti<:Integer,Tj<:Integer}(α, x::AbstractArray, rx::AbstractArray{Ti}, y::AbstractArray, ry::AbstractArray{Tj})
     if length(x) != length(y)
         throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))"))
     elseif minimum(rx) < 1 || maximum(rx) > length(x)
@@ -463,7 +463,7 @@ function axpy!{Ti<:Integer,Tj<:Integer}(alpha, x::AbstractArray, rx::AbstractArr
         throw(ArgumentError("rx has length $(length(rx)), but ry has length $(length(ry))"))
     end
     for i = 1:length(rx)
-        @inbounds y[ry[i]] += alpha * x[rx[i]]
+        @inbounds y[ry[i]] += x[rx[i]]*α
     end
     y
 end

test/linalg/generic.jl

@@ -141,3 +141,12 @@ let x = Vector{Int}[[1,2], [3,4]]
     @test norm(x, 1) ≈ sqrt(5) + 5
     @test norm(x, 3) ≈ cbrt(sqrt(125)+125)
 end
+
+# test that LinAlg.axpy! works for element type without commutative multiplication
+let
+    α = ones(Int, 2, 2)
+    x = fill([1 0; 1 1], 3)
+    y = fill(zeros(Int, 2, 2), 3)
+    @test LinAlg.axpy!(α, x, deepcopy(y)) == x .* Matrix{Int}[α]
+    @test LinAlg.axpy!(α, x, deepcopy(y)) != Matrix{Int}[α] .* x
+end