Inplace multiplicacion for Matrix{Taylor1} and Matrix{TaylorN} (#403)

* Add taylormatmul! for allocation-free matrix multiplication with Taylor1 or TaylorN elements

Co-authored-by: LuEdRaMo <73906617+LuEdRaMo@users.noreply.github.com>

* Add some tests

* Rename to matmul!

* Bump patch version

* Reduce tolerance of tests, and add new

---------

Co-authored-by: LuEdRaMo <73906617+LuEdRaMo@users.noreply.github.com>

Author: lbenet

Project.toml

@@ -1,6 +1,6 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
-version = "0.21.3"
+version = "0.21.4"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
 
 [deps]

src/arithmetic.jl

@@ -1461,3 +1461,24 @@ function lu(A::AbstractMatrix{Taylor1{T}}; check::Bool = true) where {T<:Number}
         return lu!(copy_oftype(A, S), LU_RowMaximum; check = check)
     end
 end
+
+
+# Fast allocation-free matrix multiplication
+for T in (:Taylor1, :TaylorN)
+    @eval function matmul!(C::Matrix{$T{T}},
+                           A::Matrix{$T{T}}, B::Matrix{$T{T}}) where {T}
+        mc, nc = size(C)
+        ma, na = size(A)
+        mb, nb = size(B)
+        @assert (na == mb && mc == ma && nc == nb)
+        for j in axes(C,2)
+            for i in axes(C,1)
+                TS.zero!(C[i,j])
+                for k in 1:na
+                    TS.muladd!(C[i,j], A[i,k], B[k,j])
+                end
+            end
+        end
+        return nothing
+    end
+end

test/onevariable.jl

@@ -691,6 +691,21 @@ end
     end
 end
 
+@testset "Matrix multiplication with Taylor1" begin
+    A = [Taylor1.(randn.(16)) for i=1:3, j=1:3]
+    B = [Taylor1.(randn.(16)) for i=1:3, j=1:3]
+    C = [Taylor1.(randn.(16)) for i=1:3, j=1:3]
+    D = [Taylor1.(randn.(16)) for i=1:3, j=1:3]
+    for _ = 1:4
+        A .= [Taylor1.(randn.(16)) for i=1:3, j=1:3]
+        B .= [Taylor1.(randn.(16)) for i=1:3, j=1:3]
+        TS.matmul!(C, A, B)
+        mul!(D, A, B) # usual matrix product
+        @test all(D .≈ C)
+        @test norm(D - C, Inf) < 1.0e-13
+    end
+end
+
 # @testset "Matrix multiplication for Taylor1" begin
 #     order = 30
 #     n1 = 100