Stan Math Library  2.20.0
reverse mode automatic differentiation
cholesky_decompose.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_OPENCL_CHOLESKY_DECOMPOSE_HPP
2 #define STAN_MATH_OPENCL_CHOLESKY_DECOMPOSE_HPP
3 #ifdef STAN_OPENCL
16 #include <CL/cl.hpp>
17 #include <algorithm>
18 #include <cmath>
19 
20 namespace stan {
21 namespace math {
43 inline void cholesky_decompose(matrix_cl& A) {
44  if (A.rows() == 0)
45  return;
46  // Repeats the blocked cholesky decomposition until the size of the remaining
47  // submatrix is smaller or equal to the minimum blocks size
48  // or a heuristic of 100.
49  // The Cholesky (OpenCL) algorithm only uses one local block so we need the
50  // matrix To be less than the max thread block size.
52  try {
54  cl::NDRange(A.rows()), A, A.rows());
55  } catch (const cl::Error& e) {
56  check_opencl_error("cholesky_decompose", e);
57  }
58  return;
59  }
60  // NOTE: The code in this section follows the naming conventions
61  // in the report linked in the docs.
62  const int block
64  // Subset the top left block of the input A into A_11
65  matrix_cl A_11(block, block);
66  A_11.sub_block(A, 0, 0, 0, 0, block, block);
67  // The following function either calls the
68  // blocked cholesky recursively for the submatrix A_11
69  // or calls the kernel directly if the size of the block is small enough
70  cholesky_decompose(A_11);
71  // Copies L_11 back to the input matrix
72  A.sub_block(A_11, 0, 0, 0, 0, block, block);
73 
74  const int block_subset = A.rows() - block;
75  matrix_cl A_21(block_subset, block);
76  A_21.sub_block(A, block, 0, 0, 0, block_subset, block);
77  // computes A_21*((L_11^-1)^T)
78  // and copies the resulting submatrix to the lower left hand corner of A
79  matrix_cl L_21
80  = opencl::multiply<TriangularViewCL::Entire, TriangularViewCL::Upper>(
81  A_21, transpose(tri_inverse<TriangularViewCL::Lower>(A_11)));
82  A.sub_block(L_21, 0, 0, block, 0, block_subset, block);
83  matrix_cl A_22(block_subset, block_subset);
84  A_22.sub_block(A, block, block, 0, 0, block_subset, block_subset);
85  // computes A_22 - L_21*(L_21^T)
86  matrix_cl L_22 = A_22 - multiply_transpose(L_21);
87  // copy L_22 into A's lower left hand corner
88  cholesky_decompose(L_22);
89  A.sub_block(L_22, 0, 0, block, block, block_subset, block_subset);
90 }
91 
92 } // namespace math
93 } // namespace stan
94 
95 #endif
96 #endif
The API to access the methods and values in opencl_context_base.
The matrix_cl class - allocates memory space on the OpenCL device, functions for transfering matrices...
opencl_context_base::tuning_struct & tuning_opts()
Returns the thread block size for the Cholesky Decompositions L_11.
void cholesky_decompose(matrix_cl &A)
Performs an in-place of the the lower-triangular Cholesky factor (i.e., matrix square root) of the sp...
void sub_block(const matrix_cl &A, size_t A_i, size_t A_j, size_t this_i, size_t this_j, size_t nrows, size_t ncols)
Write the contents of A into this starting at the top left of this
Definition: sub_block.hpp:28
Represents a matrix on the OpenCL device.
Definition: matrix_cl.hpp:29
const kernel_cl< in_out_buffer, int > cholesky_decompose("cholesky_decompose", {indexing_helpers, cholesky_decompose_kernel_code})
See the docs for cholesky_decompose() .
checking OpenCL error numbers
Initialization for OpenCL:
matrix_cl transpose(const matrix_cl &src)
Takes the transpose of the matrix on the OpenCL device.
Definition: transpose.hpp:20
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:12
matrix_cl multiply_transpose(const matrix_cl &A)
Computes the product of a square OpenCL matrix with its transpose.
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > block(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, size_t i, size_t j, size_t nrows, size_t ncols)
Return a nrows x ncols submatrix starting at (i-1, j-1).
Definition: block.hpp:22
void check_opencl_error(const char *function, const cl::Error &e)
Throws the domain error with specifying the OpenCL error that occured.

     [ Stan Home Page ] © 2011–2018, Stan Development Team.