Hello again, I’m solving iteratively this system of 6 equations as part of a larger simulation

## Toy data

```
using StaticArrays
using LinearAlgebra
mutable struct V{T}
u::T
A::T
A0::T
γ::T
β::T
end
struct N{T}
U::MVector{6, T}
W::MVector{3, T}
F::MVector{6, T}
J::MArray{Tuple{6, 6}, T, 2, 36}
end
function N()
U = @MArray zeros(6)
W = @MArray zeros(3)
F = @MArray zeros(6)
J = @MArray zeros(6,6)
J .+= I(6)
N(U, W, F, J)
end
function getVs(n::Int64)
u1 = zeros(n)
u2 = zeros(n)
u3 = zeros(n)
r1 = 0.018105
r2 = 0.015799
r3 = 0.008561
A1 = zeros(n) .+ pi * r1^2 .+ rand(n)
A2 = zeros(n) .+ pi * r2^2 .+ rand(n)
A3 = zeros(n) .+ pi * r3^2 .+ rand(n)
A01 = zeros(n) .+ pi * r1^2
A02 = zeros(n) .+ pi * r2^2
A03 = zeros(n) .+ pi * r3^2
σ = 0.5
Eh0 = 445.0
ρ = 1060.0
β1 = sqrt.(pi ./ A1) .* Eh0 / (1 - σ^2)
γ1 = β1 ./ (3ρ * r1 * sqrt(pi))
β2 = sqrt.(pi ./ A2) .* Eh0 / (1 - σ^2)
γ2 = β2 ./ (3ρ .* r2sqrt(pi))
β3 = sqrt.(pi ./ A3) .* Eh0 / (1 - σ^2)
γ3 = β3 ./ (3ρ .* r3 * sqrt(pi))
V(u1, A1, A01, γ1, β1), V(u2, A2, A02, γ2, β2), V(u3, A3, A03, γ3, β3)
end
```

## Utility functions

```
function w!(n::N, k)
n.W[1] = n.U[1] + 4 * k[1] * n.U[4]
n.W[2] = n.U[2] - 4 * k[2] * n.U[5]
n.W[3] = n.U[3] - 4 * k[3] * n.U[6]
end
function jac!(n::N, v1::V, v2::V, v3::V, k)
n.J[1, 4] = 4k[1]
n.J[2, 5] = -4k[2]
n.J[3, 6] = -4k[3]
n.J[4, 1] = (n.U[4] * n.U[4] * n.U[4] * n.U[4])
n.J[4, 2] = -(n.U[5] * n.U[5] * n.U[5] * n.U[5])
n.J[4, 3] = -(n.U[6] * n.U[6] * n.U[6] * n.U[6])
n.J[4, 4] = 4n.U[1] * (n.U[4] * n.U[4] * n.U[4])
n.J[4, 5] = -4n.U[2] * (n.U[5] * n.U[5] * n.U[5])
n.J[4, 6] = -4n.U[3] * (n.U[6] * n.U[6] * n.U[6])
n.J[5, 4] = 2v1.β[end] * n.U[4] / sqrt(v1.A0[end])
n.J[5, 5] = -2v2.β[1] * n.U[5] / sqrt(v2.A0[1])
n.J[6, 4] = 2v1.β[end] * n.U[4] / sqrt(v1.A0[end])
n.J[6, 6] = -2v3.β[1] * n.U[6] / sqrt(v3.A0[1])
end
function f!(n::N, v1::V, v2::V, v3::V, k)
n.F[1] = n.U[1] + 4k[1] * n.U[4] - n.W[1]
n.F[2] = n.U[2] - 4k[2] * n.U[5] - n.W[2]
n.F[3] = n.U[3] - 4k[3] * n.U[6] - n.W[3]
n.F[4] =
n.U[1] * (n.U[4] * n.U[4] * n.U[4] * n.U[4]) - n.U[2] * (n.U[5] * n.U[5] * n.U[5] * n.U[5]) -
n.U[3] * (n.U[6] * n.U[6] * n.U[6] * n.U[6])
n.F[5] =
v1.β[end] * (n.U[4] * n.U[4] / sqrt(v1.A0[end]) - 1) -
(v2.β[1] * (n.U[5] * n.U[5] / sqrt(v2.A0[1]) - 1))
n.F[6] =
v1.β[end] * (n.U[4] * n.U[4] / sqrt(v1.A0[end]) - 1) -
(v3.β[1] * (n.U[6] * n.U[6] / sqrt(v3.A0[1]) - 1))
end
function u!(n::N, v1::V, v2::V, v3::V)
n.U[1] = v1.u[end]
n.U[2] = v2.u[1]
n.U[3] = v3.u[1]
n.U[4] = sqrt(sqrt(v1.A[end]))
n.U[5] = sqrt(sqrt(v2.A[1]))
n.U[6] = sqrt(sqrt(v3.A[1]))
end
```

## Implementation

```
function alloc!(n::N, v1::V, v2::V, v3::V, k)
u!(n, v1, v2, v3)
w!(n, k)
jac!(n, v1, v2, v3, k)
f!(n, v1, v2, v3, k)
end
function NR!(n::N, v1::V, v2::V, v3::V, k)
while norm(n.F)>1e-5
n.U .+= n.J \ (-n.F)
w!(n, k)
f!(n, v1, v2, v3, k)
end
end
function update!(v1::V, v2::V, v3::V, n::N)
v1.u[end] = n.U[1]
v2.u[1] = n.U[2]
v3.u[1] = n.U[3]
v1.A[end] = n.U[4] * n.U[4] * n.U[4] * n.U[4]
v2.A[1] = n.U[5] * n.U[5] * n.U[5] * n.U[5]
v3.A[1] = n.U[6] * n.U[6] * n.U[6] * n.U[6]
end
function solve!(n::N, v1::V, v2::V, v3::V)
k = (sqrt(1.5*v1.γ[end]), sqrt(1.5*v2.γ[1]), sqrt(1.5*v3.γ[1]))
alloc!(n, v1, v2, v3, k)
NR!(n, v1, v2, v3, k)
update!(v1, v2, v3, n)
end
```

Benchmarking the solver I get

```
julia> v1,v2,v3=getVs(100);n=N();
julia> @btime solve!($n, $v1, $v2, $v3)
75.288 ns (0 allocations: 0 bytes)
```

and

```
julia> @btime alloc!($n, $v1, $v2, $v3, $k)
47.169 ns (0 allocations: 0 bytes)
julia> @btime NR!($n, $v1, $v2, $v3, $k)
14.059 ns (0 allocations: 0 bytes)
julia> @btime update!($v1, $v2, $v3, $n.U)
7.202 ns (0 allocations: 0 bytes)
```

therefore the *slow* part is the building of relevant vectors in `alloc!`

.

However this is not what I see when profiling the simulation as most of the time is spent in `NR!`

when calling `\`

in

```
n.U .+= n.J \ (-n.F)
```

Any advice on how can I go about optimising `\`

?