Description
The cholesky_decompose kernel (stan/math/opencl/kernels/cholesky_decompose.hpp) runs as one work group with one work item per row. The work items pass intermediate results to each other through the __global buffer A: work item 0 writes the diagonal entry of column j, and after a barrier every work item i > j reads that entry and row j, which work item j wrote in earlier steps. The barriers are barrier(CLK_LOCAL_MEM_FENCE). That flag orders operations on local memory only; reading buffer contents written by other work items before the barrier needs CLK_GLOBAL_MEM_FENCE.
On an AMD Radeon integrated GPU (OpenCL device gfx1035, AMD-APP 3652.0, Windows) the difference is visible. As soon as the work group has more than 64 work items, some of them read entries of A as they were before other work items overwrote them. The factor is wrong, by more than 100 % in some entries, and no error is raised. The blocked algorithm in stan/math/opencl/cholesky_decompose.hpp passes blocks of up to 256 rows to this kernel, so any matrix_cl Cholesky decomposition with more than 64 rows is affected. With a standalone copy of the kernel, 60 of 60 runs with N between 65 and 256 were wrong and 40 of 40 runs with N up to 64 were correct; with barrier(CLK_GLOBAL_MEM_FENCE) all 100 were correct.
With PoCL and with the Intel OpenCL CPU runtime the kernel gives correct results as it is.
The tests cholesky_decompose_small, _big and _big_tuning_opts in test/unit/math/opencl/cholesky_decompose_test.cpp do not run the kernel on any device: they call the Eigen overload of cholesky_decompose, which no longer dispatches to OpenCL, so they compare Eigen with Eigen.
Example
int N = 200;
Eigen::MatrixXd m = Eigen::MatrixXd::Random(N, N);
Eigen::MatrixXd a = m * m.transpose() + N * Eigen::MatrixXd::Identity(N, N);
stan::math::matrix_cl<double> a_cl(a);
stan::math::opencl::cholesky_decompose(a_cl);
Eigen::MatrixXd diff
= stan::math::from_matrix_cl(a_cl) - stan::math::cholesky_decompose(a);
std::cout << diff.cwiseAbs().maxCoeff() << std::endl;Expected Output
Expected: a difference at rounding level.
Current, on the GPU above, in five runs of this computation as a unit test (MinGW g++ 15.2): 13,852 of the 40,000 entries differ by more than 1e-8. The largest difference is 2.9 to 4.3, on a diagonal entry whose correct value is about 16.
Description
The
cholesky_decomposekernel (stan/math/opencl/kernels/cholesky_decompose.hpp) runs as one work group with one work item per row. The work items pass intermediate results to each other through the__globalbufferA: work item 0 writes the diagonal entry of columnj, and after a barrier every work itemi > jreads that entry and rowj, which work itemjwrote in earlier steps. The barriers arebarrier(CLK_LOCAL_MEM_FENCE). That flag orders operations on local memory only; reading buffer contents written by other work items before the barrier needsCLK_GLOBAL_MEM_FENCE.On an AMD Radeon integrated GPU (OpenCL device
gfx1035, AMD-APP 3652.0, Windows) the difference is visible. As soon as the work group has more than 64 work items, some of them read entries ofAas they were before other work items overwrote them. The factor is wrong, by more than 100 % in some entries, and no error is raised. The blocked algorithm instan/math/opencl/cholesky_decompose.hpppasses blocks of up to 256 rows to this kernel, so anymatrix_clCholesky decomposition with more than 64 rows is affected. With a standalone copy of the kernel, 60 of 60 runs with N between 65 and 256 were wrong and 40 of 40 runs with N up to 64 were correct; withbarrier(CLK_GLOBAL_MEM_FENCE)all 100 were correct.With PoCL and with the Intel OpenCL CPU runtime the kernel gives correct results as it is.
The tests
cholesky_decompose_small,_bigand_big_tuning_optsintest/unit/math/opencl/cholesky_decompose_test.cppdo not run the kernel on any device: they call the Eigen overload ofcholesky_decompose, which no longer dispatches to OpenCL, so they compare Eigen with Eigen.Example
Expected Output
Expected: a difference at rounding level.
Current, on the GPU above, in five runs of this computation as a unit test (MinGW g++ 15.2): 13,852 of the 40,000 entries differ by more than 1e-8. The largest difference is 2.9 to 4.3, on a diagonal entry whose correct value is about 16.