Julia crashes without reporting anything when I optimize a vector of models in parallel

[WalterMadelim]

My setting is

julia> Threads.nthreads()
4

I can’t figure out what is happening here. The behavior is

julia> parallel_CG!(B, θ, β, μ, ν)
┌ Info: before entering @threads
│   sub_j_vec =
│    4-element Vector{Int64}:
│     1
│     4
│     2
└     3
┌ Info: inside @threads
└   j = 2
┌ Info: inside @threads
└   j = 1
┌ Info: inside @threads
└   j = 3
PS K:\uc24> 

The Info are some debugging message that I intend to let julia display. As you can see, the julia crashes so that the prompt reverts to my windows powershell.

The code to reproduce this is (copy-paste to REPL)

import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import JuMP.value as ı
import LinearAlgebra.⋅ as ⋅
import LinearAlgebra.norm as norm
import Random

function prev(t, d, T) return (n = t - d; n<1 ? T+n : n) end; # used in h[:U]
function gen_a_D_house(ID, Drng) return (
    ID = ID,
    CND = .5rand(1:7), # thermal conductance of the building wall
    INR = rand(6:20), # thermal inertia of the building
    COP = rand(2:.5:4), # how many Watts cooling power from 1-Watt e-power
    Q_I = rand(3:9, T), # heat generated by indoor human activity, server etc.
    Q_BUS = rand(25:35), # RateA of the cooling water/air pipeline
    OH = rand(24:29), # the hottest indoor degree Celsius
    OΔ = rand(4:9), # the coldest indoor degree Celsius = OH - OΔ
    D = rand(Drng, T), # base demand
    P_BUS = rand(12:33), # unidirectional due to no renewables
    P_EV = [2, rand(3:7)], # MIN and MAX charging power
    P_AC = rand(6:9), # AC's cooling power is limited
    E_EV = rand(10:39),
    U = [rand(1:4, rand(2:5)) for _ = 1:rand(1:4)] # an entry is a cycle vector of _an_ uninterruptible load 
) end; function gen_a_G_house(ID, Drng) return (gen_a_D_house(ID, Drng)..., G = rand(3:17, T)) end;
function solve_compact_formulation()
    x = similar(D, NamedTuple); # store the VariableRefs in NamedTuple's
    model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # monolithic
    JuMP.set_attribute(model, "MIPGap", 0);
    JuMP.set_attribute(model, "MIPGapAbs", 0);
    for (j, block) = enumerate(D) # operational detail for each individual block
        if length(block) == 1
            h = block[1]
            bEV, pEV = JuMP.@variable(model, [1:T], Bin), JuMP.@variable(model, [1:T]); JuMP.@constraint(model, sum(pEV) ≥ h[:E_EV])
            JuMP.@constraint(model, h[:P_EV][begin]bEV ≤ pEV); JuMP.@constraint(model, pEV ≤ h[:P_EV][end]bEV)
            bU = JuMP.@variable(model, [t = 1:T, i = eachindex(h[:U])], Bin); JuMP.@constraint(model, sum(bU; dims = 1) .≥ true)
            pU = JuMP.@expression(model, [t=1:T], sum(sum(bU[prev(t,φ-1,T), i]P for (φ,P) = enumerate(v)) for (i,v) = enumerate(h[:U]))) # Uninterruptible demand
            q = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = h[:Q_BUS]) # direct cooling power supply
            pAC = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = h[:P_AC]) # we consider only cooling mode in this instance
            o = JuMP.@variable(model, [1:T], lower_bound = h[:OH]-h[:OΔ], upper_bound = h[:OH]) # indoor temperature (state variable)
            JuMP.@constraint(model, [t=1:T], (O[t]-o[t])h[:CND] + h[:Q_I][t] -q[t] -h[:COP]pAC[t] == (o[t<T ? t+1 : 1]-o[t])h[:INR]) # dynamic equation
            pBus = JuMP.@variable(model, [1:T], lower_bound = 0, upper_bound = h[:P_BUS]) # unidirectional for NO-DER household
            JuMP.@constraint(model, pBus ≥ h[:D] + pEV + pU + pAC);
        else
            bLent, pLent = JuMP.@variable(model, [1:T], Bin), JuMP.@variable(model, [1:T])
            bEV, pEV = JuMP.@variable(model, [1:T, 1:2], Bin), JuMP.@variable(model, [1:T, 1:2])
            for (i, h) = enumerate(block) # EV's part
                JuMP.@constraint(model, h[:P_EV][begin]view(bEV, :, i) ≤ view(pEV, :, i))
                JuMP.@constraint(model, view(pEV, :, i) ≤ h[:P_EV][end]view(bEV, :, i))
                if i == 1 # the household that has RES
                    JuMP.@constraint(model, (1 .- bLent)h[:P_EV][begin] ≤ pLent)
                    JuMP.@constraint(model, pLent ≤ (1 .- bLent)h[:P_EV][end])
                    JuMP.@constraint(model, h[:E_EV] ≤ sum(view(pEV, :, i)))
                else
                    JuMP.@constraint(model, h[:E_EV] ≤ sum(view(pEV, :, i)) + sum(pLent))
                end
                JuMP.@constraint(model, view(bEV, :, i) ≤ bLent) # bLent = 0 means lending
            end
            bU = JuMP.@variable(model, [i = 1:2, j = eachindex(block[i][:U]), t = 1:T], Bin) # ⚠️ This construct is a bit involved
            JuMP.@constraint(model, [i = 1:2, j = eachindex(block[i][:U])], sum(bU[i, j, :]) .≥ true)
            pU = JuMP.@expression(model, [i = 1:2, t = 1:T], sum(sum(bU[i, j, prev(t,φ-1,T)]P for (φ,P) = enumerate(v)) for (j,v) = enumerate(block[i][:U])))
            q = JuMP.@variable(model, [1:T, i = 1:2], lower_bound = 0, upper_bound = block[i][:Q_BUS])
            pAC = JuMP.@variable(model, [1:T, i = 1:2], lower_bound = 0, upper_bound = block[i][:P_AC])
            o = JuMP.@variable(model, [1:T, i = 1:2], lower_bound = block[i][:OH]-block[i][:OΔ], upper_bound = block[i][:OH])
            JuMP.@constraint(model, [t = 1:T, i = 1:2], 
                (O[t]-o[t, i])block[i][:CND] + block[i][:Q_I][t] -q[t, i] -block[i][:COP]pAC[t, i] == (o[t<T ? t+1 : 1, i]-o[t, i])block[i][:INR]
            ); # dynamic equation
            pBus = JuMP.@variable(model, [1:T, i = 1:2], lower_bound = i>1 ? 0 : -block[i][:P_BUS], upper_bound = block[i][:P_BUS])
            for (i, h) = enumerate(block)
                if i == 1
                    JuMP.@constraint(model, view(pBus, :, i) + h[:G] ≥ pLent + h[:D] + view(pEV, :, i) + view(pU, i, :) + view(pAC, :, i))
                else
                    JuMP.@constraint(model, view(pBus, :, i)         ≥         h[:D] + view(pEV, :, i) + view(pU, i, :) + view(pAC, :, i))
                end
            end
        end
        x[j] = (
            q    = q,   
            o    = o,   
            pBus = pBus,
            bEV  = bEV, 
            pEV  = pEV, 
            pAC  = pAC, 
            bU   = bU
        ); length(block) > 1 && (x[j] = (bLent = bLent, pLent = pLent, x[j]...))
    end
    JuMP.@expression(model, pAgr, sum(haskey(x[j], :bLent) ? vec(sum(x[j][:pBus]; dims = 2)) : x[j][:pBus] for j = eachindex(x)))
    JuMP.@expression(model, qAgr, sum(haskey(x[j], :bLent) ? vec(sum(x[j][:q]; dims = 2)) : x[j][:q] for j = eachindex(x)))
    JuMP.@constraint(model, β[t = 1:T], P_AGR ≥ pAgr[t]) # 🟡 Linking (note that we require positiveness thought not explicit here)
    JuMP.@constraint(model, μ[t = 1:T], Q_AGR ≥ qAgr[t]) # 🟡 Linking
    JuMP.@constraint(model, ν, E_AGR ≥ sum(qAgr)) # 🟡 Linking [ice-storage is limited]
    JuMP.@objective(model, Min, MarketPrice ⋅ pAgr)
    JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
    println("compact one-shot solved in $(JuMP.solve_time(model)) seconds")
    ı.(x[1][:bU]);
    return an_UB = JuMP.objective_value(model) + 200
