# Threads stop working when scaling computation

**URL:** <https://discourse.julialang.org/t/threads-stop-working-when-scaling-computation/114392>\
**Category:** Numerics\
**Created:** [May 17, 2024, 12:06pm UTC](https://discourse.julialang.org/t/threads-stop-working-when-scaling-computation/114392 "2024-05-17T12:06:39Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![SpuriousEigenstate](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spuriouseigenstate/32/208390_2.png) [@SpuriousEigenstate](https://discourse.julialang.org/u/SpuriousEigenstate)\
**Post date:** [May 17, 2024, 12:06pm UTC](https://discourse.julialang.org/t/threads-stop-working-when-scaling-computation/114392/1 "2024-05-17T12:06:39Z")

</div>

I have been stumped by this for a while now, so I thought I’d put it up here.

I developed a script that calculates some properties of a simulation “self-consistently”. I.e. the system is set up as a matrix, diagonalized, then some properties are calculated from the eigenvectors and plugged back in to the matrix. This process is then repeated until convergence (or some `n_max_iter`). I then want to do this for a bunch of different parameters, so I pass them to different threads using `Threads.@threads`.

The script works fine (fast & results look right) on my desktop, and for the smaller system sizes I used when developing the code. When I put the same, small system size (800x800 matrices), on our slurm cluster, (96 threads per node, I usually request 72) **it works fine**.

Increasing the system size to 1800x1800 matrices the script **still works** fine on my desktop (16 threads, 32G memory), however when I run it on the cluster and monitor with `htop`, initially all threads are running at 100% until suddenly after a minute or so, all threads go into `S` or sleep mode. They do not recover.

 ![Screenshot from 2024-05-16 10-10-40](https://global.discourse-cdn.com/julialang/original/3X/8/1/81b6600c626a7fd3738d3da6549358958045f1ff.png)

 ![Screenshot from 2024-05-16 10-12-07](https://global.discourse-cdn.com/julialang/original/3X/b/d/bd135a8135a9dd7e4501813723ad3c15f38f2a8f.png)

In fact, this already occurs at 20 threads, but **doesn’t** at 16 threads on the same machine. I’m stumped.

Things I tried include:

- Making sure it’s not a (too-little) memory issue. This is “easy” since the cluster has up to 256G memory available.
- Forcing garbage collection. I know my code allocates a lot (a consequence of running `LinearAlgebra.eigen()` hundreds of times

I made myself a sort of minimal example, where I diagonalize a bunch of random matrices in the same sort of loop structure as in my real code, and that doesn’t cause the sleeping to occur. So I conclude that the problem is somehow related to the way I update the matrix. For completeness, this is the exact code:

```julia
function update_H!(H, erray, N, E, uv, U, β, μ0, t, ϵ)
	erray[1] = 0.0
	erray[2] = 0.0
	for i in 1:N
		δi = 0+0im
		mi = μ0
		
		for (k, e) in enumerate(E)
			δi += U[i] * uv[2i-1, k] * conj(uv[2i, k]) * tanh(β * e / 2)
			mi -= U[i] * abs2(uv[2i-1, k]) * fermi_dirac(e, β)
		end

		if abs(H[2i-1, 2i] - δi) > erray[1]
			erray[1] = abs(H[2i-1, 2i] - δi)
		end
		
		if abs(H[2i-1, 2i-1] - mi) > erray[2]
			erray[2] = abs(H[2i-1, 2i-1] - (4t + ϵ[i] - mi))
		end

		H[2i-1, 2i] = δi
		H[2i, 2i-1] = conj(δi)

		H[2i-1, 2i-1] = 4t + ϵ[i] - mi
		H[2i, 2i] = -4t - ϵ[i] + mi
	end
	return nothing
end

function get_self_consistent(n, m, t, μ, μ0, ϵ, Δ, U, pbc_x, pbc_y, β; epsa=1e-4, N_iter=40)
    
    i = 0
    erra = 1
	erray = [0.0, 0.0]
	H = H_BdG(n, m, t, μ, ϵ, Δ, pbc_x=pbc_x, pbc_y=pbc_y)
	E = Array{Float64}(undef, n*m)
	uv = similar(H)

    while (i < N_iter) && (erra > epsa)

		E, uv = eigen(Hermitian(H))

		update_H!(H, erray, n*m, E, uv, U, β, μ0, t, ϵ)
		
        erra = max(erray...)

        i += 1
		print("i = ", i, " | ϵ(abs) = ", round.(erray,digits=5), "\n")
		flush(stdout)
    end

	Δ = calc_Δ(n*m, E, uv, U, β)
    μ = calc_μ(n*m, E, uv, U, β, µ0)
	# print("i = ", i, " | ϵ(abs) = ", round(erra,digits=5), "\n")
    return Δ, μ
end

```

So `H` is the matrix being diagonalized, `update_H!` updates `H` in place, and keeps track of the difference with the previous iteration (`erray`). There are also a bunch of parameters which I won’t explain, but left in for completeness. The `print `statement is only there for debugging purposes, and informed me that most workers don’t make it past the first few iterations of the self-consistent loop. The bug is also present without printing anything.

I’m sorry if this is a bit much, but I would appreciate any insight into what’s going on 🙂