Stan Math Library  2.20.0
reverse mode automatic differentiation
mdivide_left_ldlt.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_LDLT_HPP
2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_LDLT_HPP
3 
4 #include <stan/math/rev/meta.hpp>
5 #include <stan/math/rev/core.hpp>
10 
11 namespace stan {
12 namespace math {
13 namespace internal {
14 template <int R1, int C1, int R2, int C2>
16  public:
18 
24  boost::shared_ptr<Eigen::LDLT<Eigen::Matrix<double, R1, C1> > > ldltP_;
25  Eigen::Matrix<double, R2, C2> C_;
26 };
27 
38 template <int R1, int C1, int R2, int C2>
40  public:
41  int M_; // A.rows() = A.cols() = B.rows()
42  int N_; // B.cols()
47 
49  const Eigen::Matrix<var, R2, C2> &B)
50  : vari(0.0),
51  M_(A.rows()),
52  N_(B.cols()),
53  variRefB_(reinterpret_cast<vari **>(
54  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
55  * B.cols()))),
56  variRefC_(reinterpret_cast<vari **>(
57  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
58  * B.cols()))),
59  alloc_(new mdivide_left_ldlt_alloc<R1, C1, R2, C2>()),
60  alloc_ldlt_(A.alloc_) {
61  int pos = 0;
62  alloc_->C_.resize(M_, N_);
63  for (int j = 0; j < N_; j++) {
64  for (int i = 0; i < M_; i++) {
65  variRefB_[pos] = B(i, j).vi_;
66  alloc_->C_(i, j) = B(i, j).val();
67  pos++;
68  }
69  }
70 
71  alloc_ldlt_->ldlt_.solveInPlace(alloc_->C_);
72 
73  pos = 0;
74  for (int j = 0; j < N_; j++) {
75  for (int i = 0; i < M_; i++) {
76  variRefC_[pos] = new vari(alloc_->C_(i, j), false);
77  pos++;
78  }
79  }
80  }
81 
82  virtual void chain() {
83  Eigen::Matrix<double, R1, C1> adjA(M_, M_);
84  Eigen::Matrix<double, R2, C2> adjB(M_, N_);
85 
86  int pos = 0;
87  for (int j = 0; j < N_; j++)
88  for (int i = 0; i < M_; i++)
89  adjB(i, j) = variRefC_[pos++]->adj_;
90 
91  alloc_ldlt_->ldlt_.solveInPlace(adjB);
92  adjA.noalias() = -adjB * alloc_->C_.transpose();
93 
94  for (int j = 0; j < M_; j++)
95  for (int i = 0; i < M_; i++)
96  alloc_ldlt_->variA_(i, j)->adj_ += adjA(i, j);
97 
98  pos = 0;
99  for (int j = 0; j < N_; j++)
100  for (int i = 0; i < M_; i++)
101  variRefB_[pos++]->adj_ += adjB(i, j);
102  }
103 };
104 
115 template <int R1, int C1, int R2, int C2>
117  public:
118  int M_; // A.rows() = A.cols() = B.rows()
119  int N_; // B.cols()
123 
125  const Eigen::Matrix<var, R2, C2> &B)
126  : vari(0.0),
127  M_(A.rows()),
128  N_(B.cols()),
129  variRefB_(reinterpret_cast<vari **>(
130  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
131  * B.cols()))),
132  variRefC_(reinterpret_cast<vari **>(
133  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
134  * B.cols()))),
135  alloc_(new mdivide_left_ldlt_alloc<R1, C1, R2, C2>()) {
136  using Eigen::Map;
137  using Eigen::Matrix;
138 
139  int pos = 0;
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();
145  pos++;
146  }
147  }
148 
149  alloc_->ldltP_ = A.ldltP_;
150  alloc_->ldltP_->solveInPlace(alloc_->C_);
151 
152  pos = 0;
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);
156  pos++;
157  }
158  }
159  }
160 
161  virtual void chain() {
162  Eigen::Matrix<double, R2, C2> adjB(M_, N_);
163 
164  int pos = 0;
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_;
168 
169  alloc_->ldltP_->solveInPlace(adjB);
170 
171  pos = 0;
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);
175  }
176 };
177 
188 template <int R1, int C1, int R2, int C2>
190  public:
191  int M_; // A.rows() = A.cols() = B.rows()
192  int N_; // B.cols()
196 
198  const Eigen::Matrix<double, R2, C2> &B)
199  : vari(0.0),
200  M_(A.rows()),
201  N_(B.cols()),
202  variRefC_(reinterpret_cast<vari **>(
203  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
204  * B.cols()))),
205  alloc_(new mdivide_left_ldlt_alloc<R1, C1, R2, C2>()),
206  alloc_ldlt_(A.alloc_) {
207  alloc_->C_ = B;
208  alloc_ldlt_->ldlt_.solveInPlace(alloc_->C_);
209 
210  int pos = 0;
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);
214  pos++;
215  }
216  }
217  }
218 
219  virtual void chain() {
220  Eigen::Matrix<double, R1, C1> adjA(M_, M_);
221  Eigen::Matrix<double, R1, C2> adjC(M_, N_);
222 
223  int pos = 0;
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_;
227 
228  adjA = -alloc_ldlt_->ldlt_.solve(adjC * alloc_->C_.transpose());
229 
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);
233  }
234 };
235 } // namespace internal
236 
244 template <int R1, int C1, int R2, int C2>
245 inline Eigen::Matrix<var, R1, C2> mdivide_left_ldlt(
246  const LDLT_factor<var, R1, C1> &A, const Eigen::Matrix<var, R2, C2> &b) {
247  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
248 
249  check_multiplicable("mdivide_left_ldlt", "A", A, "b", b);
250 
253 
254  int pos = 0;
255  for (int j = 0; j < res.cols(); j++)
256  for (int i = 0; i < res.rows(); i++)
257  res(i, j).vi_ = baseVari->variRefC_[pos++];
258 
259  return res;
260 }
261 
269 template <int R1, int C1, int R2, int C2>
270 inline Eigen::Matrix<var, R1, C2> mdivide_left_ldlt(
271  const LDLT_factor<var, R1, C1> &A, const Eigen::Matrix<double, R2, C2> &b) {
272  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
273 
274  check_multiplicable("mdivide_left_ldlt", "A", A, "b", b);
275 
278 
279  int pos = 0;
280  for (int j = 0; j < res.cols(); j++)
281  for (int i = 0; i < res.rows(); i++)
282  res(i, j).vi_ = baseVari->variRefC_[pos++];
283 
284  return res;
285 }
286 
294 template <int R1, int C1, int R2, int C2>
295 inline Eigen::Matrix<var, R1, C2> mdivide_left_ldlt(
296  const LDLT_factor<double, R1, C1> &A, const Eigen::Matrix<var, R2, C2> &b) {
297  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
298 
299  check_multiplicable("mdivide_left_ldlt", "A", A, "b", b);
300 
303 
304  int pos = 0;
305  for (int j = 0; j < res.cols(); j++)
306  for (int i = 0; i < res.rows(); i++)
307  res(i, j).vi_ = baseVari->variRefC_[pos++];
308 
309  return res;
310 }
311 
312 } // namespace math
313 } // namespace stan
314 #endif
int rows(const Eigen::Matrix< T, R, C > &m)
Return the number of rows in the specified matrix, vector, or row vector.
Definition: rows.hpp:20
Eigen::LDLT< Eigen::Matrix< double, R, C > > ldlt_
Definition: LDLT_alloc.hpp:53
Eigen::Matrix< vari *, R, C > variA_
Definition: LDLT_alloc.hpp:54
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
mdivide_left_ldlt_alloc< R1, C1, R2, C2 > * alloc_
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.
The variable implementation base class.
Definition: vari.hpp:30
mdivide_left_ldlt_alloc< R1, C1, R2, C2 > * alloc_
LDLT_factor is a thin wrapper on Eigen::LDLT to allow for reusing factorizations and efficient autodi...
Definition: LDLT_factor.hpp:63
boost::shared_ptr< ldlt_t > ldltP_
mdivide_left_ldlt_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.
Definition: cols.hpp:20
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<double>.
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.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
mdivide_left_ldlt_vv_vari(const LDLT_factor< 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. ...
mdivide_left_ldlt_dv_vari(const LDLT_factor< double, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)
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...
Definition: vari.hpp:44
The vari for mdivide_left_ldlt(A, b) which handles the chain() call for all elements of the result...
This struct always provides access to the autodiff stack using the singleton pattern.
mdivide_left_ldlt_vd_vari(const LDLT_factor< var, R1, C1 > &A, const Eigen::Matrix< double, R2, C2 > &B)
The vari for mdivide_left_ldlt(A, b) which handles the chain() call for all elements of the result...
The vari for mdivide_left_ldlt(A, b) which handles the chain() call for all elements of the result...

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