Fix race conditions caused by missing fences in OpenCL kernels for cholesky_decompose, tridiagonalization and diag_inv

stan/math/opencl/kernels/cholesky_decompose.hpp

@@ -43,7 +43,9 @@ static constexpr const char* cholesky_decompose_kernel_code = STRINGIFY(
           }
           A(j, j) = sqrt(A(j, j) - sum);
         }
-        barrier(CLK_LOCAL_MEM_FENCE);
+        // Work items read entries of A written by other work items, so the
+        // barriers need a global memory fence.
+        barrier(CLK_GLOBAL_MEM_FENCE);
         if (local_index < j) {
           A(local_index, j) = 0.0;
         } else if (local_index > j) {
@@ -52,7 +54,7 @@ static constexpr const char* cholesky_decompose_kernel_code = STRINGIFY(
             sum = sum + A(local_index, k) * A(j, k);
           A(local_index, j) = (A(local_index, j) - sum) / A(j, j);
         }
-        barrier(CLK_LOCAL_MEM_FENCE);
+        barrier(CLK_GLOBAL_MEM_FENCE);
       }
     }
     // \cond

stan/math/opencl/kernels/diag_inv.hpp

@@ -72,6 +72,9 @@ static constexpr const char* diag_inv_kernel_code = STRINGIFY(
         }
         barrier(CLK_LOCAL_MEM_FENCE);
       }
+      // The copy below overwrites entries of A that other work items read
+      // above, so a global memory fence separates the two.
+      barrier(CLK_GLOBAL_MEM_FENCE);
       for (int j = 0; j < block_size; j++) {
         // Each thread copies one column.
         A(A_offset + j, A_offset + index) = tmp_inv[tmp_offset + j];

stan/math/opencl/kernels/tridiagonalization.hpp

@@ -55,7 +55,9 @@ static constexpr const char* tridiagonalization_householder_kernel_code
           // calculate column norm between threads
           __local double q_local[LOCAL_SIZE_];
           q_local[lid] = q;
-          barrier(CLK_LOCAL_MEM_FENCE);
+          // Work items below read and write entries of P that other work
+          // items wrote above, so this barrier also fences global memory.
+          barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
           for (int step = lsize / REDUCTION_STEP_SIZE; step > 0;
                step /= REDUCTION_STEP_SIZE) {
             if (lid < step) {

test/unit/math/opencl/cholesky_decompose_test.cpp

@@ -29,6 +29,20 @@ TEST(MathMatrixOpenCL, cholesky_decompose_cpu_vs_cl_small) {
   EXPECT_MATRIX_NEAR(m1, m1_res, 1e-8);
 }
 
+TEST(MathMatrixOpenCL, cholesky_decompose_cpu_vs_cl_one_work_group) {
+  int size = 200;
+  stan::math::matrix_d m = stan::math::matrix_d::Random(size, size);
+  stan::math::matrix_d m_pos_def
+      = m * m.transpose() + size * Eigen::MatrixXd::Identity(size, size);
+
+  stan::math::matrix_cl<double> m_cl(m_pos_def);
+  stan::math::matrix_d m_res = stan::math::cholesky_decompose(m_pos_def);
+
+  stan::math::opencl::cholesky_decompose(m_cl);
+
+  EXPECT_MATRIX_NEAR(m_res, stan::math::from_matrix_cl(m_cl), 1e-8);
+}
+
 namespace {
 inline void cholesky_decompose_test(int size) {
   stan::math::matrix_d m1 = stan::math::matrix_d::Random(size, size);