Hi there,

I am creating a function for the Bellman operator, a common tool in economics and macroeconomics. The function takes an Array of dimension n where n is the number of state variables and iterates it one period, output is therefore an array of dimension n as well.

I am new to Julia but I followed the performance tips I could understand and added `@view`

where applicable. Computation time has improved but it’s still too slow in my opinion, here’s the code

```
using Distributions
using Random
using BenchmarkTools
function T(V0::Array{Float64, 4}, cost_g::Matrix{Float64}, cost_v::Matrix{Float64}, n_grid_ω::Int64, n_grid_k::Int64, n_grid_c::Int64, P::Array{Float64, 3}, β::Float64, grid_Π::Matrix{Float64})
VP = zeros(n_grid_ω, n_grid_k, n_grid_c, n_grid_c)
EV_cg_cv = zeros(n_grid_ω, n_grid_k, n_grid_c, n_grid_c, 3)
for i in 1:n_grid_ω::Int64, r in 1:3, cg in 1:n_grid_c::Int64, cv in 1:n_grid_c::Int64
EV_cg_cv[i, :, cg, cv, r] = transpose(@view(P[i,:,r])) * @view(V0[:,:,cg,cv]) #Slow
end
EV = zeros(n_grid_ω, n_grid_k, 3)
EV = dropdims(sum(EV_cg_cv, dims = [3;4])/(n_grid_c)^2, dims = tuple(3,4))
for cg in 1:n_grid_c::Int64, cv in 1:n_grid_c::Int64, k in 1:n_grid_k::Int64
cost = [0, cost_g[cg, k], cost_v[cv, k]]
for i in 1:n_grid_ω
VP[i, k, cg, cv] = maximum(β* EV[i,k,:] - cost) .+ @view(grid_Π[i, k]) #Slow
end
end
return(VP)
end
nω, nk, nc = [8, 9, 10];
Random.seed!(2022);
V = rand(nω, nk, nc, nc);
cost = rand(nc, nk);
grid_Π = rand(nω, nk);
P = randn(nω, nω, 3);
julia> @btime T(V, cost, cost, nω, nk, nc, P, 0.9, grid_Π);
2.318 ms (33333 allocations: 3.22 MiB)
```

Does this appear reasonable for a problem of this size or is it too slow? I know from profiling that much time is spent on the lines marked as #slow. But these 2 lines are also computed on 4 iterative dimensions so it could be normal after all.

I read that the matrix representation for data (common in Matlab) is not so appropriate in Julia. Anything in that direction ?

Thanks