end;

function initialize_out(an_UB) # an_UB ⚠️ should be valid
    model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # outer in DantzigWolfe decomposition
    JuMP.set_silent(model)
    JuMP.@variable(model, β[1:T] ≥ 0) # e_power per t
    JuMP.@variable(model, μ[1:T] ≥ 0) # q_power per t
    JuMP.@variable(model, ν ≥ 0) # q_energy scalar
    JuMP.@variable(model, θ[1:J]) # This is a concise formulation---the cut is info-intensive
    JuMP.@expression(model, common, -sum(β)P_AGR -sum(μ)Q_AGR -(ν)E_AGR)
    JuMP.@expression(model, out_obj_tbMax, common + sum(θ))
    JuMP.@objective(model, Max, out_obj_tbMax)
    JuMP.@constraint(model, out_obj_tbMax ≤ an_UB)
    return model, θ, β, μ, ν
end; function initialize_inn(D)
    inn = Vector{JuMP.Model}(undef, J); # the intrablock MIP subproblems
    for (j, block) = enumerate(D)
        model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
        JuMP.set_silent(model)
        JuMP.set_attribute(model, "Threads", 1)
        JuMP.set_attribute(model, "TimeLimit", 3)
        JuMP.set_attribute(model, "MIPGap", 0)
        JuMP.set_attribute(model, "MIPGapAbs", 0)
        if length(block) == 1 # [pBus1; pBus2; pLent; vec(pEV); bLent; vec(bEV); vec(bDW)]
            h = block[1]
            JuMP.@variable(model, bEV[1:T], Bin); JuMP.@variable(model, pEV[1:T]); JuMP.@constraint(model, sum(pEV) ≥ h[:E_EV])
            JuMP.@constraint(model, h[:P_EV][begin]bEV ≤ pEV); JuMP.@constraint(model, pEV ≤ h[:P_EV][end]bEV)
            JuMP.@variable(model, bU[t = 1:T, i = eachindex(h[:U])], Bin); JuMP.@constraint(model, sum(bU; dims = 1) .≥ true)
            JuMP.@expression(model, pU[t=1:T], sum(sum(bU[prev(t,φ-1,T), i]P for (φ,P) = enumerate(v)) for (i,v) = enumerate(h[:U]))) # Uninterruptible demand
            JuMP.@variable(model, q[1:T], lower_bound = 0, upper_bound = h[:Q_BUS]) # direct cooling power supply
            JuMP.@variable(model, pAC[1:T], lower_bound = 0, upper_bound = h[:P_AC]) # we consider only cooling mode in this instance
            JuMP.@variable(model, o[1:T], lower_bound = h[:OH]-h[:OΔ], upper_bound = h[:OH]) # indoor temperature (state variable)
            JuMP.@constraint(model, [t=1:T], (O[t]-o[t])h[:CND] + h[:Q_I][t] -q[t] -h[:COP]pAC[t] == (o[t<T ? t+1 : 1]-o[t])h[:INR]) # dynamic equation
            JuMP.@variable(model, pBus[1:T], lower_bound = 0, upper_bound = h[:P_BUS]) # unidirectional for NO-DER household
            JuMP.@constraint(model, pBus ≥ h[:D] + pEV + pU + pAC);
            JuMP.@expression(model, prim_obj, MarketPrice ⋅ pBus) # some penalty terms still needed
        else
            JuMP.@variable(model, bLent[1:T], Bin); JuMP.@variable(model, pLent[1:T])
            JuMP.@variable(model, bEV[1:T, 1:2], Bin); JuMP.@variable(model, pEV[1:T, 1:2])
            for (i, h) = enumerate(block) # EV's part
                JuMP.@constraint(model, h[:P_EV][begin]view(bEV, :, i) ≤ view(pEV, :, i))
                JuMP.@constraint(model, view(pEV, :, i) ≤ h[:P_EV][end]view(bEV, :, i))
                if i == 1 # the household that has RES
                    JuMP.@constraint(model, (1 .- bLent)h[:P_EV][begin] ≤ pLent)
                    JuMP.@constraint(model, pLent ≤ (1 .- bLent)h[:P_EV][end])
                    JuMP.@constraint(model, h[:E_EV] ≤ sum(view(pEV, :, i)))
                else
                    JuMP.@constraint(model, h[:E_EV] ≤ sum(view(pEV, :, i)) + sum(pLent))
                end
                JuMP.@constraint(model, view(bEV, :, i) ≤ bLent) # bLent = 0 means lending
            end
            JuMP.@variable(model, bU[i = 1:2, j = eachindex(block[i][:U]), t = 1:T], Bin) # ⚠️ This construct is a bit involved
            JuMP.@constraint(model, [i = 1:2, j = eachindex(block[i][:U])], sum(bU[i, j, :]) .≥ true)
            JuMP.@expression(model, pU[i = 1:2, t = 1:T], sum(sum(bU[i, j, prev(t,φ-1,T)]P for (φ,P) = enumerate(v)) for (j,v) = enumerate(block[i][:U])))
            JuMP.@variable(model, q[1:T, i = 1:2], lower_bound = 0, upper_bound = block[i][:Q_BUS])
            JuMP.@variable(model, pAC[1:T, i = 1:2], lower_bound = 0, upper_bound = block[i][:P_AC])
            JuMP.@variable(model, o[1:T, i = 1:2], lower_bound = block[i][:OH]-block[i][:OΔ], upper_bound = block[i][:OH])
            JuMP.@constraint(model, [t = 1:T, i = 1:2], 
                (O[t]-o[t, i])block[i][:CND] + block[i][:Q_I][t] -q[t, i] -block[i][:COP]pAC[t, i] == (o[t<T ? t+1 : 1, i]-o[t, i])block[i][:INR]
            ); # dynamic equation
            JuMP.@variable(model, pBus[1:T, i = 1:2], lower_bound = i>1 ? 0 : -block[i][:P_BUS], upper_bound = block[i][:P_BUS])
            for (i, h) = enumerate(block)
                if i == 1
                    JuMP.@constraint(model, view(pBus, :, i) + h[:G] ≥ pLent + h[:D] + view(pEV, :, i) + view(pU, i, :) + view(pAC, :, i))
                else
                    JuMP.@constraint(model, view(pBus, :, i)         ≥         h[:D] + view(pEV, :, i) + view(pU, i, :) + view(pAC, :, i))
                end
            end
            JuMP.@expression(model, prim_obj, MarketPrice ⋅ sum(pBus; dims = 2)) # some penalty terms still needed
        end; inn[j] = model # save
    end # not include: penalty objective, non-block linking constr
    return inn
