1 #ifndef STAN_MATH_REV_MAT_FUN_CHOLESKY_DECOMPOSE_HPP 2 #define STAN_MATH_REV_MAT_FUN_CHOLESKY_DECOMPOSE_HPP 39 for (
size_type j = 0; j < L.cols(); ++j) {
40 for (
size_type i = j; i < L.cols(); ++i) {
41 L.coeffRef(i, j).vi_ = vari_ref[pos++];
44 L.coeffRef(k, j).vi_ = dummy;
53 typedef Eigen::Block<Eigen::MatrixXd>
Block_;
74 const Eigen::Matrix<double, -1, -1>& L_A)
83 block_size_ =
std::min(block_size_, 128);
86 vari_ref_A_[pos] = A.coeffRef(i, j).vi_;
87 vari_ref_L_[pos] =
new vari(L_A.coeffRef(i, j),
false);
101 using Eigen::StrictlyUpper;
103 L.transposeInPlace();
104 L_adj = (L * L_adj.triangularView<
Lower>()).eval();
105 L_adj.triangularView<StrictlyUpper>()
106 = L_adj.adjoint().triangularView<StrictlyUpper>();
107 L.triangularView<
Upper>().solveInPlace(L_adj);
108 L.triangularView<
Upper>().solveInPlace(L_adj.transpose());
120 using Eigen::MatrixXd;
121 using Eigen::StrictlyUpper;
123 auto L_adj = Eigen::MatrixXd::Zero(M_, M_).eval();
124 auto L = Eigen::MatrixXd::Zero(M_, M_).eval();
128 L_adj.coeffRef(i, j) = vari_ref_L_[pos]->
adj_;
129 L.coeffRef(i, j) = vari_ref_L_[pos]->
val_;
134 for (
int k = M_; k > 0; k -= block_size_) {
135 int j =
std::max(0, k - block_size_);
136 Block_ R = L.block(j, 0, k - j, j);
137 Block_ D = L.block(j, j, k - j, k - j);
138 Block_ B = L.block(k, 0, M_ - k, j);
139 Block_ C = L.block(k, j, M_ - k, k - j);
140 Block_ R_adj = L_adj.block(j, 0, k - j, j);
141 Block_ D_adj = L_adj.block(j, j, k - j, k - j);
142 Block_ B_adj = L_adj.block(k, 0, M_ - k, j);
143 Block_ C_adj = L_adj.block(k, j, M_ - k, k - j);
144 if (C_adj.size() > 0) {
145 C_adj = D.transpose()
146 .triangularView<
Upper>()
147 .solve(C_adj.transpose())
149 B_adj.noalias() -= C_adj * R;
150 D_adj.noalias() -= C_adj.transpose() * C;
152 symbolic_rev(D, D_adj);
153 R_adj.noalias() -= C_adj.transpose() * B;
154 R_adj.noalias() -= D_adj.selfadjointView<
Lower>() * R;
155 D_adj.diagonal() *= 0.5;
156 D_adj.triangularView<StrictlyUpper>().setZero();
161 vari_ref_A_[pos++]->adj_ += L_adj.coeffRef(i, j);
186 const Eigen::Matrix<double, -1, -1>& L_A)
194 size_t accum_i = accum;
198 size_t pos = j + accum_i;
199 vari_ref_A_[pos] = A.coeffRef(i, j).vi_;
200 vari_ref_L_[pos] =
new vari(L_A.coeffRef(i, j),
false);
221 using Eigen::RowMajor;
222 Matrix<double, -1, -1, RowMajor> adjL(M_, M_);
223 Matrix<double, -1, -1, RowMajor> LA(M_, M_);
224 Matrix<double, -1, -1, RowMajor> adjA(M_, M_);
228 adjL.coeffRef(i, j) = vari_ref_L_[pos]->
adj_;
229 LA.coeffRef(i, j) = vari_ref_L_[pos]->
val_;
235 for (
int i = M_ - 1; i >= 0; --i) {
236 for (
int j = i; j >= 0; --j) {
238 adjA.coeffRef(i, j) = 0.5 * adjL.coeff(i, j) / LA.coeff(i, j);
240 adjA.coeffRef(i, j) = adjL.coeff(i, j) / LA.coeff(j, j);
242 -= adjL.coeff(i, j) * LA.coeff(i, j) / LA.coeff(j, j);
244 for (
int k = j - 1; k >= 0; --k) {
245 adjL.coeffRef(i, k) -= adjA.coeff(i, j) * LA.coeff(j, k);
246 adjL.coeffRef(j, k) -= adjA.coeff(i, j) * LA.coeff(i, k);
248 vari_ref_A_[pos--]->
adj_ += adjA.coeffRef(i, j);
276 const Eigen::Matrix<double, -1, -1>& L_A)
286 vari_ref_A_[pos] = A.coeffRef(i, j).vi_;
287 vari_ref_L_[pos] =
new vari(L_A.coeffRef(i, j),
false);
300 L_adj = opencl::multiply<TriangularViewCL::Upper, TriangularViewCL::Entire>(
303 L =
transpose(tri_inverse<TriangularViewCL::Lower>(L));
319 const int packed_size = M_ * (M_ + 1) / 2;
320 std::vector<double> L_adj_cpu(packed_size);
321 std::vector<double> L_val_cpu(packed_size);
323 for (
size_type j = 0; j < packed_size; ++j) {
324 L_adj_cpu[j] = vari_ref_L_[j]->
adj_;
325 L_val_cpu[j] = vari_ref_L_[j]->
val_;
327 matrix_cl L = packed_copy<TriangularViewCL::Lower>(L_val_cpu, M_);
328 matrix_cl L_adj = packed_copy<TriangularViewCL::Lower>(L_adj_cpu, M_);
331 block_size =
std::max(block_size, 8);
337 for (
int k = M_; k > 0; k -= block_size) {
338 const int j =
std::max(0, k - block_size);
339 const int k_j_ind = k - j;
340 const int m_k_ind = M_ - k;
353 D.
sub_block(L, j, j, 0, 0, k_j_ind, k_j_ind);
355 C.
sub_block(L, k, j, 0, 0, m_k_ind, k_j_ind);
357 R_adj.
sub_block(L_adj, j, 0, 0, 0, k_j_ind, j);
358 D_adj.
sub_block(L_adj, j, j, 0, 0, k_j_ind, k_j_ind);
359 B_adj.
sub_block(L_adj, k, 0, 0, 0, m_k_ind, j);
360 C_adj.
sub_block(L_adj, k, j, 0, 0, m_k_ind, k_j_ind);
363 = opencl::multiply<TriangularViewCL::Entire, TriangularViewCL::Lower>(
364 C_adj, tri_inverse<TriangularViewCL::Lower>(D));
365 B_adj = B_adj - C_adj * R;
368 symbolic_rev(D, D_adj);
370 R_adj = R_adj -
transpose(C_adj) * B - D_adj * R;
374 L_adj.sub_block(R_adj, 0, 0, j, 0, k_j_ind, j);
375 L_adj.sub_block(D_adj, 0, 0, j, j, k_j_ind, k_j_ind);
376 L_adj.sub_block(B_adj, 0, 0, k, 0, m_k_ind, j);
377 L_adj.sub_block(C_adj, 0, 0, k, j, m_k_ind, k_j_ind);
379 L_adj_cpu = packed_copy<TriangularViewCL::Lower>(L_adj);
380 for (
size_type j = 0; j < packed_size; ++j) {
381 vari_ref_A_[j]->
adj_ += L_adj_cpu[j];
399 const Eigen::Matrix<var, -1, -1>& A) {
406 Eigen::LLT<Eigen::Ref<Eigen::MatrixXd>,
Eigen::Lower> L_factor(L_A);
413 Eigen::Matrix<var, -1, -1> L(A.rows(), A.cols());
414 if (L_A.rows() <= 35) {
417 size_t accum_i = accum;
418 for (
size_type j = 0; j < L.cols(); ++j) {
419 for (
size_type i = j; i < L.cols(); ++i) {
421 size_t pos = j + accum_i;
425 L.coeffRef(k, j).vi_ = dummy;
virtual void chain()
Reverse mode differentiation algorithm refernce:
int rows(const Eigen::Matrix< T, R, C > &m)
Return the number of rows in the specified matrix, vector, or row vector.
cholesky_block(const Eigen::Matrix< var, -1, -1 > &A, const Eigen::Matrix< double, -1, -1 > &L_A)
Constructor for cholesky function.
int min(const std::vector< int > &x)
Returns the minimum coefficient in the specified column vector.
double value_of_rec(const fvar< T > &v)
Return the value of the specified variable.
void symbolic_rev(Block_ &L, Block_ &L_adj)
Symbolic adjoint calculation for cholesky factor A.
cholesky_scalar(const Eigen::Matrix< var, -1, -1 > &A, const Eigen::Matrix< double, -1, -1 > &L_A)
Constructor for cholesky function.
void check_square(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is square.
The variable implementation base class.
The API to access the methods and values in opencl_context_base.
opencl_context_base::tuning_struct & tuning_opts()
Returns the thread block size for the Cholesky Decompositions L_11.
Independent (input) and dependent (output) variables for gradients.
void triangular_transpose()
Copies a lower/upper triangular of a matrix to it's upper/lower.
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...
const double val_
The value of this variable.
auto multiply(const matrix_cl &A, const matrix_cl &B)
Computes the product of the specified matrices with the option of specifying the triangularity of eit...
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double e.
void zeros()
Stores zeros in the matrix on the OpenCL device.
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
virtual void chain()
Reverse mode differentiation algorithm using OpenCL.
Represents a matrix on the OpenCL device.
int cholesky_size_worth_transfer
cholesky_opencl(const Eigen::Matrix< var, -1, -1 > &A, const Eigen::Matrix< double, -1, -1 > &L_A)
Constructor for OpenCL cholesky function.
void set_lower_tri_coeff_ref(Eigen::Matrix< var, -1, -1 > &L, vari **vari_ref)
Set the lower right triangular of a var matrix given a set of vari**.
int max(const std::vector< int > &x)
Returns the maximum coefficient in the specified column vector.
matrix_cl transpose(const matrix_cl &src)
Takes the transpose of the matrix on the OpenCL device.
int cholesky_rev_min_block_size
int cholesky_rev_block_partition
void symbolic_rev(matrix_cl &L, matrix_cl &L_adj)
Symbolic adjoint calculation for cholesky factor A.
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
virtual void chain()
Reverse mode differentiation algorithm refernce:
This struct always provides access to the autodiff stack using the singleton pattern.
void check_pos_definite(const char *function, const char *name, const Eigen::Matrix< T_y, -1, -1 > &y)
Check if the specified square, symmetric matrix is positive definite.
void check_symmetric(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is symmetric.
matrix_cl diagonal_multiply(const matrix_cl &A, const double scalar)
Multiplies the diagonal of a matrix on the OpenCL device with the specified scalar.
Eigen::Block< Eigen::MatrixXd > Block_