ForwardDiff on LazySum is failing

[AmitRotem]

using QuantumOptics
import ForwardDiff as FD
mb = NLevelBasis(5)
s = randoperator(mb)
h = [randoperator(mb),randoperator(mb)]

LZtr(p) = abs2(tr(s'*(LazySum(p, h)*s)))

LZtr_workaround(p) = abs2(tr(s'*(LazySum(p, h)*s)))
function LZtr_workaround(p::Vector{<:FD.Dual})
    LS = LazySum(p, h)
    S = sum(LS.factors.*LS.operators)
    abs2(tr(s'*(S*s)))
end

FD.gradient(LZtr_workaround, rand(2)) |> display # ok
FD.gradient(LZtr, rand(2)) # error

seems to come from mul!(R,M,A,α,β) tring to convert (M*A)*α, which is type a has FD.Dual elements, to R which have ComplexF64 elements.