# My function works on it's own but not in NonLinearSolve

**URL:** <https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204>\
**Category:** General Usage\
**Tags:** forwarddiff, nonlinearsolve\
**Created:** [April 18, 2024, 5:29pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204 "2024-04-18T17:29:31Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Laura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laura/32/207742_2.png) [@Laura](https://discourse.julialang.org/u/Laura)\
**Post date:** [April 18, 2024, 5:29pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/1 "2024-04-18T17:29:32Z")

</div>

I have a set of equations that I would like to solve using the non linear equation solver NonLinearSolve.jl. I have writen a function for this which runs fine on it’s own but when I try to solve it I get an error message.

The equation I would like to solve is the following:

```julia
cd("C:/Users/...correct directory etc")

using NonlinearSolve
using LinearAlgebra
using JLD2

mutable struct parameters

    n::Int64
    matrix::Matrix{Float64}
    tensor::Array{Float64}
    norm::Vector{Float64}

end

function equation!(y, z, p::parameters)

    @inbounds for i = 1:p.n
        for j = 1:i
            acc = zero(eltype(p.matrix))
            for k = 1:p.n
                acctmp = zero(eltype(p.matrix))
                @simd for l = 1:p.n
                    acctmp += p.tensor[l, k, j, i] * z[p.n+l]
                end
                acc += acctmp * z[k]
            end
            p.matrix[i, p.n+j] = p.matrix[j, p.n+i] = p.matrix[p.n+i, j] = p.matrix[p.n+j, i] = acc
        end
    end

    p.matrix[p.n+1:2p.n, 1:p.n] .*= p.norm / (z[1] * z[p.n+1])
    p.matrix[1:p.n, p.n+1:2p.n] .*= p.norm / (z[1] * z[p.n+1])

    y[1:2p.n] = (p.matrix - I * z[2p.n+1]) * z[1:2p.n]

end

load_object(string("example u.jld2"))
load_object(string("example p.jld2"))

prob = NonlinearProblem(equation!, u, p)

sol = solve(prob)

```

Now, I tried to upload the example u and p files to this post but that does not seem to be possible. I hope perhaps someone can just see the problem without needing those files ore perhaps knows how to upload them. 😅

When I run this I get the following error message:

```julia
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
   @ Base char.jl:50
  ...

```

With the following stacktrace:

```julia
Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11})
    @ Base .\number.jl:7
  [2] setindex!
    @ .\array.jl:971 [inlined]
  [3] equation!(y::Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, z::Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, p::parameters)
    @ Main c:\Users\laura\Documents\Universiteit\Doctoraat\Impurities in Fermi Fluids\Julia\example code.jl:29
  [4] NonlinearFunction
    @ C:\Users\laura\.julia\packages\SciMLBase\2HZ5m\src\scimlfunctions.jl:2356 [inlined]
  [5] JacobianWrapper
    @ C:\Users\laura\.julia\packages\SciMLBase\2HZ5m\src\function_wrappers.jl:99 [inlined]
  [6] chunk_mode_jacobian!(result::Matrix{Float64}, f!::SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, parameters}, y::Vector{Float64}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}}})
    @ ForwardDiff C:\Users\laura\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:183
  [7] jacobian!
    @ C:\Users\laura\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:80 [inlined]
  [8] jacobian!
    @ C:\Users\laura\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:76 [inlined]
  [9] sparse_jacobian!
    @ C:\Users\laura\.julia\packages\SparseDiffTools\CsonI\src\highlevel\forward_mode.jl:70 [inlined]
 [10] #jacobian!!#130
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\jacobian.jl:62 [inlined]
 [11] jacobian!!
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\jacobian.jl:53 [inlined]
 [12] perform_step!(cache::NonlinearSolve.BroydenCache{true, :true_jacobian, :good_broyden, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Broyden{:true_jacobian, :good_broyden, true, AutoForwardDiff{nothing, ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}}, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, parameters, SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, parameters}, Matrix{Float64}, Nothing, Vector{Float64}, Float64, Nothing, typeof(DiffEqBase.NONLINEARSOLVE_DEFAULT_NORM), Float64, Float64, Float64, NonlinearSolve.var"#118#119"{Float64}, SparseDiffTools.ForwardDiffJacobianCache{SparseDiffTools.NoMatrixColoring, ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}}}, Nothing, Vector{Float64}, Vector{Float64}}, NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.NoLineSearchCache{Float64}, DiffEqBase.NonlinearTerminationModeCache{Vector{Float64}, Float64, AbsSafeBestTerminationMode{Nothing, Int64, Float64}, Float64, Vector{Float64}, Nothing}, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\broyden.jl:171
 [13] solve!(cache::NonlinearSolve.BroydenCache{true, :true_jacobian, :good_broyden, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Broyden{:true_jacobian, :good_broyden, true, AutoForwardDiff{nothing, ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}}, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, parameters, SciMLBase.JacobianWrapper{true, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, parameters}, Matrix{Float64}, Nothing, Vector{Float64}, Float64, Nothing, typeof(DiffEqBase.NONLINEARSOLVE_DEFAULT_NORM), Float64, Float64, Float64, NonlinearSolve.var"#118#119"{Float64}, SparseDiffTools.ForwardDiffJacobianCache{SparseDiffTools.NoMatrixColoring, ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.NonlinearSolveTag, Float64}, Float64, 11}}}}, Nothing, Vector{Float64}, Vector{Float64}}, NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.NoLineSearchCache{Float64}, DiffEqBase.NonlinearTerminationModeCache{Vector{Float64}, Float64, AbsSafeBestTerminationMode{Nothing, Int64, Float64}, Float64, Vector{Float64}, Nothing}, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\NonlinearSolve.jl:149
 [14] __solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, ::Broyden{:true_jacobian, :good_broyden, true, Nothing, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:default_set, :second_time), Tuple{Bool, Bool}}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\NonlinearSolve.jl:135
 [15] __solve
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\NonlinearSolve.jl:132 [inlined]
 [16] macro expansion
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:123 [inlined]
 [17] __solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, ::NonlinearSolvePolyAlgorithm{:NLS, 7, Tuple{Broyden{:identity, :good_broyden, false, Nothing, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Broyden{:true_jacobian, :good_broyden, true, Nothing, Nothing, LineSearch{Nothing, Nothing, Bool}, Nothing}, Klement{:identity, false, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Bool}, NewtonRaphson{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}}, NewtonRaphson{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{BackTracking{Float64, Int64}, Nothing, Bool}}, TrustRegion{nothing, Nothing, Rational{Int64}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), Rational{Int64}, Nothing}, TrustRegion{nothing, Nothing, Rational{Int64}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), Rational{Int64}, Nothing}}}; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:default_set, :second_time), Tuple{Bool, Bool}}})
    @ NonlinearSolve C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:115
 [18] __solve
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:115 [inlined]
 [19] #__solve#196
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:386 [inlined]
 [20] __solve
    @ C:\Users\laura\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:383 [inlined]
 [21] #__solve#61
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1370 [inlined]
 [22] __solve
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1363 [inlined]
 [23] #solve_call#34
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:608 [inlined]
 [24] solve_call
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:566 [inlined]
 [25] #solve_up#42
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1049 [inlined]
 [26] solve_up
    @ C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1043 [inlined]
 [27] solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1037
 [28] solve(::NonlinearProblem{Vector{Float64}, true, parameters, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(equation!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem})
    @ DiffEqBase C:\Users\laura\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1027
 [29] top-level scope
    @ c:\Users\laura\Documents\Universiteit\Doctoraat\Impurities in Fermi Fluids\Julia\example code.jl:45

```

In my limited understanding of what is happening here, the non linear solver tries to calculate the jacobian of my function by propagating some object of the type _::ForewardDiff_ through my function which does not work for some reason.

Does anyone know if there is a way to change my function in such a way that it is possible for the solver to propagate this _::ForewardDiff_ object through it or should I attempt to make a seperate jacobian function?

Or perhaps the problem is something else all together?

Best regards,  
Laura

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 18, 2024, 5:32pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/2 "2024-04-18T17:32:02Z")

