Closed
Description
Using the code
using ForwardDiff: Dual
using StaticArrays
v = ((1.,2.), (2.,3.))
function f(x::Vector{T}) where {T}
x1 = SVector((x[1],)) # comment this to get correct result
la1 = SVector((x[1],))
f1 = zero(SVector{1}) # can also initialize with zero(SVector{1, T}) to get correct result
for _ in v
f1 += la1
@show f1
end
return f1
end
x = [Dual(1.0,1.0)]
@show f(x)
gives the erroneous output:
f1 = Dual{Nothing,Float64,1}[Dual{Nothing}(1.0,1.0)]
f1 = Dual{Nothing,Float64,1}[Dual{Nothing}(1.0,1.0)]
f(x) = Dual{Nothing,Float64,1}[Dual{Nothing}(1.0,1.0)]
Now, we comment out the line x1 = SVector((x[1],))
f1 = Dual{Nothing,Float64,1}[Dual{Nothing}(1.0,1.0)]
f1 = Dual{Nothing,Float64,1}[Dual{Nothing}(2.0,2.0)]
f(x) = Dual{Nothing,Float64,1}[Dual{Nothing}(2.0,2.0)]
and get the correct result.
Executing the original code in the REPL also gives the correct result
julia> let
x1 = SVector((x[1],))
la1 = SVector((x[1],))
f1 = zero(SVector{1})
for _ in v
f1 += la1
@show f1
end
f1
end
f1 = Dual{Nothing,Float64,1}[Dual{Nothing}(1.0,1.0)]
f1 = Dual{Nothing,Float64,1}[Dual{Nothing}(2.0,2.0)]
1-element SArray{Tuple{1},Dual{Nothing,Float64,1},1,1}:
Dual{Nothing}(2.0,2.0)
If we initialize f1
to be of type T
we also get correct result.
My guess it is something with the type instability in combination with broadcasting in combination with optimization but not sure...