From cad8082312a695263daddbf93bee015bf68d1be5 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Sun, 18 Nov 2018 11:28:42 -0800
Subject: [PATCH] Test multiply-add interface

---
 stdlib/LinearAlgebra/test/matmul.jl | 46 ++++++++++++++++++++++++++++-
 1 file changed, 45 insertions(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl
index 62c7d7cc8faf2..462d1269ced51 100644
--- a/stdlib/LinearAlgebra/test/matmul.jl
+++ b/stdlib/LinearAlgebra/test/matmul.jl
@@ -3,7 +3,7 @@
 module TestMatmul
 
 using Test, LinearAlgebra, Random
-using LinearAlgebra: mul!
+using LinearAlgebra: mul!, addmul!
 
 ## Test Julia fallbacks to BLAS routines
 
@@ -105,6 +105,19 @@ end
         @test mul!(C, transpose(A), transpose(B)) == A'*B'
         @test LinearAlgebra.mul!(C, adjoint(A), transpose(B)) == A'*transpose(B)
 
+        # Inplace multiply-add
+        α = rand(-10:10)
+        β = rand(-10:10)
+        rand!(C, -10:10)
+        βC = β * C
+        _C0 = copy(C)
+        C0() = (C .= _C0; C)  # reset C but don't change the container type
+        @test addmul!(C0(), A, B, α, β) == α*A*B .+ βC
+        @test addmul!(C0(), transpose(A), B, α, β) == α*A'*B .+ βC
+        @test addmul!(C0(), A, transpose(B), α, β) == α*A*B' .+ βC
+        @test addmul!(C0(), transpose(A), transpose(B), α, β) == α*A'*B' .+ βC
+        @test addmul!(C0(), adjoint(A), transpose(B), α, β) == α*A'*transpose(B) .+ βC
+
         #test DimensionMismatch for generic_matmatmul
         @test_throws DimensionMismatch LinearAlgebra.mul!(C, adjoint(A), transpose(fill(1,4,4)))
         @test_throws DimensionMismatch LinearAlgebra.mul!(C, adjoint(fill(1,4,4)), transpose(B))
@@ -113,6 +126,9 @@ end
     CC = Matrix{Int}(undef, 2, 2)
     for v in (copy(vv), view(vv, 1:2)), C in (copy(CC), view(CC, 1:2, 1:2))
         @test @inferred(mul!(C, v, adjoint(v))) == [1 2; 2 4]
+
+        C .= [1 0; 0 1]
+        @test @inferred(addmul!(C, v, adjoint(v), 2, 3)) == [5 4; 4 11]
     end
 end
 
@@ -127,11 +143,15 @@ end
     CC = Matrix{Int}(undef, 3, 3)
     for v in (copy(vv), view(vv, 1:3)), C in (copy(CC), view(CC, 1:3, 1:3))
         @test mul!(C, v, transpose(v)) == v*v'
+        C .= C0 = rand(-10:10, size(C))
+        @test addmul!(C, v, transpose(v), 2, 3) == 2v*v' .+ 3C0
     end
     vvf = map(Float64,vv)
     CC = Matrix{Float64}(undef, 3, 3)
     for vf in (copy(vvf), view(vvf, 1:3)), C in (copy(CC), view(CC, 1:3, 1:3))
         @test mul!(C, vf, transpose(vf)) == vf*vf'
+        C .= C0 = rand(eltype(C), size(C))
+        @test addmul!(C, vf, transpose(vf), 2, 3) == 2vf*vf' .+ 3C0
     end
 end
 
@@ -143,6 +163,17 @@ end
         @test LinearAlgebra.mul!(C, transpose(A), transpose(B)) == transpose(A)*transpose(B)
         @test LinearAlgebra.mul!(C, A, adjoint(B)) == A*transpose(B)
         @test LinearAlgebra.mul!(C, adjoint(A), B) == transpose(A)*B
+
+        # Inplace multiply-add
+        α = rand(Float64)
+        β = rand(Float64)
+        rand!(C)
+        βC = β * C
+        _C0 = copy(C)
+        C0() = (C .= _C0; C)  # reset C but don't change the container type
+        @test addmul!(C0(), transpose(A), transpose(B), α, β) ≈ α*transpose(A)*transpose(B) .+ βC
+        @test addmul!(C0(), A, adjoint(B), α, β) ≈ α*A*transpose(B) .+ βC
+        @test addmul!(C0(), adjoint(A), B, α, β) ≈ α*transpose(A)*B .+ βC
     end
 end
 
@@ -350,6 +381,9 @@ Transpose(x::RootInt) = x
     C = [0]
     mul!(C, a, transpose(a))
     @test C[1] == 9
+    C = [1]
+    addmul!(C, a, transpose(a), 2, 3)
+    @test C[1] == 21
     a = [RootInt(2),RootInt(10)]
     @test a*adjoint(a) == [4 20; 20 100]
     A = [RootInt(3) RootInt(5)]
@@ -360,6 +394,16 @@ function test_mul(C, A, B)
     mul!(C, A, B)
     @test Array(A) * Array(B) ≈ C
     @test A*B ≈ C
+
+    rand!(C)
+    T = promote_type(eltype.((A, B))...)
+    α = rand(T)
+    β = rand(T)
+    βArrayC = β * Array(C)
+    βC = β * C
+    addmul!(C, A, B, α, β)
+    @test α * Array(A) * Array(B) .+ βArrayC ≈ C
+    @test α * A * B .+ βC ≈ C
 end
 
 @testset "mul! vs * for special types" begin