Workaround type instability in complex _residual!

src/codim2/diffeq.jl

@@ -248,7 +248,18 @@ function _residual!(::Any, u, prob::DiffEqCodim2Problem,
 
     v = ds_eigvec(prob, u)
     iw = ds_eigval(prob, u)
+
+    # Whe v is complex and _residual! is invoked via ForwardDiff,
+    # eltype(J * v) becomes Any.  `as_reals` relies on the eltype to
+    # do the right thing so I need to fix this.  Here is a hack to
+    # workaround it.
+    H2ʹ = fixeltype(promote_type(eltype.((v, J, iw))...),  # TODO: don't
+                    J * v .- iw .* v)
+    H2 = as_reals(H2ʹ)
+    #=
     H2 = as_reals(J * v .- iw .* v)
+    =#
+
     H3 = eigvec_constraint(v, augsys_cache)
 
     return cat_outputs(

src/codim2/utils.jl

@@ -90,3 +90,25 @@ as_reals(v::Vector{<: Complex{T}}) where T = reinterpret(T, v)
 
 as_reals(s::Real) = SVector(s)
 as_reals(s::Complex) = SVector(real(s), imag(s))
+
+fixeltype(::Type{T}, v::AbstractArray{T}) where {T} = v
+
+@generated function fixeltype(::Type{E}, v::SVector{S, Any}) where {E, S}
+    values = [:(v[$i]) for i in 1:S]
+    quote
+        SVector{$S, $E}($(values...))
+    end
+end
+
+#=
+@generated function fixeltype(::Type{E}, M::SMatrix) where {E}
+    s = Size(M)
+    S1 = s[1]
+    S2 = s[2]
+    L = S1 * S2
+    values = [:(M[$i]) for i in 1:L]
+    quote
+        SMatrix{$S1, $S2, $E}($(values...))
+    end
+end
+=#