end; function reset_mj_obj!(mj, Β, Μ, Ν)
    pair = haskey(mj, :bLent)
    pBus, q = mj[:pBus], mj[:q]
    term_beta = JuMP.@expression(mj, pair ? sum((pBus[t,1] + pBus[t,2])Β[t] for t=1:T) : Β ⋅ pBus)
    term_mu = JuMP.@expression(mj, pair ? sum((q[t,1] + q[t,2])Μ[t] for t=1:T) : Μ ⋅ q)
    term_nu = JuMP.@expression(mj, sum(q)Ν)
    JuMP.@objective(mj, Min, +(mj[:prim_obj], term_beta, term_mu, term_nu))
end; function build_con(j, mj, θ, β, μ, ν)
    model = JuMP.owner_model(ν)
    pair = haskey(mj, :bLent)
    pBus, q = ı.(mj[:pBus]), ı.(mj[:q])
    term_beta = JuMP.@expression(model, pair ? sum((pBus[t,1] + pBus[t,2])β[t] for t=1:T) : β ⋅ pBus)
    term_mu = JuMP.@expression(model, pair ? sum((q[t,1] + q[t,2])μ[t] for t=1:T) : μ ⋅ q)
    term_nu = JuMP.@expression(model, sum(q)ν)
    cut_term = JuMP.@expression(model, +(ı(mj[:prim_obj]), term_beta, term_mu, term_nu))
    return JuMP.@build_constraint(θ[j] ≤ cut_term)
