From d075c864f24602a99e2d0c66c0b97f3742f99830 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sun, 4 Oct 2026 12:56:03 +0200 Subject: [PATCH 1/7] 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/7] 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/7] 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/7] 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/7] 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 From 65268281454db07e9967ecb389ea245aa2117668 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Wed, 21 Oct 2026 19:32:14 +0200 Subject: [PATCH 6/7] Share row j through local memory in the OpenCL cholesky_decompose kernel Work item j now computes the diagonal entry of row j and copies the row into a local buffer, from which the other work items read it. Each work item then accesses only its own row of A in global memory, so local memory fences suffice and the global fences, which slowed the kernel down, are no longer needed. Reading row j from local memory also makes the kernel faster than develop's on an AMD GPU. Suggested by Andrew Johnson in review. --- stan/math/opencl/cholesky_decompose.hpp | 3 ++- .../opencl/kernels/cholesky_decompose.hpp | 25 +++++++++++-------- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/stan/math/opencl/cholesky_decompose.hpp b/stan/math/opencl/cholesky_decompose.hpp index 41a5ab77b65..b03024dbb58 100644 --- a/stan/math/opencl/cholesky_decompose.hpp +++ b/stan/math/opencl/cholesky_decompose.hpp @@ -52,7 +52,8 @@ inline void cholesky_decompose(matrix_cl& A) { if (A.rows() <= opencl_context.tuning_opts().cholesky_min_L11_size) { try { opencl_kernels::cholesky_decompose(cl::NDRange(A.rows()), - cl::NDRange(A.rows()), A, A.rows()); + cl::NDRange(A.rows()), A, A.rows(), + cl::Local(sizeof(double) * A.rows())); } catch (const cl::Error& e) { check_opencl_error("cholesky_decompose", e); } diff --git a/stan/math/opencl/kernels/cholesky_decompose.hpp b/stan/math/opencl/kernels/cholesky_decompose.hpp index 6e5736e6472..5cf95e71591 100644 --- a/stan/math/opencl/kernels/cholesky_decompose.hpp +++ b/stan/math/opencl/kernels/cholesky_decompose.hpp @@ -24,37 +24,42 @@ static constexpr const char* cholesky_decompose_kernel_code = STRINGIFY( * @param[in, out] A The input matrix and the result of the cholesky * decomposition * @param rows The number of rows for A and B. + * @param row_j Local memory buffer of rows doubles that holds row j. * @note Code is a const char* held in * cholesky_decompose_kernel_code. * Used in math/opencl/cholesky_decompose.hpp. * This kernel uses the helper macros available in helpers.cl. * */ - __kernel void cholesky_decompose(__global double* A, int rows) { + __kernel void cholesky_decompose(__global double* A, int rows, + __local double* row_j) { const int local_index = get_local_id(0); // The following code is the sequential version of the inplace // cholesky decomposition. Only the innermost loops are parallelized. The // rows are processed sequentially. This loop process all the rows: for (int j = 0; j < rows; j++) { - if (local_index == 0) { + if (local_index == j) { double sum = 0.0; for (int k = 0; k < j; k++) { - sum = sum + A(j, k) * A(j, k); + const double a = A(j, k); + row_j[k] = a; + sum = sum + a * a; } A(j, j) = sqrt(A(j, j) - sum); + row_j[j] = A(j, j); } - // Work items read entries of A written by other work items, so the - // barriers need a global memory fence. - barrier(CLK_GLOBAL_MEM_FENCE); + // Each work item accesses only its own row of A and reads row j from + // local memory, so the barriers need only a local memory fence. + barrier(CLK_LOCAL_MEM_FENCE); if (local_index < j) { A(local_index, j) = 0.0; } else if (local_index > j) { double sum = 0.0; for (int k = 0; k < j; k++) - sum = sum + A(local_index, k) * A(j, k); - A(local_index, j) = (A(local_index, j) - sum) / A(j, j); + sum = sum + A(local_index, k) * row_j[k]; + A(local_index, j) = (A(local_index, j) - sum) / row_j[j]; } - barrier(CLK_GLOBAL_MEM_FENCE); + barrier(CLK_LOCAL_MEM_FENCE); } } // \cond @@ -65,7 +70,7 @@ static constexpr const char* cholesky_decompose_kernel_code = STRINGIFY( * See the docs for \link kernels/cholesky_decompose.hpp cholesky_decompose() * \endlink */ -const kernel_cl cholesky_decompose( +const kernel_cl cholesky_decompose( "cholesky_decompose", {indexing_helpers, cholesky_decompose_kernel_code}); } // namespace opencl_kernels From a12718b149103cdfbe3a4e8074a35a242706a01d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Wed, 21 Oct 2026 22:26:29 +0200 Subject: [PATCH 7/7] Use += for the sums in the OpenCL cholesky_decompose kernel --- stan/math/opencl/kernels/cholesky_decompose.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/opencl/kernels/cholesky_decompose.hpp b/stan/math/opencl/kernels/cholesky_decompose.hpp index 5cf95e71591..ff8c3d7be70 100644 --- a/stan/math/opencl/kernels/cholesky_decompose.hpp +++ b/stan/math/opencl/kernels/cholesky_decompose.hpp @@ -43,7 +43,7 @@ static constexpr const char* cholesky_decompose_kernel_code = STRINGIFY( for (int k = 0; k < j; k++) { const double a = A(j, k); row_j[k] = a; - sum = sum + a * a; + sum += a * a; } A(j, j) = sqrt(A(j, j) - sum); row_j[j] = A(j, j); @@ -56,7 +56,7 @@ static constexpr const char* cholesky_decompose_kernel_code = STRINGIFY( } else if (local_index > j) { double sum = 0.0; for (int k = 0; k < j; k++) - sum = sum + A(local_index, k) * row_j[k]; + sum += A(local_index, k) * row_j[k]; A(local_index, j) = (A(local_index, j) - sum) / row_j[j]; } barrier(CLK_LOCAL_MEM_FENCE);