# Alternative permutation/reducing for loops

**URL:** <https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554>\
**Category:** Performance\
**Tags:** question, package, linearalgebra\
**Created:** [September 5, 2023, 9:29pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554 "2023-09-05T21:29:02Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Leticia-maria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leticia-maria/32/30981_2.png) [@Leticia-maria](https://discourse.julialang.org/u/Leticia-maria)\
**Post date:** [September 5, 2023, 9:29pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/1 "2023-09-05T21:29:02Z")

</div>

Hello, I am writing this function where I need to combine all terms in this structure (called Basis):

```julia
struct GaussianBasisSet <: AbstractBasisSet
    R
    α::Matrix{Float64}
    d::Matrix{Float64}
    ℓ::Int
    m::Int
    n::Int
end

```

Example:

```julia
13-element Vector{BasisSets.GaussianBasisSet}:
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [5484.67166 825.234946 … 16.8975704 5.79963534], [0.00183107443 0.0139501722 … 0.470192898 0.358520853], 0, 0, 0)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [15.53961625 3.599933586 1.01376175], [-0.1107775495 -0.1480262627 1.130767015], 0, 0, 0)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [15.53961625 3.599933586 1.01376175], [0.07087426823 0.3397528391 0.7271585773], 1, 0, 0)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [15.53961625 3.599933586 1.01376175], [0.07087426823 0.3397528391 0.7271585773], 0, 1, 0)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [15.53961625 3.599933586 1.01376175], [0.07087426823 0.3397528391 0.7271585773], 0, 0, 1)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [0.2700058226;;], [1.0;;], 0, 0, 0)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [0.2700058226;;], [1.0;;], 1, 0, 0)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [0.2700058226;;], [1.0;;], 0, 1, 0)
 BasisSets.GaussianBasisSet([0.0 -0.007156 0.965491], [0.2700058226;;], [1.0;;], 0, 0, 1)
 BasisSets.GaussianBasisSet([0.0 0.001486 -0.003471], [18.73113696 2.825394365 0.6401216923], [0.03349460434 0.2347269535 0.8137573261], 0, 0, 0)
 BasisSets.GaussianBasisSet([0.0 0.001486 -0.003471], [0.1612777588;;], [1.0;;], 0, 0, 0)
 BasisSets.GaussianBasisSet([0.0 0.931026 1.207929], [18.73113696 2.825394365 0.6401216923], [0.03349460434 0.2347269535 0.8137573261], 0, 0, 0)
 BasisSets.GaussianBasisSet([0.0 0.931026 1.207929], [0.1612777588;;], [1.0;;], 0, 0, 0)

```

```julia
function overlap(basis)
    K = length(basis)
    S = zeros(K, K)

    for (i, basisᵢ) in enumerate(basis)
        for (j, basisⱼ) in enumerate(basis)
            for (αᵢ, dᵢ) in zip(basisᵢ.α, basisᵢ.d)
                for (αⱼ, dⱼ) in zip(basisⱼ.α, basisⱼ.d)

                    Rᵢ = basisᵢ.R
                    Rⱼ = basisⱼ.R

                    ℓᵢ, mᵢ, nᵢ = basisᵢ.ℓ, basisᵢ.m, basisᵢ.n
                    ℓⱼ, mⱼ, nⱼ = basisⱼ.ℓ, basisⱼ.m, basisⱼ.n

                    S[i, j] += (
                        exp(-αᵢ * αⱼ * distance(Rᵢ, Rⱼ) / (αᵢ + αⱼ)) *
                        normalization(αᵢ, ℓᵢ, mᵢ, nᵢ) *
                        normalization(αⱼ, ℓⱼ, mⱼ, nⱼ) *
                        dᵢ *
                        dⱼ *
                        Sxyz(Rᵢ, Rⱼ, αᵢ, αⱼ, ℓᵢ, ℓⱼ, mᵢ, mⱼ, nᵢ, nⱼ)
                    )
                end
            end
        end
    end

    return S
end

```

However, this is expensive and I would like to reduce the number of for loops to make it more performant. Any tips?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 5, 2023, 9:42pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/2 "2023-09-05T21:42:07Z")

</div>

First things to do are to put a type on your `R` parameter which will make this more type stable. The next is to move your computation of `distance(Rᵢ, Rⱼ)` out two loops.

---

<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:** [September 5, 2023, 9:43pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/3 "2023-09-05T21:43:02Z")

</div>

And use a `SVector{3,Float64}` for `R`

---

<div class="post-metadata">

**Author:** ![Leticia-maria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leticia-maria/32/30981_2.png) [@Leticia-maria](https://discourse.julialang.org/u/Leticia-maria)\
**Post date:** [September 5, 2023, 9:46pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/4 "2023-09-05T21:46:22Z")

</div>

Does `SVector()` allocate less memory?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 5, 2023, 9:53pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/5 "2023-09-05T21:53:40Z")

</div>

No (but also yes). You still have to allocate space for the 3 numbers, but in a type stable context, you can often inline that allocation into the parent object which will reduce the amount of space the object takes up.

---

<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:** [September 5, 2023, 10:46pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/6 "2023-09-05T22:46:56Z")

</div>

Also arithmetic operations with SVectors return SVectors, and these resulting vectors then do not allocate (on the heap, which is what matters). Thus it may be possible that you get speedups depending on what you have inside `Sxyz`.

---

<div class="post-metadata">

**Author:** ![Leticia-maria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leticia-maria/32/30981_2.png) [@Leticia-maria](https://discourse.julialang.org/u/Leticia-maria)\
**Post date:** [September 7, 2023, 2:22pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/7 "2023-09-07T14:22:19Z")

</div>

I have moved the `distance(R_i, R_j)` out of the loops, but nothing has changed in terms of allocations or computing time.

```julia
julia> @benchmark overlap_2(basis, hydrogen)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 254.334 μs … 26.683 ms ┊ GC (min … max): 0.00% … 98.85%
 Time (median): 276.459 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 300.024 μs ± 558.434 μs ┊ GC (mean ± σ): 6.98% ± 3.78%

       ▂▅▇▆▃ ▂▄▆▇█▆▅▅▂▁
  ▂▂▃▅█████████████████████▇█▆▆▅▆▄▄▄▃▃▄▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂ ▄
  254 μs Histogram: frequency by time 340 μs <

 Memory estimate: 317.41 KiB, allocs estimate: 13672.

julia> @benchmark overlap(basis, hydrogen)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 252.375 μs … 31.246 ms ┊ GC (min … max): 0.00% … 98.93%
 Time (median): 278.708 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 300.581 μs ± 616.198 μs ┊ GC (mean ± σ): 6.59% ± 3.26%

     ▂ ▂▄▅█▇▆▄▃▃▃▂▂▁▁▁
  ▁▂▇█▇▄▅▇███████████████▇▆▅▄▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  252 μs Histogram: frequency by time 366 μs <

 Memory estimate: 317.41 KiB, allocs estimate: 13672.

```

Any other suggestions?

---

<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:** [September 7, 2023, 3:50pm UTC](https://discourse.julialang.org/t/alternative-permutation-reducing-for-loops/103554/8 "2023-09-07T15:50:07Z")

</div>

Did you fix the type of `R` in the struct? It will be easier to help if you provide the updated code and a running example.