Performance improvemnts

 - got rid of all allocations in calc_chigamma_impr
 - imrpoved performance of threaded version
 - loops reordered

Author: Atomtomate

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.5"
 manifest_format = "2.0"
-project_hash = "9ea45be61937a755ee6757eec725af4882f38692"
+project_hash = "17cb083765589e55555c6f9f8a12f6c841814a84"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -126,11 +126,11 @@ version = "1.11.0"
 
 [[deps.BSE_Asymptotics]]
 deps = ["ColorSchemes", "Documenter", "JLD2", "LinearAlgebra", "OffsetArrays", "Test", "TimerOutputs"]
-git-tree-sha1 = "7bda0c84f387e2ef73e0523b25cb144346389396"
+git-tree-sha1 = "df23f5ec4ec8bc5456e1734ca539954237ffffa0"
 repo-rev = "master"
 repo-url = "https://github.com/Atomtomate/BSE_Asymptotics.jl"
 uuid = "646ec391-26e4-4808-a965-02314d3a2269"
-version = "3.0.2"
+version = "3.3.0"
 
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
@@ -309,9 +309,9 @@ version = "0.9.4"
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
-git-tree-sha1 = "9d733459cea04dcf1c41522ec25c31576387be8a"
+git-tree-sha1 = "b7af952d4701252dc45d3b0025693e9cb4dedcd8"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.10.1"
+version = "1.10.2"
 
 [[deps.DocumenterTools]]
 deps = ["AbstractTrees", "Base64", "DocStringExtensions", "Documenter", "FileWatching", "Git", "Gumbo", "LibGit2", "OpenSSH_jll", "Sass"]
@@ -749,9 +749,9 @@ version = "0.12.3"
 
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
-git-tree-sha1 = "44f6c1f38f77cafef9450ff93946c53bd9ca16ff"
+git-tree-sha1 = "7d2f8f21da5db6a806faf7b9b292296da42b2810"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.8.2"
+version = "2.8.3"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"]

Project.toml

@@ -2,7 +2,7 @@ name = "LadderDGA"
 uuid = "78e40beb-bf89-4c0e-9d2b-bee278912f2b"
 authors = ["Julian Stobbe <Atomtomate@gmx.de> and contributors"]
 
-version = "1.6.0"
+version = "1.7.0"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
@@ -36,7 +36,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
-BSE_Asymptotics = "3.0.2"
+BSE_Asymptotics = "3.3.0"
 DataFrames = "1.7.0"
 Dispersions = "1.0.1"
 EquivalenceClassesConstructor = "0.5.5"

src/BSETools.jl

