# What's a good and easy way to parallelize this?

**URL:** <https://discourse.julialang.org/t/whats-a-good-and-easy-way-to-parallelize-this/74893>\
**Category:** Performance\
**Created:** [January 20, 2022, 12:55am UTC](https://discourse.julialang.org/t/whats-a-good-and-easy-way-to-parallelize-this/74893 "2022-01-20T00:55:00Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [January 20, 2022, 12:55am UTC](https://discourse.julialang.org/t/whats-a-good-and-easy-way-to-parallelize-this/74893/1 "2022-01-20T00:55:00Z")

</div>

Ultimately, I want to speed up `fullsys!`, and as you can see there are functions that fill arrays. I thought those filling operations could be parallelized. How can I get started with that?

```julia
const szAge = 64 - 19 + 1
const szE = 2
const β = 0.92
const δ = 0.10

maxT(am, af; Z=realage(szAge)) = Z - max(am, af) + 1

realage(x) = x + 18
realageinv(x) = x - 18

function insideproduct_num(am, em, af, ef; μ0m=μ0m, μ0f=μ0f)
    exponent = 0.5*(β*(1 - δ))
    power = 1.0
    p = one(promote_type(eltype(μ0m), eltype(μ0f)))
    for k in 0:maxT(realage(am), realage(af))-1
        p *= abs(μ0m[am + k, em] * μ0f[af + k, ef])^power
        power *= exponent
    end
    return p
end

function insideproduct_den(am, em, af, ef; Sm=Sm, Sf=Sf)
    exponent = 0.5*(β*(1 - δ))
    power = 1.0
    p = one(promote_type(eltype(Sm), eltype(Sf)))
    for k in 0:maxT(realage(am), realage(af))-1
        p *= abs(Sm[am + k, em] * Sf[af + k, ef])^power
        power *= exponent
    end
    return p
end

function fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f)
    for ef in 1:szE, af in 1:szAge, em in 1:szE, am in 1:szAge
        auxmat[am, em, af, ef] = Π[am, em, af, ef] * sqrt(Sm[am, em]*Sf[af, ef]) /
                insideproduct_den(am, em, af, ef; Sm=Sm, Sf=Sf)
    end
    return auxmat
end

function MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f)
    for em in 1:szE, am in 1:szAge
        s = zero(eltype(μ0m))
        for ef in 1:szE, af in 1:szAge
            s += auxmat[am, em, af, ef] * insideproduct_num(am, em, af, ef; μ0m=μ0m, μ0f=μ0f)
        end
        res[am, em] = Sm[am, em] - μ0m[am, em] - s
    end
    return res
end

function MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m)
    for ef in 1:szE, af in 1:szAge
        s = zero(eltype(μ0f))
        for em in 1:szE, am in 1:szAge
            s += auxmat[am, em, af, ef] * insideproduct_num(am, em, af, ef; μ0m=μ0m, μ0f=μ0f)
        end
        res[af, ef] = Sf[af, ef] - μ0f[af, ef] - s
    end
    return res
end

function fullsys!(resm, resf; auxmat=auxmat, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f)
    fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f)

    MMsysM!(resm, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f)
    MMsysF!(resf, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m)
    return resm, resf
end

Π = 1000 .* rand(szAge, szE, szAge, szE)
auxmat = similar(Π)

Sm = 1000 .* rand(szAge, szE)
Sf = 1000 .* rand(szAge, szE)

μ0m = rand(szAge, szE)
μ0f = rand(szAge, szE)

resm = similar(μ0m)
resf = similar(μ0f)

```

for reference, here’s the benchmark for this version:

```julia
julia> @benchmark fullsys!($resm, $resf; auxmat=$auxmat, Sm=$Sm, Sf=$Sf, μ0m=$μ0m, μ0f=$μ0f)
BenchmarkTools.Trial: 625 samples with 1 evaluation.
 Range (min … max): 6.852 ms … 13.761 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 7.932 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 7.996 ms ± 515.543 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

                  ▄▂▃▄█▁▃▅▅▂▅▃▂ ▁ ▂                           
  ▂▂▂▂▃▃▃▄▄▃▄▇▆▄▆█████████████████▇█▇▆▅▅▅▆▅▃▃▃▃▃▃▂▃▂▄▃▃▁▁▂▁▁▂ ▄
  6.85 ms Histogram: frequency by time 9.47 ms <

 Memory estimate: 96 bytes, allocs estimate: 2.

```

---

<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:** [January 20, 2022, 1:05am UTC](https://discourse.julialang.org/t/whats-a-good-and-easy-way-to-parallelize-this/74893/2 "2022-01-20T01:05:45Z")

</div>

Before trying to speed your code up via paralellization, you should try to optimize it. For example, `insideproduct_num` can be calculated in the log domain, which would remove the floating point powers which would make your code roughly 20x faster.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [January 20, 2022, 1:25am UTC](https://discourse.julialang.org/t/whats-a-good-and-easy-way-to-parallelize-this/74893/3 "2022-01-20T01:25:30Z")

</div>

Good idea, thank you!

```julia
julia> @benchmark insideproduct_den(30, 2, 25, 2; Sm=$Sm, Sf=$Sf)
BenchmarkTools.Trial: 10000 samples with 243 evaluations.
 Range (min … max): 283.642 ns … 939.449 ns ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 322.564 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 328.856 ns ± 44.206 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

      ▄ ▆▇█▄▂ ▁ ▁ ▁
  ▇▂█▄█▆███████▇▇▇█▆▆█▇█▆▆▇▇▇▆▆█▆▅▆▆▆▆▆▅▆▆▆▅▆▆▆▆▅▆▅▅▅▅▆▄▅▃▄▅▄▅▅ █
  284 ns Histogram: log(frequency) by time 550 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark insideproduct_den_log(30, 2, 25, 2; Sm=$Sm, Sf=$Sf)
BenchmarkTools.Trial: 10000 samples with 873 evaluations.
 Range (min … max): 126.515 ns … 1.086 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 143.638 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 147.881 ns ± 24.177 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

   ▁▂ ▄▆▇█▃▂▂▁▁▁ ▁ ▁ ▂
  ▆██▁▇██████████████████▇█▇███▇▇█▇▇▇▇▇▇▇▆▇▇▆▇▆▆▆▄▆▄▅▆▅▄▅▅▆▅▅▅ █
  127 ns Histogram: log(frequency) by time 240 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

so now it’s

```julia
const szAge = 64 - 19 + 1
const szE = 2
const β = 0.92
const δ = 0.10

maxT(am, af; Z=realage(szAge)) = Z - max(am, af) + 1

realage(x) = x + 18
realageinv(x) = x - 18

function insideproduct_num_log(am, em, af, ef; μ0m=μ0m, μ0f=μ0f)
    exponent = 0.5*(β*(1 - δ))
    power = 1.0
    p = zero(promote_type(eltype(μ0m), eltype(μ0f)))
    for k in 0:maxT(realage(am), realage(af))-1
        p += power * log(abs(μ0m[am + k, em] * μ0f[af + k, ef]))
        power *= exponent
    end
    return exp(p)
end

function insideproduct_den_log(am, em, af, ef; Sm=Sm, Sf=Sf)
    exponent = 0.5*(β*(1 - δ))
    power = 1.0
    p = zero(promote_type(eltype(Sm), eltype(Sf)))
    for k in 0:maxT(realage(am), realage(af))-1
        p += power * log(abs(Sm[am + k, em] * Sf[af + k, ef]))
        power *= exponent
    end
    return exp(p)
end

function fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f)
    for ef in 1:szE, af in 1:szAge, em in 1:szE, am in 1:szAge
        auxmat[am, em, af, ef] = Π[am, em, af, ef] * sqrt(Sm[am, em]*Sf[af, ef]) /
                insideproduct_den_log(am, em, af, ef; Sm=Sm, Sf=Sf)
    end
    return auxmat
end

function MMsysM!(res, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f)
    for em in 1:szE, am in 1:szAge
        s = zero(eltype(μ0m))
        for ef in 1:szE, af in 1:szAge
            s += auxmat[am, em, af, ef] * insideproduct_num_log(am, em, af, ef; μ0m=μ0m, μ0f=μ0f)
        end
        res[am, em] = Sm[am, em] - μ0m[am, em] - s
    end
    return res
end

function MMsysF!(res, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m)
    for ef in 1:szE, af in 1:szAge
        s = zero(eltype(μ0f))
        for em in 1:szE, am in 1:szAge
            s += auxmat[am, em, af, ef] * insideproduct_num_log(am, em, af, ef; μ0m=μ0m, μ0f=μ0f)
        end
        res[af, ef] = Sf[af, ef] - μ0f[af, ef] - s
    end
    return res
end

function fullsys!(resm, resf; auxmat=auxmat, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f)
    fillauxmat!(auxmat; Π=Π, Sm=Sm, Sf=Sf, μ0m=μ0m, μ0f=μ0f)

    MMsysM!(resm, μ0m; auxmat=auxmat, Sm=Sm, μ0f=μ0f)
    MMsysF!(resf, μ0f; auxmat=auxmat, Sf=Sf, μ0m=μ0m)
    return resm, resf
end

Π = 1000 .* rand(szAge, szE, szAge, szE)
auxmat = similar(Π)

Sm = 1000 .* rand(szAge, szE)
Sf = 1000 .* rand(szAge, szE)

μ0m = rand(szAge, szE)
μ0f = rand(szAge, szE)

resm = similar(μ0m)
resf = similar(μ0f)

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [January 20, 2022, 5:35am UTC](https://discourse.julialang.org/t/whats-a-good-and-easy-way-to-parallelize-this/74893/4 "2022-01-20T05:35:21Z")

</div>

Why is this a problem?

In other words: I don’t believe the overhead of parallelization is justified by this workload.