end; function get_full_num_tuple(mj)
    x = (
        q = ı.(mj[:q]),
        o = ı.(mj[:o]),
        pBus = ı.(mj[:pBus]),
        bEV = round.(Bool, ı.(mj[:bEV])),
        pEV = ı.(mj[:pEV]),
        pAC = ı.(mj[:pAC]),
        bU = round.(Bool, ı.(mj[:bU]))
    );  haskey(mj, :bLent) && (x = (
        bLent = round.(Bool, ı.(mj[:bLent])),
        pLent = ı.(mj[:pLent]),
        x...
    )); return x
end; function batch(J, TH) return (J-1) ÷ TH + 1 end;
function parallel_CG!(B, θ, β, μ, ν)
    COT = 1e-11
    model = JuMP.owner_model(ν)
    for ite = Base.OneTo(typemax(Int))
        JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false, dual = true)
        bound_of_out = JuMP.objective_bound(model) # Record
        Β, Μ, Ν = ı.(β), ı.(μ), ı(ν) # ✂️
        rand_j_vec, BATCH = Random.shuffle(Base.OneTo(J)), batch(J, TH)
        cut_found_flag = false
        for ba = Base.OneTo(BATCH) # do one partial parallel MIPs solving at an iteration
            sub_j_vec = rand_j_vec[1+(ba-1)TH : min(J, (ba)TH)]
            @info "before entering @threads" sub_j_vec
            Threads.@threads for j = sub_j_vec # strictly This is NOT a loop
                mj = inn[j]
                reset_mj_obj!(mj, Β, Μ, Ν)
                JuMP.optimize!(mj); # JuMP.assert_is_solved_and_feasible(mj; allow_local = false) # easy MIPs
                @info "inside @threads" j
            end
            error("#1 can you see this error??")
            Δ_vec = [ı(θ[j]) - JuMP.objective_value(inn[j]) for j = sub_j_vec]
            Δ, ii = findmax(Δ_vec)
            if Δ > COT
                j = sub_j_vec[ii]
                mj = inn[j]
                intFullvertex = get_full_num_tuple(mj)
                if !haskey(B[j], intFullvertex)
                    con = build_con(j, mj, θ, β, μ, ν)
                    cut = JuMP.add_constraint(model, con)
                    B[j][intFullvertex] = cut
                    @info "ite = $ite, cut's_maxvio = $Δ, out's bound = $bound_of_out"
                    cut_found_flag = true
                    break
                end
            end
        end
        cut_found_flag && continue
        break
    end
