fix alpha != 1

Project.toml

@@ -1,7 +1,7 @@
 name = "Octavian"
 uuid = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4"
 authors = ["Mason Protter", "Chris Elrod", "Dilum Aluthge", "contributors"]
-version = "0.3.6"
+version = "0.3.7"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

src/forward_diff.jl

@@ -68,8 +68,12 @@ _view1(B::AbstractArray{<:Any,3}) = @view(B[1,:,:])
     end
   end
   Pstatic = static(P)
-  @tturbo for n ∈ indices((B,C),3), m ∈ indices((A,C),2), p ∈ 1:Pstatic, k ∈ indices((A,B),(3,2))
-    C[p+1,m,n] += A[1,m,k] * B[p+1,k,n]
+  @tturbo for n ∈ indices((B,C),3), m ∈ indices((A,C),2), p ∈ 1:Pstatic
+    Cₚₘₙ = zero(eltype(C))
+    for k ∈ indices((A,B),(3,2))
+      Cₚₘₙ += A[1,m,k] * B[p+1,k,n]
+    end
+    C[p+1,m,n] = C[p+1,m,n] + α*Cₚₘₙ
   end
   _C
 end

src/matmul.jl

@@ -220,11 +220,13 @@ end
 
 Multiply matrices `A` and `B`.
 """
-@inline function matmul(A::AbstractMatrix, B::AbstractVecOrMat)
+@inline function matmul(A::AbstractMatrix, B::AbstractVecOrMat, α, β)
   C, (M,K,N) = alloc_matmul_product(A, B)
-  matmul!(C, A, B, One(), Zero(), nothing, (M,K,N), ArrayInterface.contiguous_axis(C))
+  matmul!(C, A, B, α, β, nothing, (M,K,N), ArrayInterface.contiguous_axis(C))
   return C
 end
+@inline matmul(A::AbstractMatrix, B::AbstractVecOrMat) = matmul(A, B, One(), Zero())
+@inline matmul(A::AbstractMatrix, B::AbstractVecOrMat, α) = matmul(A, B, α, Zero())
 
 """
     matmul!(C, A, B[, α, β, max_threads])

test/forward_diff.jl

@@ -65,7 +65,7 @@ randdual(x, ::Val{N}=Val(3)) where {N} = ForwardDiff.Dual(x, ntuple(_ -> randn()
   @testset "two dual arrays" begin
     A1d = randdual.(A1)
     B1d = randdual.(B1)
-    @test reinterpret(Float64, Octavian.matmul(A1d, B1d)) ≈ reinterpret(Float64, A1d * B1d)
+    @test reinterpret(Float64, Octavian.matmul(A1d, B1d, 1.3)) ≈ reinterpret(Float64, (A1d * B1d) .* 1.3)
     @test reinterpret(Float64, Octavian.matmul(@view(A1d[begin:end-1,:]), B1d)) ≈ reinterpret(Float64, @view(A1d[begin:end-1,:]) * B1d)
   end
 end