faster * and / of dual arr by constant mtx

Author: alyst

src/dual.jl

@@ -773,6 +773,60 @@ function SpecialFunctions.gamma_inc(a::Real, d::Dual{T,<:Real}, ind::Integer) wh
     return (Dual{T}(p, ∂p), Dual{T}(q, -∂p))
 end
 
+# Efficient left multiplication/division of #
+# Dual array by a constant matrix           #
+#-------------------------------------------#
+
+# creates the copy of x and applies fvalue!(values) to its values, and fpartial!(partial, ix) to its partials
+function _map_dual_components(fvalue!, fpartial!, x::AbstractArray{DT}) where DT <: Dual{<:Any, T} where T
+    N = npartials(DT)
+    @assert Base.elsize(x) == sizeof(T) * (N + 1) # check that accessing x as Vector{T} is possible
+    res = similar(x) # result
+    t = similar(x, T) # temporary Array{T} for fvalue!/fpartial! application
+    ystride = N + 1
+    # y allows res to be accessed as Vector{T}
+    y = Base.unsafe_wrap(Array, Base.unsafe_convert(Ptr{T}, res),
+                         length(res) * ystride, own=false)
+
+    # calculate res values
+    @inbounds for (j, v) in enumerate(x)
+        t[j] = value(v)
+    end
+    fvalue!(t)
+    k = 1
+    @inbounds for tt in t
+        y[k] = tt
+        k += ystride
+    end
+
+    # calculate each res partial
+    for i in 1:N
+        @inbounds for (j, v) in enumerate(x)
+            t[j] = partials(v, i)
+        end
+        fpartial!(t, i)
+        k = i + 1
+        @inbounds for tt in t
+            y[k] = tt
+            k += ystride
+        end
+    end
+
+    return res
+end
+
+Base.:\(m::Union{LowerTriangular{<:LinearAlgebra.BlasFloat},
+                 UpperTriangular{<:LinearAlgebra.BlasFloat},
+                 Matrix{<:LinearAlgebra.BlasFloat}},
+        x::AbstractVecOrMat{<:Dual}) =
+    _map_dual_components(Base.Fix1(ldiv!, m), (x, _) -> ldiv!(m, x), x)
+
+Base.:*(m::Union{LowerTriangular{<:LinearAlgebra.BlasFloat},
+                 UpperTriangular{<:LinearAlgebra.BlasFloat},
+                 Matrix{<:LinearAlgebra.BlasFloat}},
+        x::AbstractVecOrMat{<:Dual}) =
+    _map_dual_components(Base.Fix1(lmul!, m), (x, _) -> ldiv!(m, x), x)
+
 ###################
 # Pretty Printing #
 ###################