end;

seed = abs(rand(Int))
seed = 2004686255641684212
Random.seed!(seed)

TH = Threads.nthreads(:default)
T = 24; # global data
MarketPrice = [2, 1, 1, 1, 1, 1, 3, 5, 6, 7, 8, 9, 9, 8, 8, 9, 10, 11, 9, 6, 5, 3, 2, 2];
P_AGR = 47;
Q_AGR = 27;
E_AGR = 200;
O = [27.0, 27.2, 27.5, 27.8, 28.5, 30.0, 32.5, 35.0, 37.5, 39.5, 41.0, 42.0, 41.5, 40.5, 39.0, 37.0, 34.5, 32.0, 30.0, 29.0, 28.0, 27.5, 27.2, 27.0]; # global data: outdoor temperature

D = Vector{Vector{NamedTuple}}(undef, 4); # generate all the households
D[1] = [gen_a_G_house(1, 1:3), gen_a_D_house(2, 0:5)];
D[2] = [gen_a_D_house(3, 0:5)];
D[3] = [gen_a_D_house(4, 0:6)];
D[4] = [gen_a_D_house(5, 0:5)];
J = length(D);

an_UB = solve_compact_formulation();

model, θ, β, μ, ν = initialize_out(an_UB);
inn = initialize_inn(D); # haskey(inn[1], :bLent) == true
# my_lock = Threads.ReentrantLock();

