1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_LDLT_HPP
2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_LDLT_HPP
13 template <
int R1,
int C1,
int R2,
int C2>
14 class mdivide_left_ldlt_alloc :
public chainable_alloc {
16 virtual ~mdivide_left_ldlt_alloc() {}
23 boost::shared_ptr<Eigen::LDLT<Eigen::Matrix<double, R1, C1> > >
ldltP_;
24 Eigen::Matrix<double, R2, C2>
C_;
37 template <
int R1,
int C1,
int R2,
int C2>
38 class mdivide_left_ldlt_vv_vari :
public vari {
44 mdivide_left_ldlt_alloc<R1, C1, R2, C2> *
alloc_;
47 mdivide_left_ldlt_vv_vari(
const LDLT_factor<var, R1, C1> &A,
48 const Eigen::Matrix<var, R2, C2> &B)
52 variRefB_(reinterpret_cast<vari**>
54 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
55 variRefC_(reinterpret_cast<vari**>
57 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
58 alloc_(new mdivide_left_ldlt_alloc<R1, C1, R2, C2>()),
59 alloc_ldlt_(A.alloc_) {
61 alloc_->C_.resize(M_, N_);
62 for (
int j = 0; j <
N_; j++) {
63 for (
int i = 0; i <
M_; i++) {
64 variRefB_[pos] = B(i, j).vi_;
65 alloc_->C_(i, j) = B(i, j).val();
70 alloc_ldlt_->ldlt_.solveInPlace(alloc_->C_);
73 for (
int j = 0; j <
N_; j++) {
74 for (
int i = 0; i <
M_; i++) {
75 variRefC_[pos] =
new vari(alloc_->C_(i, j),
false);
81 virtual void chain() {
82 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
83 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
86 for (
int j = 0; j <
N_; j++)
87 for (
int i = 0; i <
M_; i++)
88 adjB(i, j) = variRefC_[pos++]->adj_;
90 alloc_ldlt_->ldlt_.solveInPlace(adjB);
91 adjA.noalias() = -adjB * alloc_->C_.transpose();
93 for (
int j = 0; j <
M_; j++)
94 for (
int i = 0; i <
M_; i++)
95 alloc_ldlt_->variA_(i, j)->adj_ += adjA(i, j);
98 for (
int j = 0; j <
N_; j++)
99 for (
int i = 0; i <
M_; i++)
100 variRefB_[pos++]->adj_ += adjB(i, j);
114 template <
int R1,
int C1,
int R2,
int C2>
115 class mdivide_left_ldlt_dv_vari :
public vari {
121 mdivide_left_ldlt_alloc<R1, C1, R2, C2> *
alloc_;
123 mdivide_left_ldlt_dv_vari(
const LDLT_factor<double, R1, C1>
125 const Eigen::Matrix<var, R2, C2> &B)
129 variRefB_(reinterpret_cast<vari**>
131 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
132 variRefC_(reinterpret_cast<vari**>
134 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
135 alloc_(new mdivide_left_ldlt_alloc<R1, C1, R2, C2>()) {
140 alloc_->C_.resize(M_, N_);
141 for (
int j = 0; j <
N_; j++) {
142 for (
int i = 0; i <
M_; i++) {
143 variRefB_[pos] = B(i, j).vi_;
144 alloc_->C_(i, j) = B(i, j).val();
149 alloc_->ldltP_ = A.ldltP_;
150 alloc_->ldltP_->solveInPlace(alloc_->C_);
153 for (
int j = 0; j <
N_; j++) {
154 for (
int i = 0; i <
M_; i++) {
155 variRefC_[pos] =
new vari(alloc_->C_(i, j),
false);
161 virtual void chain() {
162 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
165 for (
int j = 0; j < adjB.cols(); j++)
166 for (
int i = 0; i < adjB.rows(); i++)
167 adjB(i, j) = variRefC_[pos++]->adj_;
169 alloc_->ldltP_->solveInPlace(adjB);
172 for (
int j = 0; j < adjB.cols(); j++)
173 for (
int i = 0; i < adjB.rows(); i++)
174 variRefB_[pos++]->adj_ += adjB(i, j);
188 template <
int R1,
int C1,
int R2,
int C2>
189 class mdivide_left_ldlt_vd_vari :
public vari {
194 mdivide_left_ldlt_alloc<R1, C1, R2, C2> *
alloc_;
197 mdivide_left_ldlt_vd_vari(
const LDLT_factor<var, R1, C1> &A,
198 const Eigen::Matrix<double, R2, C2> &B)
202 variRefC_(reinterpret_cast<vari**>
204 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
205 alloc_(new mdivide_left_ldlt_alloc<R1, C1, R2, C2>()),
206 alloc_ldlt_(A.alloc_) {
208 alloc_ldlt_->ldlt_.solveInPlace(alloc_->C_);
211 for (
int j = 0; j <
N_; j++) {
212 for (
int i = 0; i <
M_; i++) {
213 variRefC_[pos] =
new vari(alloc_->C_(i, j),
false);
219 virtual void chain() {
220 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
221 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
224 for (
int j = 0; j < adjC.cols(); j++)
225 for (
int i = 0; i < adjC.rows(); i++)
226 adjC(i, j) = variRefC_[pos++]->adj_;
228 adjA = -alloc_ldlt_->ldlt_.solve(adjC*alloc_->C_.transpose());
230 for (
int j = 0; j < adjA.cols(); j++)
231 for (
int i = 0; i < adjA.rows(); i++)
232 alloc_ldlt_->variA_(i, j)->adj_ += adjA(i, j);
244 template <
int R1,
int C1,
int R2,
int C2>
245 inline Eigen::Matrix<var, R1, C2>
247 const Eigen::Matrix<var, R2, C2> &b) {
248 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
254 mdivide_left_ldlt_vv_vari<R1, C1, R2, C2> *baseVari
255 =
new mdivide_left_ldlt_vv_vari<R1, C1, R2, C2>(A, b);
258 for (
int j = 0; j < res.cols(); j++)
259 for (
int i = 0; i < res.rows(); i++)
260 res(i, j).vi_ = baseVari->variRefC_[pos++];
272 template <
int R1,
int C1,
int R2,
int C2>
273 inline Eigen::Matrix<var, R1, C2>
275 const Eigen::Matrix<double, R2, C2> &b) {
276 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
282 mdivide_left_ldlt_vd_vari<R1, C1, R2, C2> *baseVari
283 =
new mdivide_left_ldlt_vd_vari<R1, C1, R2, C2>(A, b);
286 for (
int j = 0; j < res.cols(); j++)
287 for (
int i = 0; i < res.rows(); i++)
288 res(i, j).vi_ = baseVari->variRefC_[pos++];
300 template <
int R1,
int C1,
int R2,
int C2>
301 inline Eigen::Matrix<var, R1, C2>
303 const Eigen::Matrix<var, R2, C2> &b) {
304 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
310 mdivide_left_ldlt_dv_vari<R1, C1, R2, C2> *baseVari
311 =
new mdivide_left_ldlt_dv_vari<R1, C1, R2, C2>(A, b);
314 for (
int j = 0; j < res.cols(); j++)
315 for (
int i = 0; i < res.rows(); i++)
316 res(i, j).vi_ = baseVari->variRefC_[pos++];
boost::shared_ptr< Eigen::LDLT< Eigen::Matrix< double, R1, C1 > > > ldltP_
This share_ptr is used to prevent copying the LDLT factorizations for mdivide_left_ldlt(ldltA, b) when ldltA is a LDLT_factor.
int rows(const Eigen::Matrix< T, R, C > &m)
Return the number of rows in the specified matrix, vector, or row vector.
Eigen::Matrix< fvar< T2 >, R1, C2 > mdivide_left_ldlt(const LDLT_factor< double, R1, C1 > &A, const Eigen::Matrix< fvar< T2 >, R2, C2 > &b)
Returns the solution of the system Ax=b given an LDLT_factor of A.
mdivide_left_ldlt_alloc< R1, C1, R2, C2 > * alloc_
bool check_multiplicable(const char *function, const char *name1, const T1 &y1, const char *name2, const T2 &y2)
Return true if the matrices can be multiplied.
const LDLT_alloc< R1, C1 > * alloc_ldlt_
int cols(const Eigen::Matrix< T, R, C > &m)
Return the number of columns in the specified matrix, vector, or row vector.
Eigen::Matrix< double, R2, C2 > C_
AutodiffStackStorage< vari, chainable_alloc > ChainableStack