Stream: helpdesk (published)

Topic: Parallelisation changes result (unexpectadly)

Timothy (Jul 06 2021 at 10:17):

Hello,

In the following code, I have tried parallelising the for (j, pair) loop, however although the results of one iteration should be independent of other iterations, whenever I add Threads.@threads the result I get changes :frown:

    tightnesses = Vector{Float64}(undef, length(𝐑))
    besttightness = Inf
    @simpleprog length(ixy) * length(𝐑) for (i, x, y) ∈ ixy
        D₁ = copy(D)
        pairs = D[[x,y],:]
        pairs′ = [R * pairs for R ∈ 𝐑]
        for (j, pair) ∈ (enumerate(pairs′) |> collect)
            proginc!()
            D₁[[x,y],:] = pair
            𝚡 = basisparameterise(D₁, 𝐱cent𝛟, 𝐱𝛟)
            tightnesses[j] = basistightness(𝚡, 𝐱; k, krepr)
        end
        if minimum(tightnesses) < besttightness
            besttightness = minimum(tightnesses)
            D[[x,y],:] = pairs′[argmin(tightnesses)]
        end
        𝐃[:,:,i+1] = copy(D)
    end

If anyone has any ideas what's happening and/or how I can get this to work, that would be much appreciated!

Timothy (Jul 06 2021 at 10:18):

For reference length(𝐑) is 360.

Andrey Oskin (Jul 06 2021 at 10:40):

Many questions arise.

What is ixy? What is @simpleprog? What is proginc!? What is D? Is there any difference between D and 𝐃? What is basisparametrise? What is 𝐱cent𝛟? And so on....

Any of this things can potentially make loop iterations-dependent, and it is impossible to say which one just by looking at this code.

One thing, though - you make D1 as a copy of D, and at the same time you change D inside iteration (in if minimum... branch). So, it looks like iterations are not independent.

Timothy (Jul 06 2021 at 11:07):

Unfortunately I the full context is rather ... large. However, I'll attempt to give an idea of what's going on here:

• ixy is a vector of tuples [(1, 1, 2), (2, 1, 3), (3, 1, 4), ..., (26, 2, 13), ..., (240, 15, 16)]
• 𝐑 is a vector of 2D rotation matrices
• @simpleprog n body is a simple macro that generates a progress output for body which is incremented by proginc!()
• D is just some ortonormal matrix
• 𝐃 is a tensor recording the evolution of D
• basisparameterise is a pure function (i.e. no state) which is completely determanistic
• basistightness is also a pure function, but has a very small random element
• 𝐱cent𝛟 and 𝐱𝛟 are another can of worms, but they're constant values

Regarding your comment:

One thing, though - you make D1 as a copy of D, and at the same time you change D inside iteration (in if minimum... branch). So, it looks like iterations are not independent.

D is changed in the outer loop for (i, x, y), I'm talking about the inner loop for (j, pair).

To further clarify regarding that inner loop:

• proginc!() has no effect on the result
• D₁ and tightnesses[j] are the only values affected by it
• all other values are (within that loop) constants

Timothy (Jul 06 2021 at 11:10):

Perhaps it could help if I make a "reduced" version?

besttightness = Inf
for x ∈ xs, y ∈ ys
    D₁ = copy(D)
    pairs = D[[x,y],:]
    pairs′ = makerotations(pairs)
    for (j, pair) ∈ (enumerate(pairs′) |> collect)
        D₁[[x,y],:] = pair
        tightnesses[j] = somepurefunction(D₁)
    end
    if minimum(tightnesses) < besttightness
        besttightness = minimum(tightnesses)
        D[[x,y],:] = pairs′[argmin(tightnesses)]
    end
end

Timothy (Jul 06 2021 at 11:12):

I also tried adding D₁ = copy(D) before D₁[[x,y],:] = pair at some point, but that didn't improve things.

Argel Ramírez Reyes (Jul 08 2021 at 18:28):

a minimal working example would be better but definitely it seems that that line is the one causing you problems.

By the way, it seems like the collect in there is not necessary:

julia> a = [(1,2) (1,3) (1,4)]
1×3 Matrix{Tuple{Int64, Int64}}:
 (1, 2)  (1, 3)  (1, 4)

julia> for (i,j) in enumerate(a)
       println(i," ",j)
       end
1 (1, 2)
2 (1, 3)
3 (1, 4)

julia> for (i,j) in (enumerate(a) |> collect)
       println(i," ",j)
       end
1 (1, 2)
2 (1, 3)
3 (1, 4)

Timothy (Jul 11 2021 at 04:46):

Which line are you referring to? Oh, and collect is needed when using @threads.

Gunnar Farnebäck (Jul 11 2021 at 09:22):

    for (j, pair) ∈ (enumerate(pairs′) |> collect)
        D₁[[x,y],:] = pair
        tightnesses[j] = somepurefunction(D₁)
    end

I don't see how this possibly could be parallelized. All threads will try to update the contents of D₁ and will overwrite each other in a data race. You will need to set things up so each thread has its own copy of D₁ to play with.

Last updated: Oct 02 2023 at 04:34 UTC