@@ -214,14 +214,15 @@ This method solves the following equation:
 ``
 """
 function calc_χγ(type::Symbol, h::Union{lDΓAHelper,AlDΓAHelper}, χ₀::χ₀T; verbose=true, ω_symmetric::Bool=false, use_threads=false)
-    calc_χγ(type, getfield(h, Symbol("Γ_$(type)")), χ₀, h.kG, h.mP, h.sP, verbose=verbose, ω_symmetric=ω_symmetric, use_threads=use_threads)
+    calc_χγ(type, getfield(h, Symbol("Γ_$(type)")), χ₀, h.kG, h.mP, h.sP; verbose=verbose, ω_symmetric=ω_symmetric, use_threads=use_threads)
 end
 
-function calc_χγ(type::Symbol, Γr::ΓT, χ₀::χ₀T, kG::KGrid, mP::ModelParameters, sP::SimulationParameters; verbose=true, ω_symmetric::Bool=false, use_threads=false)
+function calc_χγ(type::Symbol, Γr::ΓT, χ₀::χ₀T, kG::KGrid, mP::ModelParameters, sP::SimulationParameters;
+         verbose=true, ω_symmetric::Bool=false, use_threads=false)
     s = type == :d ? -1 : 1
     !(type in (:d, :m)) && error("Unkown type $type")
 
-    NT = Threads.nthreads()
+    
     Nν = 2 * sP.n_iν
     Nq = length(kG.kMult)
     Nω = size(χ₀.data, χ₀.axis_types[:ω])
@@ -232,31 +233,36 @@ function calc_χγ(type::Symbol, Γr::ΓT, χ₀::χ₀T, kG::KGrid, mP::ModelPa
     qi_range = 1:Nq
     χ_ω = Vector{Float64}(undef, Nω)
 
-
     if use_threads
+        bthreads = BLAS.get_num_threads()
+        BLAS.set_num_threads(1)
+        NT = Threads.nthreads()
         χννpω = [Matrix{ComplexF64}(undef, Nν, Nν) for ti in 1:NT]
         ipiv = [Vector{Int}(undef, Nν) for ti in 1:NT]
-        work = [_gen_inv_work_arr(χννpω[1], ipiv[1]) for ti in 1:NT]
+        work = [_gen_inv_work_arr(χννpω[Threads.threadid()], ipiv[Threads.threadid()]) for ti in 1:NT]
         λ_cache = [Vector{ComplexF64}(undef, Nν) for ti in 1:NT]
+        diag_cache = [similar(sP.χ_helper.diag_asym_buffer) for ti in 1:NT]
         Threads.@threads for qi in qi_range
             for ωm in ωm_range
                 ωi =  ωm + sP.n_iω + 1
                 invert_BSE!(χ, γ, χννpω[Threads.threadid()], ipiv[Threads.threadid()], work[Threads.threadid()], 
-                        sP.χ_helper, λ_cache[Threads.threadid()], type, ω_symmetric, s, Γr, χ₀.data, χ₀.asym, 
+                        sP.χ_helper, λ_cache[Threads.threadid()], diag_cache[Threads.threadid()], type, ω_symmetric, s, Γr, χ₀.data, χ₀.asym, 
                         χ₀.ν_shell_size, qi, sP.n_iω, ωm, ωi, mP.U, mP.β)
             end
         end
+        BLAS.set_num_threads(bthreads)
     else
         χννpω = Matrix{ComplexF64}(undef, Nν, Nν)
         ipiv = Vector{Int}(undef, Nν)
         work = _gen_inv_work_arr(χννpω, ipiv)
         λ_cache = Vector{ComplexF64}(undef, Nν)
 
-        for qi in qi_range
-            for ωm in ωm_range
-                ωi =  ωm + sP.n_iω + 1
-                invert_BSE!(χ, γ, χννpω, ipiv, work, sP.χ_helper, λ_cache, type, ω_symmetric, s, Γr, χ₀.data, χ₀.asym,
-                            χ₀.ν_shell_size, qi, sP.n_iω, ωm, ωi, mP.U, mP.β)
+        for ωm in ωm_range
+            ωi =  ωm + sP.n_iω + 1
+            for qi in qi_range
+                invert_BSE!(χ, γ, χννpω, ipiv, work, 
+                    sP.χ_helper, λ_cache, sP.χ_helper.diag_asym_buffer, type, ω_symmetric, s, Γr, χ₀.data, χ₀.asym, 
+                    χ₀.ν_shell_size, qi, sP.n_iω, ωm, ωi, mP.U, mP.β)
             end
         end
     end
@@ -272,23 +278,25 @@ function calc_χγ(type::Symbol, Γr::ΓT, χ₀::χ₀T, kG::KGrid, mP::ModelPa
 end
 
 function invert_BSE!(χ::AbstractArray{Float64,2}, γ::AbstractArray{ComplexF64,3}, χννpω::AbstractArray{ComplexF64,2}, 
-                    ipiv::Vector{Int}, work::Vector{ComplexF64}, χ_helper, λ_cache, type::Symbol, ω_symmetric::Bool, s::Int, 
+                    ipiv::Vector{Int}, work::Vector{ComplexF64}, χ_helper::BSE_Asym_Helper, λ_cache::Vector{ComplexF64}, diag_cache::Vector{ComplexF64}, 
+                    type::Symbol, ω_symmetric::Bool, s::Int, 
                     Γr::AbstractArray{ComplexF64,3}, χ₀_data::AbstractArray{ComplexF64,3}, χ₀_asym::AbstractArray{ComplexF64,2},
                     ν_shell_size::Int, qi::Int, ωm_max::Int, ωm::Int, ωi::Int, U::Float64, β::Float64)
-    χννpω[:, :] = deepcopy(Γr[:, :, ωi])
+    
+    copy!(χννpω, view(Γr,:,:, ωi))
     for l in axes(χννpω, 1)
-        χννpω[l, l] += 1.0 / χ₀_data[qi, ν_shell_size+l, ωi]
+        @inbounds χννpω[l, l] += 1.0 / χ₀_data[qi, ν_shell_size+l, ωi]
     end
-    
     inv!(χννpω, ipiv, work)
     if typeof(χ_helper) <: BSE_Asym_Helpers
         χ[qi, ωi] = real(
-            calc_χλ_impr!(λ_cache, type, ωm, χννpω, 
-                          view(χ₀_data, qi, :, ωi), U, β,
-                          χ₀_asym[qi, ωi],  χ_helper,
-            ),
-        )
-        γ[qi, :, ωi] = (1 .- s * λ_cache) ./ (1 .+ s * U .* χ[qi, ωi])
+            calc_χλ_impr!(λ_cache, diag_cache, type, qi, ωi, ωm, χννpω, 
+                          χ₀_data, U, β,
+                          χ₀_asym[qi, ωi], χ_helper,
+            ))
+        for νi in axes(γ,2)
+            @inbounds γ[qi, νi, ωi] = (1 - s * λ_cache[νi]) / (1 + s * U * χ[qi, ωi])
+        end
     else
         if typeof(χ_helper) === BSE_SC_Helper
             improve_χ!(type, ωi, view(χννpω, :, :, ωi), view(χ₀, qi, :, ωi), U, β, χ_helper)
@@ -306,7 +314,6 @@ function invert_BSE!(χ::AbstractArray{Float64,2}, γ::AbstractArray{ComplexF64,
 end
 
 
-
 """
     calc_gen_χ(Γr::ΓT, χ₀::χ₀T, kG::KGrid)