# Julia code becomes slower on running on supercomputers and does not scale well when parallelizing with Base.Threads

**URL:** <https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439>\
**Category:** Julia at Scale\
**Tags:** fortran, parallel, linearalgebra, threads\
**Created:** [January 6, 2024, 2:11pm UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439 "2024-01-06T14:11:51Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 17, 2024, 1:52am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/22 "2024-01-17T01:52:02Z")

</div>

Could this be a related issue: [[OpenBLAS\_jll] Use new build compiled with GCC 11 by giordano · Pull Request #52928 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/52928) ?

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [January 17, 2024, 3:41am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/23 "2024-01-17T03:41:58Z")

</div>

> Could this be a related issue: [[OpenBLAS\_jll] Use new build compiled with GCC 11 by giordano · Pull Request #52928 · JuliaLang/julia · GitHub ](https://github.com/JuliaLang/julia/pull/52928) ?

Thank you for your reply. I am sorry but I am extremely new to such developments. Could you please paraphrase what you are referring to ?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 17, 2024, 7:12am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/24 "2024-01-17T07:12:53Z")

</div>

> [@Uranium238](#):
>
> Could you please paraphrase what you are referring to ?

New in the master branch of Julia is support for the [sapphire rapids](https://en.wikipedia.org/wiki/Sapphire_Rapids) CPU family from Intel, which covers the latest XEON CPUs that where released in 2023. You still did not report which exact CPU your super computer is using, but if it should be one of the latest XEON CPUs trying out Julia from master might be worth a try…

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [January 17, 2024, 7:27am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/25 "2024-01-17T07:27:57Z")

</div>

> [@ufechner7](#):
>
> You still did not report which exact CPU your super computer is using

I’m sorry for that, I did not realize the the importance of quoting the exact CPU. As per the website [here](http://www.hpc.iitkgp.ac.in/HPCF/systemConfig) we are using `Intel Xeon SKL G-6148`.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [January 17, 2024, 7:34am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/26 "2024-01-17T07:34:24Z")

</div>

I’m unsure if this is related, but there are known issues with garbage collection when running multithreaded code:

> <https://github.com/JuliaLang/julia/issues/40626>
>
> When returning values from tasks executed on another thread, those values are ke…pt alive even though it should be possible to collect them:
> 
> \`\`\`julia
> julia\> a = \[2,2\];
> julia\> finalizer(a) do \_
> Core.println("finalizing a")
> end
> 
> julia\> t = Threads.@spawn (println("returning a from thread $(Threads.threadid())"); a)
> Task (runnable) @0x00007fffb73e8e80
> returning a from thread 2
> 
> julia\> dump(t)
> Task
> next: Nothing nothing
> queue: Nothing nothing
> storage: Nothing nothing
> donenotify: Base.GenericCondition{Base.Threads.SpinLock}
> waitq: Base.InvasiveLinkedList{Task}
> head: Nothing nothing
> tail: Nothing nothing
> lock: Base.Threads.SpinLock
> owned: Int64 0
> result: Array{Int64}((2,)) \[2, 2\] # \<----------------- the array a is being kept alive by the task
> logstate: Nothing nothing
> code: #3 (function of type var"#3#4")
> \_state: UInt8 0x01
> sticky: Bool false
> \_isexception: Bool false
> 
> julia\> wait(t)
> julia\> a = nothing
> julia\> t = nothing
> julia\> GC.gc(true)
> 
> \# nothing is collected
> \`\`\`
> 
> Running under \`gdb\` reveals that the task object is being kept alive in the thread's local storage:
> 
> \`\`\`
> $ gdb --args julia -t2
> (gdb) run
> 
> julia\> code from above
> 
> (gdb) call jl\_(jl\_all\_tls\_states\[1\].current\_task)
> Task(next=nothing, queue=nothing, storage=nothing,
> donenotify=Base.GenericCondition{Base.Threads.SpinLock}(waitq=Base.InvasiveLinkedList{Task}(head=nothing, tail=nothing), lock=Base.Threads.SpinLock(owned=0)),
> result=Array{Int64, (2,)}\[2, 2\], # \<------------------- our array again
> logstate=nothing, code=Main.var"#3", \_state=0x01, sticky=false, \_isexception=false)
> \`\`\`
> 
> Running another task on the same thread replaces that sate and allows collection of the array:
> 
> \`\`\`julia
> julia\> # code from above
> 
> julia\> t = Threads.@spawn println("doing nothing from thread $(Threads.threadid())")
> Task (runnable) @0x00007fffb5b99e40
> doing nothing from thread 2
> 
> julia\> GC.gc(true)
> finalizing a
> \`\`\`
> 
> As observed in https://github.com/JuliaGPU/CUDA.jl/issues/866.

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [January 17, 2024, 9:05am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/27 "2024-01-17T09:05:03Z")

</div>

> [@jishnub](#):
>
> I’m unsure if this is related, but there are known issues with garbage collection when running multithreaded code:

The massive slowdown on using Xeon CPUs occurs regardless of whether I am running it on single or multi-threaded.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [January 17, 2024, 9:16am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/28 "2024-01-17T09:16:48Z")

</div>

Isn’t that the Fortran code is faster?

Personally I would focus in improving the performance of the Julia code. Do you have profiled it to see where is it taking most of the time? If it is is in the LA calls, then one can focus in improving that.

You have some slices and non-constant global variables there, and they might be a problem. But, without a profile, these are just guesses.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 17, 2024, 9:45am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/29 "2024-01-17T09:45:38Z")

</div>

> [@lmiq](#):
>
> Personally I would focus in improving the performance of the Julia code.

This is not the topic of this thread. The topic of this thread is why the same code is faster than Fortran on a laptop CPU and much slower than Fortran on a super computer CPU.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [January 17, 2024, 10:09am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/30 "2024-01-17T10:09:22Z")

</div>

But that can be only about the Fortran code and it’s compilers. There is no clear evidence that that has anything to do with Julia.

Also the thread is about why the Julia code does not scale well, and that is something do investigate optimizing the Julia code.

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [January 17, 2024, 10:38am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/31 "2024-01-17T10:38:31Z")

</div>

@lmiq As you rightly pointed out in your first reply, the primary focus should be to optimize the Julia code as much as possible. We have tried to do that as described in this thread [here](https://discourse.julialang.org/t/how-to-improve-performance-in-a-function-that-repeatedly-defines-and-multiplies-matrices/108153/59). We have significantly improved the performance of the serial Julia code. On observing that the Julia code is running faster than the serial Fortran code on i7 x86 architechture, we concluded that we have reached the end of “general programming optimization” , that we might have done without rewriting the matrix formalism that we used.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 17, 2024, 10:40am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/32 "2024-01-17T10:40:32Z")

</div>

> [@lmiq](#):
>
> There is no clear evidence that that has anything to do with Julia.

Would it be possible to provide this evidence?

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [January 17, 2024, 10:45am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/33 "2024-01-17T10:45:11Z")

</div>

The CPUs you are using are Skylake  
[Intel Xeon Gold 6148 Processor 27.5M Cache 2.40 GHz Product Specifications](https://ark.intel.com/content/www/us/en/ark/products/120489/intel-xeon-gold-6148-processor-27-5m-cache-2-40-ghz.html)

You still refer to massive slowdowns. Have you tried running htop during a job run as I propose?  
Honestly, some simple diagnostics may show things you don’t expect.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [January 17, 2024, 10:46am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/34 "2024-01-17T10:46:00Z")

</div>

Also you refer to MPI and threaded code. Are you trying to compare these?

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [January 17, 2024, 10:49am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/35 "2024-01-17T10:49:02Z")

</div>

My goal is to compare the parallel performance of the Fortran code with the Julia one. For Julia I am using `Base.Threads` to gain a thread level parallelization as it is the one with the least amount of overhead (as told in the Julia manual).

For Fortran we are using MPI to parallelize the code since as far as I know it does not have any native parallelization packages like Julia does. The comparison might not be “fair” in the truest sense but we are more concerned about getting the Julia code to run faster than the Fortran code in both serial and parallel (which we achieved successfully on i7 x86 but is failing with the Xeon chip).

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [January 17, 2024, 10:52am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/36 "2024-01-17T10:52:05Z")

</div>

Also your cluster has two sockets of 20 physical cores. The documentation says there are 80 cores per server, so hyperthreading is enabled.  
I would advise the following:  
You will not be able to have hyperthreading switched off in the BIOS (your sys admins will probably refuse to do this)  
You can disable all the even numbered cores at run time - again sys admins will probably refuse this

Could you try running on only odd numbered cores - I think hyperthreading may be having an influence here  
[GitHub - carstenbauer/ThreadPinning.jl: Readily pin Julia threads to CPU processors](https://github.com/carstenbauer/ThreadPinning.jl)

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [January 17, 2024, 10:58am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/37 "2024-01-17T10:58:55Z")

</div>

> [@johnh](#):
>
> Could you try running on only odd numbered cores - I think hyperthreading may be havign an influence here

I shall try my best to do this and provide the result (as stated above, I do not have direct access to the cluster in question). However, the slowdown occurs even when I am comparing the serial codes. Besides, when running for `20000` points, the Julia code comes of as about `5` minutes slower than the Fortran code (the Fortran code isn’t scaling perfectly either), so I suspect that this is the result of the accumulated serial code slowdown over multiple iterations.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [January 17, 2024, 11:30am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/38 "2024-01-17T11:30:18Z")

</div>

> [@Uranium238](#):
>
> For Fortran we are using MPI to parallelize the code since as far as I know it does not have any native parallelization packages like Julia does.

Please look at OpenMP  
[Using OpenMP with Fortran — Research Computing University of Colorado Boulder documentation (curc.readthedocs.io)](https://curc.readthedocs.io/en/latest/programming/OpenMP-Fortran.html)

[OpenMP | LLNL HPC Tutorials](https://hpc-tutorials.llnl.gov/openmp/)

You are comparing threaded code with MPI code.  
Fortran can used MPI or can use OpenMP For threaded code.  
Julia can use MPI or can use threads.  
I would say you chould compare MPI to MPI or threaded to threaded.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [January 17, 2024, 11:39am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/39 "2024-01-17T11:39:44Z")

</div>

I would really advise you to look at thread pinning here  
[ThreadPinning · ThreadPinning.jl (carstenbauer.github.io)](https://carstenbauer.github.io/ThreadPinning.jl/stable/)

[Autochecking BLAS Thread Settings · ThreadPinning.jl (carstenbauer.github.io)](https://carstenbauer.github.io/ThreadPinning.jl/stable/examples/ex_blas/)

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [January 17, 2024, 11:41am UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/40 "2024-01-17T11:41:23Z")

</div>

Also you say you do not have direct access to the cluster. Try the following:

ssh login to a login node  
Submit your job  
use qstat to see which compute node is in use  
ssh to that compute node  
Often HPC clusters stop you logging in to compute nodes - but you CAN do this when you are running a job on a node

Also look at running an interactive job - which batch queue system are you using?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [January 17, 2024, 12:16pm UTC](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439/41 "2024-01-17T12:16:33Z")

</div>

After profiling the code (with `@profview` in VSCode), I can see that `vcat` operations were taking a lot of time, in these lines:

```julia
# X = [X11 X12
# X21 X22]
and
# Y = [Y1
# Y2]

```

By allocating these intermediate `X` and `Y` outside the loop and updating their values instead of recreating the arrays, the number of allocations is cut in half:  
Before:

```julia
leandro@m3g:~/Downloads/uranium% julia -t 1 --project pal-rv6.jl 
Number of threads = 1
 31.559063 seconds (9.76 M allocations: 42.744 GiB, 1.54% gc time, 14.28% compilation time)
GKS: Possible loss of precision in routine SET_WINDOW
The rate constant is 2.1751344328561788e9

```

After:

```julia
leandro@m3g:~/Downloads/uranium% julia -t 1 --project pal-rv6_2.jl 
Number of threads = 1
 29.952532 seconds (9.86 M allocations: 21.572 GiB, 1.35% gc time, 15.68% compilation time)
GKS: Possible loss of precision in routine SET_WINDOW
The rate constant is 2.1751344328561788e9

```

Next I note that a line inside GHT which performs matrix multiplications is taking most of the time, and calling `similar` internally. So I preallocate the result of the matrix multiplication to use `mul!(WinX, Win, X)` instead, and get:

```julia
leandro@m3g:~/Downloads/uranium% julia -t 1 --project pal-rv6_2.jl 
Number of threads = 1
 30.032092 seconds (9.74 M allocations: 847.634 MiB, 0.86% gc time, 16.09% compilation time)
GKS: Possible loss of precision in routine SET_WINDOW
The rate constant is 2.1751344328561788e9

```

note that allocations are now in 874Mib. And mostly because of compilation. In a second run I get:

```julia
julia> @time calc(t, b, Delta, sysize, omegaf, omegai, J, D, P; nchunks=1)
 25.172865 seconds (311.39 k allocations: 228.236 MiB, 0.05% gc time)
The rate constant is 2.1751344328561788e9

```

and with 8 threads:

```julia
@time calc(t, b, Delta, sysize, omegaf, omegai, J, D, P; nchunks=8)
  4.791953 seconds (311.52 k allocations: 230.700 MiB)
The rate constant is 2.1751344328561788e9

```

I think this a better starting point to check the scalability in the cluster. Excessive allocations can really be a problem for that.

Here is the modified code I’m running. It has some other small differences which didn´t make a difference.

> **code**
>
> ```julia
> using MKL
> using DelimitedFiles
> using LinearAlgebra
> using LaTeXStrings
> using TimerOutputs
> using Plots
> using Plots.PlotMeasures
> using FFTW
> using Base.Threads
> using SharedArrays
> using ChunkSplitters
> BLAS.set_num_threads(1)
> 
> println("Number of threads = $(Threads.nthreads())")
> open("THREAD-PARALLEL-OUTPUT-RV6.TXT", "a") do f
> write(f, "\n\n\n====================================================\n\n\n")
> end
> 
> function GFC(
> t::Float64,
> b::Float64,
> Δ::Float64,
> sysize::Int64,
> Ωf::Diagonal{Float64},
> Ωi::Diagonal{Float64},
> J::Matrix{Float64},
> D::Matrix{Float64},
> ρ::Diagonal{Float64},
> Si::Diagonal{ComplexF64},
> Sf::Diagonal{ComplexF64},
> Bi::Diagonal{ComplexF64},
> Bf::Diagonal{ComplexF64}
> )
> t = t * 4.1341373335e16 # converting time to atomic units
> tp = -im * b - t
> for i = 1:sysize
> Sf[i, i] = Ωf[i, i] / sin(Ωf[i, i] * t)
> Si[i, i] = Ωi[i, i] / sin(Ωi[i, i] * tp)
> Bf[i, i] = Ωf[i, i] / tan(Ωf[i, i] * t)
> Bi[i, i] = Ωi[i, i] / tan(Ωi[i, i] * tp)
> end
> BB = (Bf + J' * Bi * J)
> AA = (Sf + J' * Si * J)
> 
> W = [BB -AA
> -AA BB] #We wont calculate the 2Nx2N determinant directly.
> 
> # Wapp = BB*(BB- AA*(BB^-1)*AA )
> Wapp = BB * (BB - AA * (BB \ AA))
> U = Bi - Si
> V = [J' * U * D
> J' * U * D]
> # F1 = sqrt(det(Si*Sf*ρ^2*Wapp^-1))
> # F1 = sqrt(det(Si*Sf*ρ^2)/det(Wapp))
> F1 = sqrt(det(Si * Sf * (Wapp \ ρ^2)))
> F2 = (-(im / 2) * (transpose(V) * (W \ V))) + ((im) * (transpose(D) * U * D))
> g = F1 * exp(F2)
> Δ = Δ * 0.0367493 # in atomic units
> g = g * exp(im * Δ * t)
> η = 2 * 4.55633e-6 # 10 cm^-1 to hatree energy
> #g = g * exp(-η * abs(t)) # Adding the damping factor Lorrentzian
> g = g * exp(-η * t^2) # Adding the damping factor Gaussian
> # println("t = $t was evaluated by process $(myid())")
> # println("GFC at t = $t is $g")
> # println("F1 at t = $t = $F1")
> return g[1]
> end
> 
> @views function GHT(
> t::Float64,
> b::Float64,
> sysize::Int64,
> Ωf::Diagonal{Float64},
> Ωi::Diagonal{Float64},
> J::Matrix{Float64},
> D::Matrix{Float64},
> Si::Diagonal{ComplexF64},
> Sf::Diagonal{ComplexF64},
> Bi::Diagonal{ComplexF64},
> Bf::Diagonal{ComplexF64},
> X11::Matrix{ComplexF64},
> X12::Matrix{ComplexF64},
> X21::Matrix{ComplexF64},
> X22::Matrix{ComplexF64},
> Y1::Matrix{ComplexF64},
> Y2::Matrix{ComplexF64},
> gfcpart::ComplexF64,
> T::Matrix{Float64}
> )
> 
> ################################## Preliminary Part , Common for both Frank Condon and Herzberg Teller ##############################
> t = t * 4.1341373335e16 # converting time to atomic units
> tp = -im * b - t
> for i = 1:sysize
> Sf[i, i] = Ωf[i, i] / sin(Ωf[i, i] * t)
> Si[i, i] = Ωi[i, i] / sin(Ωi[i, i] * tp)
> Bf[i, i] = Ωf[i, i] / tan(Ωf[i, i] * t)
> Bi[i, i] = Ωi[i, i] / tan(Ωi[i, i] * tp)
> end
> W = [(Bf+J'*Bi*J) -(Sf + J' * Si * J)
> -(Sf + J' * Si * J) (Bf+J'*Bi*J)]
> U = Bi - Si
> V = [J' * U * D
> J' * U * D]
> 
> ########## Defining the X and Y Matrices ##################################
> 
> ################# Final Calculation and returning value ##################
> s = zero(ComplexF64)
> MJSJ = J' * Si * J
> MJBJ = J' * Bi * J
> MJUD = J' * U * D
> Win = inv(W)
> WinV = Win * V
> X = zeros(ComplexF64, 2 * sysize, 2 * sysize)
> Y = zeros(ComplexF64, 2 * sysize, 1)
> WinX = similar(Win)
> for m in 1:sysize
> for mp in 1:sysize
> X11[m, :] .= (MJSJ)[mp, :] .* (-Bf[m, m])
> X12[m, :] .= (MJBJ)[mp, :] .* (Bf[m, m])
> X21[m, :] .= (MJSJ)[mp, :] .* (Sf[m, m])
> X22[m, :] .= (MJBJ)[mp, :] .* (-Sf[m, m])
> # X = [X11 X12
> # X21 X22]
> X[1:sysize, 1:sysize] .= X11
> X[1:sysize, sysize+1:2*sysize] .= X12
> X[sysize+1:2*sysize, 1:sysize] .= X21
> X[sysize+1:2*sysize, sysize+1:2*sysize] .= X22
> 
> Y1[m, :] .= Bf[m, m] * (MJUD)[mp, :]
> Y2[m, :] .= -Sf[m, m] * (MJUD)[mp, :]
> Y[1:sysize, 1] .= Y1
> Y[sysize+1:2*sysize, 1] .= Y2
> # Y = [Y1
> # Y2]
> # g = im*tr(W\X)+ (transpose(W\V)*X*(W\V))[1] - (transpose(Y)*(W\V))[1]
> # g = im*tr(Win*X)+ (transpose(Win*V)*X*(Win*V))[1] - (transpose(Y)*(Win*V))[1]
> g = im * tr(mul!(WinX, Win, X)) + (transpose(WinV)*X*(WinV))[1] - (transpose(Y)*(WinV))[1]
> s = s + g * gfcpart * T[m, mp]
> end
> X11[m, 1:sysize] .= 0
> X12[m, 1:sysize] .= 0
> X21[m, 1:sysize] .= 0
> X22[m, 1:sysize] .= 0
> Y1[m, 1] = 0
> Y2[m, 1] = 0
> end
> return s
> end
> 
> const Delta = 2.02 # in eV
> const dd = 11 # Change this to modify number of points
> const Temp = 300
> kT_temp = Temp * 1.380649e-23 # Boltzmann constant times temperature in Joules
> const kT = kT_temp * 2.29371227840e17 # in atomic units
> const b = 1 / kT
> const omegaf = Diagonal(vec(readdlm("wf.txt", '\t', Float64, '\n'))) * 4.55633e-6 # Final state frequencies in atomic units
> const omegai = Diagonal(vec(readdlm("wi.txt", '\t', Float64, '\n'))) * 4.55633e-6 # Initial state frequencies in atomic units
> const sysize = size(omegai)[1]
> const P = @. 2 * sinh(b * omegai / 2)
> T_temp = readdlm("nac.txt", '\t', Float64, '\n') # NAC matrix
> const T = T_temp * T_temp'
> const D = readdlm("d.txt", Float64) # Displacement Matrix
> const J = readdlm("j.txt", Float64) # Duchinsky Matrix
> const t = collect(range(-5e-12, stop=5e-12, length=dd))
> t[div(dd + 1, 2)] = 10e-25
> 
> function calc(t, b, Δ, sysize, Ωf, Ωi, J, D, P; nchunks=Threads.nthreads())
> # @sync tells Julia to wait for any tasks started inside it
> open("THREAD-PARALLEL-OUTPUT-RV6.TXT", "a") do f
> write(f, "STARTING LOOP WITH NUMBER OF THREADS = $(Threads.nthreads()) \n====================================================\n")
> end
> x = zeros(ComplexF64, dd)
> @time @sync for (xrange, ichunk) in chunks(range(1, dd), nchunks)
> # xrange are the indices this chunk contains, ichunk is just the index of the chunk, which we don't need in this scenario
> Threads.@spawn begin # this line tells Julia to start a new thread with the contents of the begin block
> Sf = Diagonal(zeros(ComplexF64, sysize))
> Si = Diagonal(zeros(ComplexF64, sysize))
> Bf = Diagonal(zeros(ComplexF64, sysize))
> Bi = Diagonal(zeros(ComplexF64, sysize))
> X11 = zeros(ComplexF64, sysize, sysize)
> X12 = zeros(ComplexF64, sysize, sysize)
> X21 = zeros(ComplexF64, sysize, sysize)
> X22 = zeros(ComplexF64, sysize, sysize)
> Y1 = zeros(ComplexF64, sysize, 1)
> Y2 = zeros(ComplexF64, sysize, 1)
> for i in xrange
> gfcpart = GFC(t[i], b, Δ, sysize, Ωf, Ωi, J, D, P, Si, Sf, Bi, Bf)
> gtotal = GHT(t[i], b, sysize, Ωf, Ωi, J, D, Si, Sf, Bi, Bf, X11, X12, X21, X22, Y1, Y2, gfcpart, T)
> x[i] = gtotal
> open("THREAD-PARALLEL-OUTPUT-RV6.TXT", "a") do f
> write(f, "ITERATION $i COMPLETED BY THREAD $(Threads.threadid()) \n")
> end
> end # loop
> end # thread work load
> end # loop over chunks, Julia will wait until all task are completed
> open("THREAD-PARALLEL-OUTPUT-RV6.TXT", "a") do f
> write(f, "LOOP COMPLETE \n====================================================")
> end
> yr = real.(x)
> yi = imag.(x)
> p = plot(t, yr, linewidth=2, color="blue", title="Δ = $Δ", label="Real part", grid=false)
> p = plot!(t, yi, linewidth=2, color="red", label="Imaginary part", size=(1000, 800), xlims=(-30e-15, 30e-15))
> savefig(p, "gt.svg")
> begin # This part does FFTW and finds final rate constant
> 
> fs = abs(t[1] - t[2])^-1
> signal = yr
> p = plan_fft(signal)
> F = fftshift(p * signal)
> freqs = fftshift(fftfreq(length(t), fs)) # X axis has frequencies in s^-1
> #Covert freqs to eV
> freqs = freqs * 4.13558e-15
> F = F * 6.57e15
> time_domain = plot(t, signal, title=L"G(t)", legend=false, color="orange", xlabel=L"\mathrm{Time}\ (s)", ylabel=L"\mathrm{A.u}")
> freq_domain = plot(freqs, abs.(F), title=L"\mathrm{Fourier\ transformed\ g(\bar{\nu}))}", legend=false, xlabel=L"\mathrm{Energy\ (eV)}", ylabel=L"\mathrm{k_{IC}\ (s^{-1})}", xlims=(0, 4))
> p = plot(time_domain, freq_domain, layout=(2, 1), size=(1000, 800), right_margin=10mm)
> savefig(p, "gt-fft.svg")
> ee = Δ# * 8065.73 #Converting eV to cm^-1
> index = findmin(abs.(freqs .- ee))[2] # Find the index of freq closest to ee
> freqs[index]
> println("The rate constant is ", abs.(F)[index])
> rc = abs.(F)[index]
> end
> begin # This part writes the output to a file
> open("THREAD-PARALLEL-OUTPUT-RV6.TXT", "a") do f
> write(f, "-------------------------\n")
> write(f, "NUMBER OF CORES = $(Threads.nthreads())\n")
> write(f, "The value of Δ is $Δ\n")
> write(f, "Number of time points = $dd\n")
> write(f, "The rate constant is $rc\n")
> write(f, "-------------------------\n")
> write(f, "The values of t,Re(g(t)),Im(g(t)) are \n\n")
> for i in 1:dd
> write(f, "$(t[i]) $(yr[i]) $(yi[i])\n")
> end
> end
> end
> 
> end
> 
> #timeneed = @elapsed calc(t, b, Delta, sysize, omegaf, omegai, J, D, P; nchunks=Threads.nthreads())
> #open("THREAD-PARALLEL-OUTPUT-RV6.TXT", "a") do f
> # write(f, "\n\n\n===============================================================\nThe time elapsed for main loop [INCLUDING GC + C + RC]= $timeneed seconds \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n")
> #end
> 
> ```

The performance of the final code improved a little, but the matrix multiplications is still the bottleneck, and I don’t think there’s nothing that you can improve much further there in the serial version. But these better allocation profiles can improve the parallel version, effectively.

A different question is why the Fortran code in your machine is so slow. I don’t have an answer to that, but from your original post it seemed that the Fortran code compiled in the cluster had a roughly similar performance than the Julia one (serial version), which is consistent with the matrix multiplication being the bottleneck.

edit: If you go to larger and larger problems, of course those repetitive allocations can become a problem, and may explain why your Julia code is getting slower when you increase the size.

[Previous page](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439.md?page=1)

[Next page](https://discourse.julialang.org/t/julia-code-becomes-slower-on-running-on-supercomputers-and-does-not-scale-well-when-parallelizing-with-base-threads/108439.md?page=3)