B = [Dict{NamedTuple, JuMP.ConstraintRef}() for j = 1:J]; # a vec of Bijections::NamedTuple

parallel_CG!(B, θ, β, μ, ν)

[WalterMadelim]

The crux is at this block

@WalterMadelim

Threads.@threads for j = sub_j_vec # strictly This is NOT a loop
                mj = inn[j]
                reset_mj_obj!(mj, Β, Μ, Ν)
                JuMP.optimize!(mj); # JuMP.assert_is_solved_and_feasible(mj; allow_local = false) # easy MIPs
                @info "inside @threads" j
            end

I’m unsure if this style is correct.
But here is one more weird observation: if I replace this parallelization with

Threads.@threads for i = eachindex(sub_j_vec)
                j = sub_j_vec[i]
                mj = inn[j]
                reset_mj_obj!(mj, Β, Μ, Ν)
                JuMP.optimize!(mj); # JuMP.assert_is_solved_and_feasible(mj; allow_local = false) # easy MIPs
                @info "inside @threads" j
            end

I see

julia> parallel_CG!(B, θ, β, μ, ν)
┌ Info: before entering @threads
│   sub_j_vec =
│    4-element Vector{Int64}:
│     1
│     4
│     2
└     3
┌ Info: inside @threads
└   j = 4
┌ Info: inside @threads
└   j = 4
┌ Info: inside @threads
└   j = 4
ERROR: TaskFailedException

    nested task error: Gurobi Error 10017: Attempted to access a model while optimization is in progress
    Stacktrace:
      [1] _check_ret
        @ K:\judepot1115\packages\Gurobi\qCzRR\src\MOI_wrapper\MOI_wrapper.jl:417 [inlined]
      [2] _update_if_necessary(model::Gurobi.Optimizer; force::Bool, check_attribute_change::Bool)
        @ Gurobi K:\judepot1115\packages\Gurobi\qCzRR\src\MOI_wrapper\MOI_wrapper.jl:604
      [3] _update_if_necessary
        @ K:\judepot1115\packages\Gurobi\qCzRR\src\MOI_wrapper\MOI_wrapper.jl:575 [inlined]
      [4] optimize!(model::Gurobi.Optimizer)
        @ Gurobi K:\judepot1115\packages\Gurobi\qCzRR\src\MOI_wrapper\MOI_wrapper.jl:2793
      [5] optimize!
        @ K:\judepot1115\packages\MathOptInterface\JdLrc\src\Bridges\bridge_optimizer.jl:367 [inlined]
      [6] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{…})
        @ MathOptInterface.Utilities K:\judepot1115\packages\MathOptInterface\JdLrc\src\Utilities\cachingoptimizer.jl:379
      [7] optimize!(model::JuMP.Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})  
        @ JuMP K:\judepot1115\packages\JuMP\RGIK3\src\optimizer_interface.jl:595
      [8] optimize!(model::JuMP.Model)
        @ JuMP K:\judepot1115\packages\JuMP\RGIK3\src\optimizer_interface.jl:546
      [9] macro expansion
        @ .\REPL[10]:23 [inlined]
     [10] (::var"#1339#threadsfor_fun#127"{var"#1339#threadsfor_fun#125#128"{…}})(tid::Int64; onethread::Bool)
        @ Main .\threadingconstructs.jl:253
     [11] #1339#threadsfor_fun
        @ .\threadingconstructs.jl:220 [inlined]
     [12] (::Base.Threads.var"#1#2"{var"#1339#threadsfor_fun#127"{var"#1339#threadsfor_fun#125#128"{…}}, Int64})()
        @ Base.Threads .\threadingconstructs.jl:154
