Reductions with CUDA are sometimes wrong

[Sbozzolo]

Starting mid December, sometimes, some of my MPI reductions produce wrong results when I am working with CUDA (and only when working with CUDA).

I came to this conclusion via printf-debugging:

println("PRE", mypid, " ", Array(remapper._interpolated_values)[end, end, end, 1]); flush(stdout)
MPI.Barrier(remapper.comms_ctx.mpicomm)
MPI.Reduce!(remapper._interpolated_values[:, :, :, 1], dest, +, mpicomm)
MPI.Barrier(remapper.comms_ctx.mpicomm)
println("POST", mypid, " ", Array(dest)[end, end, end]); flush(stdout)

that produces:

PRE1 0.0
PRE2 500.0
POST1 249.99999999999997 
POST2 0.0 

The value 249.99999999999997 is wrong and it should be 500 (0.0 + 500.0 = 500.0). I tried forcing CUDA.synchronize and/or CUDA.sync, and sometimes results change (for example, I can get 450.0 instead of 250). Note also that 250 could be the correct value for another point.

The code above is embedded in a more complex test case, and I can reproduce the problem reliably when running in a relatively isolated environment. When I work interactively or when I change the test, the test becomes very flaky and extremely hard to reproduce. To ensure consistent results, I have to clean the depot.

I initially thought it was a buggy MPI implementation, but it seems that this problem also appears in production on a different machine (with simulations that run for multiple days). On that machine, the same test never triggers the bug (as checked by running the test 100 times)

Considering the bug is hard to reproduce, this is what I think I concluded:

• With CUDA 5.6, MPI 0.20.21 is fine MPI 0.20.22 is not (build)
• With CUDA 5.6, commit aac9688 is not fine (build) but commit 9584ac8 is fine (build)
• Adding local in front of datatype as suggested on Slack does not fix the issue (build)
• With CUDA 5.5, MPI 0.20.22 is fine (build)
• Both Julia 1.10 and Julia 1.11 have the problem

[giordano]

Would you be able to craft a reproducer? That would help a lot with tracking down the issue.

[Sbozzolo]

Would you be able to craft a reproducer? That would help a lot with tracking down the issue.

The only way I found to reliably reproduce the bug is through our buildkite on our cluster with a clean depot. If I run the same test on another machine, the test always passes. If I run from an interactive shell, any change seem to affect the outcome. When running interactively, the test could go from always failing to always passing seemingly only adding a print statament.

I know it's not the best, but I can probably give you a way to run something in the same way.

I also just ran this test: reverting 877, the build passes.

It's hard to see how #877 could be the problem, so I don't know if I am reading the tea leaves here. I can also bisect #877 further.

[giordano]

I can also bisect #877 further.

Yeah, that could be useful. I'm also struggling to see how that change would have any effect on platforms like x86-64 where the MPI extent and Julia's aligned size agreed already before that PR. Also, the bug observed in #853 was that we were reading garbage memory (because of misaligned sizes of array elements), a factor of 2 of difference is a bit too specific.

[eschnett]

In your example, remapper._interpolated_values[:, :, :, 1] makes a copy of part of the array while slicing it. Can you make the bug appear without slicing? It might be that the bug is related to this copying – either the copy is made too late (missing sync), or the copy is destroyed too early (not rooted for GC), etc. What happens if you assign this copy to a variable? What happens if you use an array view instead?

[vchuravy]

I agree with Erik. That code piece is the most suspect.

One must synchronize the CUDA stream before using CUDA memory in an MPI call. The allocation of the copy and the copy itself are all stream ordered.

You could try adding a synchronize before the barrier and use a @view to avoid the extra CUDA operations.

[giordano]

@Sbozzolo I see you're working around this issue uptream in CliMA/ClimaCore.jl#2169, but did you get an intuition of what could the problem? Is there anything we can do here to prevent similar situations elsewhere?

[Sbozzolo]

(Sorry for the lack of communication over the past days. I was travelling and wanted to write a longer summary for this issue).

I found a few more things:

• My previous bisection was incorrect because I found out that the model of the GPU matters. For my test to be reliable, only P100s can to be used. When anything but a P100 is used, (even if just one of the two GPUs is not a P100), my test does not trigger the problem.
• Armed with this knowledge, I verified that MPI 0.20.22 is fine with CUDA 5.5, but no longer fine with CUDA 5.6
• I changed the line that @eschnett and @vchuravy pointed out to using view and indeed the problem disappeared.

Let me also point out that I did try to synchronize the CUDA streams before the MPI call. I tested something like:

    MPI.barrier(remapper.comms_ctx.mpicomm)
    CUDA.synchronize()
    MPI.Reduce!(
            remapper._interpolated_values[:, :, :, begin],
            dest,
            +,
            remapper.comms_ctx.mpicomm;
            root = 0
        )

This did not solve the issue. (In fact, most of my tests were running with this)

It is still not clear to me whether this pattern is supposed to work or not. If not, I can add a note about it to the docs of MPI.jl

[simonbyrne]

I think what @eschnett and @vchuravy were suggesting was that remapper._interpolated_values[:, :, :, begin] creates a new array, and then copies the contents asynchronously: it is this copy operation that isn't being synchronized.

You could fix it by:

    src = remapper._interpolated_values[:, :, :, begin]
    CUDA.synchronize()
    MPI.Reduce!(
            src,
            dest,
            +,
            remapper.comms_ctx.mpicomm;
            root = 0
        )

or by using a view:

    MPI.Reduce!(
            @view(remapper._interpolated_values[:, :, :, begin]),
            dest,
            +,
            remapper.comms_ctx.mpicomm;
            root = 0
        )

[Sbozzolo]

I think what @eschnett and @vchuravy were suggesting was that remapper._interpolated_values[:, :, :, begin] creates a new array, and then copies the contents asynchronously: it is this copy operation that isn't being synchronized.

You could fix it by:

src = remapper._interpolated_values[:, :, :, begin]
CUDA.synchronize()
MPI.Reduce!(
        src,
        dest,
        +,
        remapper.comms_ctx.mpicomm;
        root = 0
    )

or by using a view:

MPI.Reduce!(
        @view(remapper._interpolated_values[:, :, :, begin]),
        dest,
        +,
        remapper.comms_ctx.mpicomm;
        root = 0
    )

Yes, sorry if I didn't make it clear in my previous message. I changed the code to use view, which solves the problem.

Do I understand correctly that MPI.Reduce!(remapper._interpolated_values[:, :, :, begin], ...) has no guarantee of working correctly (and in fact, I should expect it not to work correctly)?

[giordano]

Do I understand correctly that MPI.Reduce!(remapper._interpolated_values[:, :, :, begin], ...) has no guarantee of working correctly (and in fact, I should expect it not to work correctly)?

I presume the problem is the fact that this is data which lives on the GPU rather than the CPU? In julia, slicing a regular Array would happen before calling the MPI.Reduce! function and you shouldn't have any problems, but for GPU arrays syncing may be different, that's why CUDA.synchronize() would help.

But also note that slicing always copies, so depending on what you want to do you may be doing something wrong.