Hi

I am using a small code to generate some results. Now I want to take the average of different realizations. So I simply add a more for loop with loop number `Nr`

. But I find that when I increase `Nr`

, the code speed becomes slow. When `Nr = 100`

, the total time obtained by `@time`

is roughly 1500 seconds. But when I increase `Nr = 500`

, the speed is suddenly much slower (it runs more than four hours and haven’t finished yet.)

The code is attached below. I want to know the reason why the code speed decreases.

```
using LinearAlgebra
using StatsBase
using Statistics
using Random
using Plots
using RandomMatrices
function main()
theta = 1.39
#println(theta)
Nk = 800
Nr = 500
ps = 3
sum_z = zeros(Nk,1)
sum_zlog = zeros(Nk,1)
sum_s = zeros(Nk,1)
sum_slog = zeros(Nk,1)
for i in 1:Nr
Npool = Int(floor(10^(ps)))
pool = zeros(Npool, 1)
pool[:,1] .= 0.5
pool_sample = zeros(Npool,1)
z_list = zeros(Nk,1)
s_list = zeros(Nk,1)
z_log_list = zeros(Nk,1)
s_log_list = zeros(Nk,1)
Kup= [cos(theta) 0; 0 sin(theta)]
Kdown = [sin(theta) 0; 0 cos(theta)]
P2up = kron([1. 0.;0. 1.], [1 0; 0 0])
P2down = kron([1 0;0 1],[0 0;0 1])
h = Haar(2)
for k in 1:Nk
poolcount =0
while poolcount < Npool
z1 = pool[rand(1:Npool),1]
rho1 = diagm(0=>[z1,1-z1])
z2 = pool[rand(1:Npool),1]
rho2 = diagm(0=>[z2,1-z2])
u1 = rand(h,2)
u2 = rand(h,2)
K1up = u1*Kup*u1'
K1down = u1*Kdown*u1'
K2up = u2*Kup*u2'
K2down = u2*Kdown*u2'
rho1p = K1up*rho1*K1up'
rho2p = K2up*rho2*K2up'
p1 = real(tr(rho1p+rho1p'))/2
p2 = real(tr(rho2p+rho2p'))/2
rho1p = (rho1p+rho1p')/(2*p1)
rho2p = (rho2p+rho2p')/(2*p2)
rho = kron(rho1p,rho2p)
U = rand(h,4)
rho_p = P2up*U*rho*U'*P2up'
p = real(tr(rho_p+rho_p'))/2
temp =(rho_p+rho_p')/2
rho_f = [tr(temp[1:2,1:2]) tr(temp[1:2,3:4]);tr(temp[3:4,1:2]) tr(temp[3:4,3:4])]/(p)
rho_f = (rho_f+rho_f')/2
t = abs(tr(rho_f*rho_f))
z = (1-t)/(1+sqrt(2*t-1))
if abs(z)>0
poolcount = poolcount+1
pool_sample[poolcount,1] = abs(z)
end
end
pool = copy(pool_sample)
pool_sample[:,1] .= 0
z_log_list[k,1] = mean(log.(pool),dims =1)[1]
s_log_list[k,1] = std(log.(pool),dims = 1)[1]/sqrt(Npool)
z_list[k,1] = mean(pool,dims =1)[1]
s_list[k,1] = std(pool,dims = 1)[1]/sqrt(Npool)
end
sum_z += z_list
sum_zlog += z_log_list
sum_s += s_list.^2
sum_slog += s_log_list.^2
end
return vec(sum_z)/Nr, sqrt.(vec(sum_s))/Nr,vec(sum_zlog)/Nr, sqrt.(vec(sum_slog))/Nr
end
@time z0,s0, z1,s1 = main();
```