Hi
I am thinking about using SharedArrays in the following problem. Suppose I have an index k which increases from 1 to 800 (for instance). And for each k, I have a pool p which has large size and many numbers stored in it. I want to get p_k recursively. The protocol is like, if I know p_{k-1}, then I can randomly choose two values z^{1,2}_{k-1} from it and get a new value through a function f like z_k = f(z^1_{k-1},z^2_{k-1}). Then I store it into p_k and repeated this many times until p_k is full and then I try to get p_{k+1} from p_k.
Due to the large size of the pool, I try to use parallel computation to speed up my code. I am trying to use pmap and use a SharedArray as my p_{k-1} within each k. So I write the following code
using Distributed
addprocs(10)
@everywhere using LinearAlgebra
@everywhere using StatsBase
@everywhere using Statistics
@everywhere using DoubleFloats
@everywhere using StaticArrays
@everywhere using SharedArrays
@everywhere using JLD
@everywhere using Dates
@everywhere using Random
@everywhere using Printf
@everywhere function rand_haar2(::Val{n}) where n
M = @SMatrix randn(ComplexDF64, n,n)
q = qr(M).Q
L = cispi.(2 .* @SVector(rand(Double64,n)))
return q*diagm(L)
end
@everywhere function pool_calc(theta,pool::SharedArray,Np)
Random.seed!(myid())
pool_store = zeros(Double64,Np)
Kup= @SMatrix[Double64(cos(theta)) 0; 0 Double64(sin(theta))]
Kdown = @SMatrix[Double64(sin(theta)) 0; 0 Double64(cos(theta))]
P2up = kron(@SMatrix[Double64(1.) 0.;0. 1.], @SMatrix[1 0; 0 0])
P2down = kron(@SMatrix[Double64(1) 0;0 1],@SMatrix[0 0;0 1])
poolcount = 0
poolsize = length(pool)
while poolcount < Np
z1 = pool[rand(1:poolsize)]
rho1 = diagm(@SVector[z1,1-z1])
z2 = pool[rand(1:poolsize)]
rho2 = diagm(@SVector[z2,1-z2])
u1 = rand_haar2_slower(Val{2}())
u2 = rand_haar2_slower(Val{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
if rand()<p1
rho1p = (rho1p+rho1p')/(2*p1)
else
rho1p = K1down*rho1*K1down'/((1-p1))
end
if rand()<p2
rho2p = (rho2p+rho2p')/(2*p2)
else
rho2p = K2down*rho2*K2down'/((1-p2))
end
rho = kron(rho1p,rho2p)
U = rand_haar2_slower(Val{4}())
rho_p = P2up*U*rho*U'*P2up'
p = real(tr(rho_p+rho_p'))/2
if rand()<p
temp =(rho_p+rho_p')/2
rho_f = @SMatrix[temp[1,1]+temp[2,2] temp[1,3]+temp[2,4]; temp[3,1]+temp[4,2] temp[3,3]+temp[4,4]]/(p)
else
temp = P2down*U*rho*U'*P2down'
rho_f = @SMatrix[temp[1,1]+temp[2,2] temp[1,3]+temp[2,4]; temp[3,1]+temp[4,2] temp[3,3]+temp[4,4]]/(1-p)
end
rho_f = (rho_f+rho_f')/2
t = abs(tr(rho_f*rho_f))
z = (1-t)/(1+abs(sqrt(2*t-1)))
if !iszero(abs(z))
poolcount = poolcount+1
pool_store[poolcount] = abs(z)
end
end
return pool_store
end
function main()
theta = parse(Double64,ARGS[1])
Nk = parse(Int,ARGS[2])
S_curve = zeros(Double64,Nk)
S_var = zeros(Double64,Nk)
Npool = Int(floor(10^6))
pool = SharedArray{Double64}(Npool)
pool_sample = zeros(Double64,Npool)
spool = zeros(Double64,Npool)
pool .=0.5
for k =1:800
ret = pmap(Np->pool_calc(theta = theta,pool=pool,Np=Np),fill(10^5,10))
pool_target = reduce(vcat,[ret[i][1] for i = 1:10])
spool .=-pool_target .*log.(pool_target).-(1.0 .- pool_target).*log1p.(-pool_target)
S_curve[k] = mean(spool)
S_var[k] = (std(spool)/sqrt(Npool))^2
pool = pool_target
end
label = @sprintf "%.3f" Float32(theta)
save("entropy_real_128p_$(label)_ps6.jld","s", S_curve, "t", S_var)
end
main();
But I faced an error
How to solve this problem?
Thanks