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>
17 class mdivide_left_spd_alloc :
public chainable_alloc {
19 virtual ~mdivide_left_spd_alloc() {}
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>
26 class mdivide_left_spd_vv_vari :
public vari {
33 mdivide_left_spd_alloc<R1, C1, R2, C2> *
_alloc;
35 mdivide_left_spd_vv_vari(
const Eigen::Matrix<var, R1, C1> &A,
36 const Eigen::Matrix<var, R2, C2> &B)
40 _variRefA(reinterpret_cast<vari**>
42 .alloc(sizeof(vari*) * A.
rows() * A.
cols()))),
43 _variRefB(reinterpret_cast<vari**>
45 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
46 _variRefC(reinterpret_cast<vari**>
48 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
49 _alloc(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
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);
86 virtual void chain() {
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>
113 class mdivide_left_spd_dv_vari :
public vari {
119 mdivide_left_spd_alloc<R1, C1, R2, C2> *
_alloc;
121 mdivide_left_spd_dv_vari(
const Eigen::Matrix<double, R1, C1> &A,
122 const Eigen::Matrix<var, R2, C2> &B)
126 _variRefB(reinterpret_cast<vari**>
128 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
129 _variRefC(reinterpret_cast<vari**>
131 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
132 _alloc(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
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);
158 virtual void chain() {
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>
178 class mdivide_left_spd_vd_vari :
public vari {
184 mdivide_left_spd_alloc<R1, C1, R2, C2> *
_alloc;
186 mdivide_left_spd_vd_vari(
const Eigen::Matrix<var, R1, C1> &A,
187 const Eigen::Matrix<double, R2, C2> &B)
191 _variRefA(reinterpret_cast<vari**>
193 .alloc(sizeof(vari*) * A.
rows() * A.
cols()))),
194 _variRefC(reinterpret_cast<vari**>
196 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
197 _alloc(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
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);
224 virtual void chain() {
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 Eigen::Matrix<var, R1, C2>
249 const Eigen::Matrix<var, R2, C2> &b) {
250 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
261 mdivide_left_spd_vv_vari<R1, C1, R2, C2> *baseVari
262 =
new mdivide_left_spd_vv_vari<R1, C1, R2, C2>(A, b);
265 for (
size_type j = 0; j < res.cols(); j++)
266 for (
size_type i = 0; i < res.rows(); i++)
267 res(i, j).vi_ = baseVari->_variRefC[pos++];
272 template <
int R1,
int C1,
int R2,
int C2>
274 Eigen::Matrix<var, R1, C2>
276 const Eigen::Matrix<double, R2, C2> &b) {
277 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
288 mdivide_left_spd_vd_vari<R1, C1, R2, C2> *baseVari
289 =
new mdivide_left_spd_vd_vari<R1, C1, R2, C2>(A, b);
292 for (
size_type j = 0; j < res.cols(); j++)
293 for (
size_type i = 0; i < res.rows(); i++)
294 res(i, j).vi_ = baseVari->_variRefC[pos++];
299 template <
int R1,
int C1,
int R2,
int C2>
301 Eigen::Matrix<var, R1, C2>
303 const Eigen::Matrix<var, R2, C2> &b) {
304 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
315 mdivide_left_spd_dv_vari<R1, C1, R2, C2> *baseVari
316 =
new mdivide_left_spd_dv_vari<R1, C1, R2, C2>(A, b);
319 for (
size_type j = 0; j < res.cols(); j++)
320 for (
size_type i = 0; i < res.rows(); i++)
321 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.
Eigen::LLT< Eigen::Matrix< double, R1, C1 > > _llt
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.
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.
mdivide_left_spd_alloc< R1, C1, R2, C2 > * _alloc
int cols(const Eigen::Matrix< T, R, C > &m)
Return the number of columns in the specified matrix, vector, or row vector.
bool check_square(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Return true if the specified matrix is square.
AutodiffStackStorage< vari, chainable_alloc > ChainableStack
Eigen::Matrix< double, R2, C2 > C_