1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_SPD_HPP 2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_SPD_HPP 16 template <
int R1,
int C1,
int R2,
int C2>
21 Eigen::LLT<Eigen::Matrix<double, R1, C1> >
llt_;
22 Eigen::Matrix<double, R2, C2>
C_;
25 template <
int R1,
int C1,
int R2,
int C2>
36 const Eigen::Matrix<var, R2, C2> &B)
40 variRefA_(reinterpret_cast<
vari **>(
43 variRefB_(reinterpret_cast<
vari **>(
46 variRefC_(reinterpret_cast<
vari **>(
53 Matrix<double, R1, C1> Ad(A.rows(), A.cols());
58 variRefA_[pos] = A(i, j).vi_;
59 Ad(i, j) = A(i, j).val();
65 alloc_->
C_.resize(M_, N_);
68 variRefB_[pos] = B(i, j).vi_;
69 alloc_->
C_(i, j) = B(i, j).val();
74 alloc_->
llt_ = Ad.llt();
75 alloc_->
llt_.solveInPlace(alloc_->
C_);
80 variRefC_[pos] =
new vari(alloc_->
C_(i, j),
false);
89 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
90 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
95 adjB(i, j) = variRefC_[pos++]->
adj_;
97 alloc_->
llt_.solveInPlace(adjB);
98 adjA.noalias() = -adjB * alloc_->
C_.transpose();
103 variRefA_[pos++]->adj_ += adjA(i, j);
108 variRefB_[pos++]->adj_ += adjB(i, j);
112 template <
int R1,
int C1,
int R2,
int C2>
122 const Eigen::Matrix<var, R2, C2> &B)
126 variRefB_(reinterpret_cast<
vari **>(
129 variRefC_(reinterpret_cast<
vari **>(
137 alloc_->
C_.resize(M_, N_);
140 variRefB_[pos] = B(i, j).vi_;
141 alloc_->
C_(i, j) = B(i, j).val();
146 alloc_->
llt_ = A.llt();
147 alloc_->
llt_.solveInPlace(alloc_->
C_);
152 variRefC_[pos] =
new vari(alloc_->
C_(i, j),
false);
161 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
164 for (
size_type j = 0; j < adjB.cols(); j++)
165 for (
size_type i = 0; i < adjB.rows(); i++)
166 adjB(i, j) = variRefC_[pos++]->
adj_;
168 alloc_->
llt_.solveInPlace(adjB);
171 for (
size_type j = 0; j < adjB.cols(); j++)
172 for (
size_type i = 0; i < adjB.rows(); i++)
173 variRefB_[pos++]->adj_ += adjB(i, j);
177 template <
int R1,
int C1,
int R2,
int C2>
187 const Eigen::Matrix<double, R2, C2> &B)
191 variRefA_(reinterpret_cast<
vari **>(
194 variRefC_(reinterpret_cast<
vari **>(
201 Matrix<double, R1, C1> Ad(A.rows(), A.cols());
206 variRefA_[pos] = A(i, j).vi_;
207 Ad(i, j) = A(i, j).val();
212 alloc_->
llt_ = Ad.llt();
213 alloc_->
C_ = alloc_->
llt_.solve(B);
218 variRefC_[pos] =
new vari(alloc_->
C_(i, j),
false);
227 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
228 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
231 for (
size_type j = 0; j < adjC.cols(); j++)
232 for (
size_type i = 0; i < adjC.rows(); i++)
233 adjC(i, j) = variRefC_[pos++]->
adj_;
235 adjA = -alloc_->
llt_.solve(adjC * alloc_->
C_.transpose());
238 for (
size_type j = 0; j < adjA.cols(); j++)
239 for (
size_type i = 0; i < adjA.rows(); i++)
240 variRefA_[pos++]->adj_ += adjA(i, j);
245 template <
int R1,
int C1,
int R2,
int C2>
247 const Eigen::Matrix<var, R1, C1> &A,
const Eigen::Matrix<var, R2, C2> &b) {
248 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
261 for (
size_type j = 0; j < res.cols(); j++)
262 for (
size_type i = 0; i < res.rows(); i++)
263 res(i, j).vi_ = baseVari->
variRefC_[pos++];
268 template <
int R1,
int C1,
int R2,
int C2>
270 const Eigen::Matrix<var, R1, C1> &A,
271 const Eigen::Matrix<double, R2, C2> &b) {
272 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
285 for (
size_type j = 0; j < res.cols(); j++)
286 for (
size_type i = 0; i < res.rows(); i++)
287 res(i, j).vi_ = baseVari->
variRefC_[pos++];
292 template <
int R1,
int C1,
int R2,
int C2>
294 const Eigen::Matrix<double, R1, C1> &A,
295 const Eigen::Matrix<var, R2, C2> &b) {
296 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
309 for (
size_type j = 0; j < res.cols(); j++)
310 for (
size_type i = 0; i < res.rows(); i++)
311 res(i, j).vi_ = baseVari->
variRefC_[pos++];
int rows(const Eigen::Matrix< T, R, C > &m)
Return the number of rows in the specified matrix, vector, or row vector.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
mdivide_left_spd_alloc< R1, C1, R2, C2 > * alloc_
virtual ~mdivide_left_spd_alloc()
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.
mdivide_left_spd_vd_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< double, R2, C2 > &B)
Eigen::Matrix< typename boost::math::tools::promote_args< T1, T2 >::type, R1, C2 > mdivide_left_spd(const Eigen::Matrix< T1, R1, C1 > &A, const Eigen::Matrix< T2, R2, C2 > &b)
Returns the solution of the system Ax=b where A is symmetric positive definite.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double e.
mdivide_left_spd_vv_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Eigen::LLT< Eigen::Matrix< double, R1, C1 > > llt_
Eigen::Matrix< double, R2, C2 > C_
int cols(const Eigen::Matrix< T, R, C > &m)
Return the number of columns in the specified matrix, vector, or row vector.
void check_multiplicable(const char *function, const char *name1, const T1 &y1, const char *name2, const T2 &y2)
Check if the matrices can be multiplied.
mdivide_left_spd_alloc< R1, C1, R2, C2 > * alloc_
mdivide_left_spd_alloc< R1, C1, R2, C2 > * alloc_
A chainable_alloc is an object which is constructed and destructed normally but the memory lifespan i...
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
This struct always provides access to the autodiff stack using the singleton pattern.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
mdivide_left_spd_dv_vari(const Eigen::Matrix< double, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)