Stacktrace:
 [1] threading_run(fun::var"#1339#threadsfor_fun#127"{var"#1339#threadsfor_fun#125#128"{Float64, Vector{…}, Vector{…}, Vector{…}, Base.OneTo{…}}}, static::Bool)
   @ Base.Threads .\threadingconstructs.jl:173
 [2] macro expansion
   @ .\threadingconstructs.jl:190 [inlined]
 [3] parallel_CG!(B::Vector{Dict{…}}, θ::Vector{JuMP.VariableRef}, β::Vector{JuMP.VariableRef}, μ::Vector{JuMP.VariableRef}, ν::JuMP.VariableRef)
   @ Main .\REPL[10]:19
 [4] top-level scope
   @ REPL[32]:1
Some type information was truncated. Use `show(err)` to see complete types.

I think this is even weirder—why is j = 4 printed for 3 times?

[odow]

Gurobi environments are not threadsafe:

GitHub - jump-dev/Gurobi.jl: A Julia interface to the Gurobi Optimizer

[WalterMadelim]

There are two issues:

1. the hyperlink under this section GitHub - jump-dev/Gurobi.jl: A Julia interface to the Gurobi Optimizer gives 404 when it is clicked.
2. it seems that the same problem persists:

If you copy my original code to a text editor, and do this replacement
(ctrl + H)
replace
JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
with
JuMP.Model(Gurobi.Optimizer)
(there are 3 changes to be made)

And then execute the revised code, I observe this

julia> parallel_CG!(B, θ, β, μ, ν)
┌ Info: before entering @threads
│   sub_j_vec =
│    4-element Vector{Int64}:
│     1
│     4
│     2
└     3

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x7ffa226fe3b0 -- GRBdelconstrs at D:\gurobi1202\win64\bin\gurobi120.DLL (unknown line)
in expression starting at REPL[32]:1
GRBdelconstrs at D:\gurobi1202\win64\bin\gurobi120.DLL (unknown line)
GRBsetdblparam at K:\judepot1115\packages\Gurobi\qCzRR\src\gen120\libgrb_api.jl:2803 [inlined]
_set_param at K:\judepot1115\packages\Gurobi\qCzRR\src\MOI_wrapper\MOI_wrapper.jl:725
set at K:\judepot1115\packages\Gurobi\qCzRR\src\MOI_wrapper\MOI_wrapper.jl:714
jfptr_set_4372 at K:\judepot1115\compiled\v1.11\Gurobi\do9v6_pSn93.dll (unknown line)
empty! at K:\judepot1115\packages\Gurobi\qCzRR\src\MOI_wrapper\MOI_wrapper.jl:491
empty! at K:\judepot1115\packages\MathOptInterface\JdLrc\src\Bridges\bridge_optimizer.jl:378
default_copy_to at K:\judepot1115\packages\MathOptInterface\JdLrc\src\Utilities\copy.jl:385
copy_to at K:\judepot1115\packages\MathOptInterface\JdLrc\src\Bridges\bridge_optimizer.jl:443 [inlined]
optimize! at K:\judepot1115\packages\MathOptInterface\JdLrc\src\MathOptInterface.jl:121 [inlined]
optimize! at K:\judepot1115\packages\MathOptInterface\JdLrc\src\Utilities\cachingoptimizer.jl:370
jfptr_optimizeNOT._5034 at K:\judepot1115\compiled\v1.11\Gurobi\do9v6_pSn93.dll (unknown line)
#optimize!#104 at K:\judepot1115\packages\JuMP\RGIK3\src\optimizer_interface.jl:595
optimize! at K:\judepot1115\packages\JuMP\RGIK3\src\optimizer_interface.jl:546
jfptr_optimizeNOT._14010 at K:\judepot1115\compiled\v1.11\JuMP\DmXqY_pSn93.dll (unknown line)
macro expansion at .\REPL[10]:16 [inlined]
#1339#threadsfor_fun#119 at .\threadingconstructs.jl:253
#1339#threadsfor_fun at .\threadingconstructs.jl:220 [inlined]
#1 at .\threadingconstructs.jl:154
unknown function (ip: 000002aea736040b)
jl_apply at C:/workdir/src\julia.h:2157 [inlined]
start_task at C:/workdir/src\task.c:1202
Allocations: 46996111 (Pool: 46982645; Big: 13466); GC: 29

