From d075c864f24602a99e2d0c66c0b97f3742f99830 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sun, 4 Oct 2026 12:56:03 +0200 Subject: [PATCH 1/5] Test OpenCL cholesky_decompose with one work group of 200 work items --- .../math/opencl/cholesky_decompose_test.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/test/unit/math/opencl/cholesky_decompose_test.cpp b/test/unit/math/opencl/cholesky_decompose_test.cpp index c056bc91b61..2a1af66e406 100644 --- a/test/unit/math/opencl/cholesky_decompose_test.cpp +++ b/test/unit/math/opencl/cholesky_decompose_test.cpp @@ -29,6 +29,23 @@ 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) { + // One call of the cholesky_decompose kernel with 200 work items, unless + // cholesky_min_L11_size is smaller. An AMD GPU returned a wrong factor from + // the kernel with more than 64 work items. + 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 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); From d7d1d2ea51570eb0e3d2ebd294cdeb2be83ff55d Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sun, 4 Oct 2026 12:56:03 +0200 Subject: [PATCH 2/5] Fix OpenCL cholesky_decompose kernel: barriers need a global memory fence --- stan/math/opencl/kernels/cholesky_decompose.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/stan/math/opencl/kernels/cholesky_decompose.hpp b/stan/math/opencl/kernels/cholesky_decompose.hpp index cac42c651bd..6e5736e6472 100644 --- a/stan/math/opencl/kernels/cholesky_decompose.hpp +++ b/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 From 0d308b63c58f4bcc48a02c27c5d499f88e40bb83 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sun, 4 Oct 2026 18:12:14 +0200 Subject: [PATCH 3/5] Fix OpenCL tridiagonalization_householder kernel: barrier needs a global memory fence --- stan/math/opencl/kernels/tridiagonalization.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/stan/math/opencl/kernels/tridiagonalization.hpp b/stan/math/opencl/kernels/tridiagonalization.hpp index c8369dbd2fc..265801366bc 100644 --- a/stan/math/opencl/kernels/tridiagonalization.hpp +++ b/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) { From a8fb654c86138d5edf64c0ce507cddb19b473ade Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sun, 4 Oct 2026 18:12:14 +0200 Subject: [PATCH 4/5] Fix OpenCL diag_inv kernel: global memory fence before the final copy --- stan/math/opencl/kernels/diag_inv.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/stan/math/opencl/kernels/diag_inv.hpp b/stan/math/opencl/kernels/diag_inv.hpp index b8af6a9aa5e..7f29ee4ad7f 100644 --- a/stan/math/opencl/kernels/diag_inv.hpp +++ b/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]; From 1af67a2f8a383684830ad611729a37cd692ecdc1 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sun, 4 Oct 2026 21:52:36 +0200 Subject: [PATCH 5/5] Remove the comment in the OpenCL cholesky_decompose one work group test --- test/unit/math/opencl/cholesky_decompose_test.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/unit/math/opencl/cholesky_decompose_test.cpp b/test/unit/math/opencl/cholesky_decompose_test.cpp index 2a1af66e406..ee9a868f3c2 100644 --- a/test/unit/math/opencl/cholesky_decompose_test.cpp +++ b/test/unit/math/opencl/cholesky_decompose_test.cpp @@ -30,9 +30,6 @@ TEST(MathMatrixOpenCL, cholesky_decompose_cpu_vs_cl_small) { } TEST(MathMatrixOpenCL, cholesky_decompose_cpu_vs_cl_one_work_group) { - // One call of the cholesky_decompose kernel with 200 work items, unless - // cholesky_min_L11_size is smaller. An AMD GPU returned a wrong factor from - // the kernel with more than 64 work items. int size = 200; stan::math::matrix_d m = stan::math::matrix_d::Random(size, size); stan::math::matrix_d m_pos_def