Add missing synchronization calls in CUDA version of LINCS

In CUDA, the same shared memory buffer is used to store different
intermediate values. To avoid overwriting these values prematurely,
extra blocking synchronizations are needed in the GPU kernel.

This bug was not exposed neither by regression tests nor by unit tests
because it was affecting values rarely and the deviation were within
tolerances.

diff --git a/src/gromacs/mdlib/lincs_gpu_internal.cu b/src/gromacs/mdlib/lincs_gpu_internal.cu
index ccd9e0a6a7c20041e860be1757cfb74b6c9003c7..487aa8c01ee490bc6ee9c52e83c5c2db71412b18 100644 (file)
--- a/src/gromacs/mdlib/lincs_gpu_internal.cu
+++ b/src/gromacs/mdlib/lincs_gpu_internal.cu
@@ -200,6 +200,9 @@ __launch_bounds__(c_maxThreadsPerBlock) __global__
      *  Inverse matrix using a set of expansionOrder matrix multiplications
      */
 
+    // Make sure that we don't overwrite the sm_r[..] array.
+    __syncthreads();
+
     // This will use the same memory space as sm_r, which is no longer needed.
     extern __shared__ float sm_rhs[];
     // Save current right-hand-side vector in the shared memory
@@ -333,7 +336,8 @@ __launch_bounds__(c_maxThreadsPerBlock) __global__
         // 6 values are saved. Dummy threads will have zeroes in their virial: targetLength,
         // lagrangeScaled and rc are all set to zero for them in the beginning of the kernel.
         // The sm_threadVirial[..] will overlap with the sm_r[..] and sm_rhs[..], but the latter
-        // two are no longer in use.
+        // two are no longer in use, which we make sure by waiting for all threads in block.
+        __syncthreads();
         extern __shared__ float sm_threadVirial[];
         float                   mult                  = targetLength * lagrangeScaled;
         sm_threadVirial[0 * blockDim.x + threadIdx.x] = mult * rc.x * rc.x;