[WalterMadelim]

I revise my “crux” part code as

sub_j_vec = rand_j_vec[1+(ba-1)TH : min(J, (ba)TH)]
@info "before entering @threads" sub_j_vec
Threads.@threads for i = eachindex(sub_j_vec) # strictly This is NOT a loop
    j = sub_j_vec[i]
    @info "inside @threads" i j
end
error("#1 can you see this error??")

Then I observe this weird thing (why is there repetition of j = 3 ??)

julia> parallel_CG!(B, θ, β, μ, ν)
┌ Info: before entering @threads
│   sub_j_vec =
│    4-element Vector{Int64}:
│     1
│     4
│     2
└     3
┌ Info: inside @threads
│   i = 1
└   j = 2
┌ Info: inside @threads
│   i = 2
└   j = 3
┌ Info: inside @threads
│   i = 4
└   j = 3
┌ Info: inside @threads
│   i = 3
└   j = 3
ERROR: #1 can you see this error??

[WalterMadelim]

I give a short recap here:

There is one subtle thing we need to be wary of:
when writing ordinary for, we can write code like

for _ = 1:1
    for i = eachindex(v)
        j = v[i]
        println(j)
    end
    j = 9
end

This style is safe in that we can ensure the validity of j in println(j) by the just-in-time assignment j = v[i] preceding it.
This style is problematic when the inner for is equipped with the Threads.@threads prefix.

[WalterMadelim]

@WalterMadelim

I’m unsure if this style is correct.

Therefore it was incorrect. To fix, write

Threads.@threads for j = sub_j_vec
    local mj = inn[j]
    # do stuff with `mj`
end

[sgaure]

This is easily the most annoying feature of julia, incorrectness when a variable is inadvertently shared when annotating a loop with @threads.

Because of this, I have made a habit to make every variable I create inside a @threads or @spawn a local variable.

This unfortunate sharing is caught if you check your code for type instabilities. This is anyway very important for performance.

The reason the mj is shared is because it exists outside the loop, and the scoping rules says that the same variable should be used inside, unless it’s marked local.

function ec(N)
    s = 0.0
    Threads.@threads for i in 1:N
        s += 1/i
    end
    return s - log(N)
end

julia> @code_warntype ec(10^6)
MethodInstance for ec(::Int64)
  from ec(N) @ Main ~/foo.jl:1
Arguments
  #self#::Core.Const(Main.ec)
  N::Int64
Locals
  threadsfor_fun::var"#58#threadsfor_fun#11"{var"#58#threadsfor_fun#10#12"{UnitRange{Int64}}}
  s@_4::Core.Box
  threadsfor_fun#10::var"#58#threadsfor_fun#10#12"{UnitRange{Int64}}
  range::UnitRange{Int64}
Body::Any
1 ─       (s@_4 = Core.Box())
│   %2  = s@_4::Core.Box
...

The Core.Box under the Locals header shines in bright red when doing this (the colors are wrong here in Discourse). This signifies a type instability, i.e. the compiler does not know what type the variable has, a side effect of sharing variables in this fashion. It’s very bad for performance (and in this case, for correctness).

[WalterMadelim]

@odow

Gurobi environments are not threadsafe:

There is a statement therein, but I don’t know where it stems from.

The main reason to use JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) instead of JuMP.Model(Gurobi.Optimizer) is that I don’t want to see the redundant Gurobi logging lines, which floods my screen.

It appears that if I build a vector of models using the ordinary for, but solve them in parallel with @threads for, it can work normally and yield the expected result (I can try larger instances in the future to see if this still holds).