Julia中用SharedArray与pmap实现分布式计算遇错求助

核心错误点与修复 #

1. 函数名不匹配 #

定义的Haar随机矩阵生成函数为rand_haar2，但代码中调用的是未定义的rand_haar2_slower，直接触发未定义函数错误。

• 修复：将所有rand_haar2_slower替换为rand_haar2。

2. 结果拼接逻辑错误 #

pmap返回的每个元素是长度为1e5的数组，但代码错误地仅取每个数组的第一个元素ret[i][1]，导致pool_target长度仅为10，远小于预期的1e6，后续计算会出现维度不匹配。

• 修复：直接拼接完整数组：

  pool_target = reduce(vcat, ret)
  

3. SharedArray类型不匹配 #

循环中pool = pool_target将SharedArray替换为普通Array，但pool_calc参数要求为SharedArray，下一轮循环会触发类型错误。

• 修复：将生成的pool_target转换为SharedArray：

  pool = SharedArray(pool_target)
  

4. 冗余变量与越界风险 #

• Nk参数定义后未使用，循环固定为1到800，需将循环改为for k in 1:Nk避免S_curve/S_var越界；
• pool_sample变量定义后未使用，直接删除；

5. 数值稳定性问题 #

当pool_target中出现0值时，log(0)会产生-Inf，破坏熵计算结果。

• 修复：添加微小偏移量限制数值范围：

  eps_val = eps(Double64)
  pool_target_clamped = clamp.(pool_target, eps_val, 1-eps_val)
  spool .= -pool_target_clamped .* log.(pool_target_clamped) .- (1.0 .- pool_target_clamped) .* log1p.(-pool_target_clamped)
  

6. 随机种子优化 #

Random.seed!(myid())可能导致不同进程的随机序列关联性过高，建议使用更独立的种子：

• 修复：替换为：

  Random.seed!(myid() + rand(UInt32))
  

修正后的完整代码 #

using Distributed
addprocs(10)

@everywhere using LinearAlgebra
@everywhere using Statistics
@everywhere using DoubleFloats
@everywhere using StaticArrays
@everywhere using SharedArrays
@everywhere using JLD
@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() + rand(UInt32))

    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(Val{2}())
        u2 = rand_haar2(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(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 += 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)
    spool = zeros(Double64,Npool)

    pool .= 0.5
    eps_val = eps(Double64)

    for k in 1:Nk
        # 10个进程各生成1e5样本，合计1e6个
        ret = pmap(Np->pool_calc(theta=theta, pool=pool, Np=Np), fill(10^5, 10))
        pool_target = reduce(vcat, ret)
        
        # 数值稳定性处理：避免对数输入为0或1
        pool_target_clamped = clamp.(pool_target, eps_val, 1-eps_val)
        spool .= -pool_target_clamped .* log.(pool_target_clamped) .- (1.0 .- pool_target_clamped) .* log1p.(-pool_target_clamped)
        
        S_curve[k] = mean(spool)
        S_var[k] = (std(spool)/sqrt(Npool))^2

        # 转换为SharedArray供下一轮使用
        pool = SharedArray(pool_target)
    end

    label = @sprintf "%.3f" Float32(theta)
    save("entropy_real_128p_$(label)_ps6.jld", "s", S_curve, "t", S_var)
end

main();

额外优化建议 #

• 移除未使用的包：StatsBase、Dates未在代码中用到，可从@everywhere using中删除，减少进程初始化负载；
• 预编译加速：在@everywhere块中添加@precompile_all_calls，加快首次运行速度；
• 内存优化：pool_store可改用Vector{Double64}(undef, Np)直接赋值，避免初始化零数组的开销；
• 并行策略调整：CPU密集型任务可考虑使用@distributed或多线程Threads.@threads替代pmap，降低进程间通信开销。

内容的提问来源于stack exchange，提问作者jisutich