</div>

It’s the same issue as:

> **[Frequently Asked Questions · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers)**
>
> Documentation for DifferentialEquations.jl.

You’ll want to make use of PreallocationTools.jl if you’re precaching values that are then being autodiff’d:

> **[GitHub - SciML/PreallocationTools.jl: Tools for building non-allocating...](https://github.com/SciML/PreallocationTools.jl)**
>
> Tools for building non-allocating pre-cached functions in Julia, allowing for GC-free usage of automatic differentiation in complex codes - SciML/PreallocationTools.jl

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [April 18, 2024, 6:06pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/3 "2024-04-18T18:06:24Z")

</div>

> [@Laura](#):
>
> ```julia
> mutable struct parameters
> n::Int64
> matrix::Matrix{Float64}
> tensor::Array{Float64}
> norm::Vector{Float64}
> end
> 
> ```

In this case your code seems fully generic with respect to element types, the only thing you need to do is allow the same genericity in your struct fields:

```julia
mutable struct parameters{T1,T2,T3}
    n::Int64
    matrix::Matrix{T1}
    tensor::Array{T2}
    norm::Vector{T3}
end

```

You probably also want to specify the order of the tensor to make the field type-inferrable (at the moment Julia doesn’t know how many dimensions the `Array` has).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 18, 2024, 6:56pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/4 "2024-04-18T18:56:54Z")

</div>

> [@gdalle](#):
>
> In this case your code seems fully generic with respect to element types, the only thing you need to do is allow the same genericity in your struct fields:
> 
> ```julia
> 
> ```

That doesn’t solve it because the `p` is constructed outside of the function, and so it needs to be generic to the multiple cache sizes. `matrix`, `tensor`, and `norm` need to be PreallocationTools.DualCache objects.

---

<div class="post-metadata">

**Author:** ![Laura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laura/32/207742_2.png) [@Laura](https://discourse.julialang.org/u/Laura)\
**Post date:** [April 19, 2024, 11:21am UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/5 "2024-04-19T11:21:54Z")

</div>

Making the the element types generic and specifying the array dimensions does indeed not seem to solve this. I’ll have a look at the Autodifferentiation and Dual Numbers.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [April 19, 2024, 1:26pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/6 "2024-04-19T13:26:37Z")

</div>

The concrete problem is the `acc` and `acctmp` variables, which you set to `zero(eltype(p.matrix))`, i.e. a zero of type `Float64`. In the loop you create the `acctmp` by multiplying with `z[...]`. Now, when differentiating, both `y` and `z` contains `Number`s of type `ForwardDiff.Dual` (which can be seen near the top of the stack dump). Now, a `Dual` can be multiplied by `Float64` to yield a new `Dual`. This makes `acctmp` a `Dual`, and then `acc` becomes a `Dual`. Then you try to put it into `p.matrix[...]`, which contains `Float64`. This fails, that’s the error message, because a `Dual` can’t be converted to a `Float64`.

To work with autodiff, your function can’t assume the inputs are `Float64`, it must return a `ForwardDiff.Dual` when the inputs are `ForwardDiff.Dual`.

---

<div class="post-metadata">

**Author:** ![Laura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laura/32/207742_2.png) [@Laura](https://discourse.julialang.org/u/Laura)\
**Post date:** [April 19, 2024, 3:33pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/7 "2024-04-19T15:33:12Z")

</div>

Ah, I think I understand what you mean.

So if for example, my function made a new temporary matrix every itteration instead of storing those values in the p.matrix then it would all work. The reason I made it store the temporary values in p.matrix was to save on allocations if that makes sense.

Or alternatively I turn p.matrix into a Dual numer. But won’t that cause problems when I just try to evaluate the function on it’s own?

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [April 19, 2024, 5:47pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/8 "2024-04-19T17:47:04Z")

</div>

> [@Laura](#):
>
> So if for example, my function made a new temporary matrix every itteration instead of storing those values in the p.matrix then it would all work. The reason I made it store the temporary values in p.matrix was to save on allocations if that makes sense.
> 
> Or alternatively I turn p.matrix into a Dual numer. But won’t that cause problems when I just try to evaluate the function on it’s own?

Yes, you’re right. Allocating would solve it, and making `p.matrix` a `Dual` will create problems with ordinary evaluations.

This is what Chris suggested a solution to above. The `PreallocationTools` package (which I haven’t used) apparently can help with this problem.

---

<div class="post-metadata">

**Author:** ![JonasWickman](https://avatars.discourse-cdn.com/v4/letter/j/9de0a6/32.png) [@JonasWickman](https://discourse.julialang.org/u/JonasWickman)\
**Post date:** [April 19, 2024, 6:15pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/9 "2024-04-19T18:15:26Z")

</div>

It’s a little hard to tell exactly what your code does without any example data, but here is a (rather nonsensical) example of how to use PreallocationTools.jl with NonlinearSolve.jl:

```julia
using NonlinearSolve
using LinearAlgebra
using PreallocationTools

struct Parameters{T}
    matrix::T
end

function equation!(y, z, p::Parameters)

    # Cached Matrix 
    cmatrix = PreallocationTools.get_tmp(p.matrix, z)

    for i ∈ axes(cmatrix, 1), j ∈ axes(cmatrix, 2)
        cmatrix[i,j] = exp( z[i] + 2*z[j] )
    end

    for i ∈ eachindex(y)
        y[i] = 2.5 - sum( cmatrix[i,j]*z[j] for j ∈ eachindex(z) )
    end
    
    return nothing
end

u = zeros(3)
m = zeros(3,3)
cmatrix = PreallocationTools.DiffCache(m)
p = Parameters(cmatrix)

prob = NonlinearProblem(equation!, u, p)

sol = solve(prob)

```

---

<div class="post-metadata">

**Author:** ![Laura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laura/32/207742_2.png) [@Laura](https://discourse.julialang.org/u/Laura)\
**Post date:** [April 22, 2024, 8:03am UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/10 "2024-04-22T08:03:30Z")

</div>

This ended up almost working for me except before converting p.matrix using get\_tmp I first had to turn it into a DiffCache object.

> [@JonasWickman](#):
>
> It’s a little hard to tell exactly what your code does without any example data, but here is a (rather nonsensical) example of how to use [PreallocationTools.jl](https://juliahub.com/ui/Packages/General/PreallocationTools) with [NonlinearSolve.jl](https://juliahub.com/ui/Packages/General/NonlinearSolve):

Using my own example it would look like the code below. I will probably still fiffle around with it a bit because I have the feeling I could further optimize this by passing the p.matrix as a DiffCache from the beginning.

```julia
using NonlinearSolve
using LinearAlgebra
using JLD2
using PreallocationTools

mutable struct parameters{T}

    n::Int64
    norm::Vector{Float64}
    matrix::Matrix{T}
    tensor::Array{Float64,4}

end

function equation!(y, z, p::parameters)

    cache_matrix = PreallocationTools.get_tmp(DiffCache(p.matrix), z)

    @inbounds for i = 1:p.n
        for j = 1:i
            acc = zero(eltype(p.matrix))
            for k = 1:p.n
                acctmp = zero(eltype(p.matrix))
                @simd for l = 1:p.n
                    acctmp += p.tensor[l, k, j, i] * z[p.n+l]
                end
                acc += acctmp * z[k]
            end
            cache_matrix[i, p.n+j] = cache_matrix[j, p.n+i] = cache_matrix[p.n+i, j] = cache_matrix[p.n+j, i] = acc
        end
    end

    cache_matrix[p.n+1:2p.n, 1:p.n] .*= p.norm / (z[1] * z[p.n+1])
    cache_matrix[1:p.n, p.n+1:2p.n] .*= p.norm / (z[1] * z[p.n+1])

    y[1:2p.n] = (cache_matrix - I * z[2p.n+1]) * z[1:2p.n]

end

u = load_object(string("example u.jld2"))
p = load_object(string("example p.jld2"))

prob = NonlinearProblem(equation!, u, p)

sol = solve(prob)

```

---

<div class="post-metadata">

**Author:** ![JonasWickman](https://avatars.discourse-cdn.com/v4/letter/j/9de0a6/32.png) [@JonasWickman](https://discourse.julialang.org/u/JonasWickman)\
**Post date:** [April 22, 2024, 3:13pm UTC](https://discourse.julialang.org/t/my-function-works-on-its-own-but-not-in-nonlinearsolve/113204/11 "2024-04-22T15:13:41Z")

</div>

> [@Laura](#):
>
> Using my own example it would look like the code below. I will probably still fiffle around with it a bit because I have the feeling I could further optimize this by passing the p.matrix as a DiffCache from the beginning.

Right, you won’t get much benefit out of PreallocationTools if you’re recreating the DiffCache at every iteration. The following might work for you:

```julia
using Accessors
# Functions for the problem
# [...]
u = load_object(string("example u.jld2"))
p_in = load_object(string("example p.jld2"))
p = @set p_in.matrix = PreallocationTools.DiffCache(p_in.matrix)
prob = NonlinearProblem(equation!, u, p)
sol = solve(prob)

```

and changing the line getting the cache to:

```julia
cache_matrix = PreallocationTools.get_tmp(p.matrix, z)

```