Stan Math Library  2.6.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros
mdivide_left_tri.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_TRI_HPP
2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_TRI_HPP
3 
8 #include <stan/math/rev/core.hpp>
10 #include <vector>
11 
12 namespace stan {
13  namespace math {
14 
15  namespace {
16  template <int TriView, int R1, int C1, int R2, int C2>
17  class mdivide_left_tri_vv_vari : public vari {
18  public:
19  int M_; // A.rows() = A.cols() = B.rows()
20  int N_; // B.cols()
21  double* A_;
22  double* C_;
23  vari** _variRefA;
24  vari** _variRefB;
25  vari** _variRefC;
26 
27  mdivide_left_tri_vv_vari(const Eigen::Matrix<var, R1, C1> &A,
28  const Eigen::Matrix<var, R2, C2> &B)
29  : vari(0.0),
30  M_(A.rows()),
31  N_(B.cols()),
32  A_(reinterpret_cast<double*>
33  (stan::math::ChainableStack::memalloc_
34  .alloc(sizeof(double) * A.rows() * A.cols()))),
35  C_(reinterpret_cast<double*>
36  (stan::math::ChainableStack::memalloc_
37  .alloc(sizeof(double) * B.rows() * B.cols()))),
38  _variRefA(reinterpret_cast<vari**>
39  (stan::math::ChainableStack::memalloc_
40  .alloc(sizeof(vari*) * A.rows() * (A.rows() + 1) / 2))),
41  _variRefB(reinterpret_cast<vari**>
42  (stan::math::ChainableStack::memalloc_
43  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
44  _variRefC(reinterpret_cast<vari**>
45  (stan::math::ChainableStack::memalloc_
46  .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
47  using Eigen::Matrix;
48  using Eigen::Map;
49 
50  size_t pos = 0;
51  if (TriView == Eigen::Lower) {
52  for (size_type j = 0; j < M_; j++)
53  for (size_type i = j; i < M_; i++)
54  _variRefA[pos++] = A(i, j).vi_;
55  } else if (TriView == Eigen::Upper) {
56  for (size_type j = 0; j < M_; j++)
57  for (size_type i = 0; i < j+1; i++)
58  _variRefA[pos++] = A(i, j).vi_;
59  }
60 
61  pos = 0;
62  for (size_type j = 0; j < M_; j++) {
63  for (size_type i = 0; i < M_; i++) {
64  A_[pos++] = A(i, j).val();
65  }
66  }
67 
68  pos = 0;
69  for (size_type j = 0; j < N_; j++) {
70  for (size_type i = 0; i < M_; i++) {
71  _variRefB[pos] = B(i, j).vi_;
72  C_[pos++] = B(i, j).val();
73  }
74  }
75 
76  Matrix<double, R1, C2> C(M_, N_);
77  C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
78 
79  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
80  .template triangularView<TriView>().solve(C);
81 
82  pos = 0;
83  for (size_type j = 0; j < N_; j++) {
84  for (size_type i = 0; i < M_; i++) {
85  C_[pos] = C(i, j);
86  _variRefC[pos] = new vari(C_[pos], false);
87  pos++;
88  }
89  }
90  }
91 
92  virtual void chain() {
93  using Eigen::Matrix;
94  using Eigen::Map;
95  Matrix<double, R1, C1> adjA(M_, M_);
96  Matrix<double, R2, C2> adjB(M_, N_);
97  Matrix<double, R1, C2> adjC(M_, N_);
98 
99  size_t pos = 0;
100  for (size_type j = 0; j < adjC.cols(); j++)
101  for (size_type i = 0; i < adjC.rows(); i++)
102  adjC(i, j) = _variRefC[pos++]->adj_;
103 
104  adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
105  .template triangularView<TriView>().transpose().solve(adjC);
106  adjA.noalias() = -adjB
107  * Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose();
108 
109  pos = 0;
110  if (TriView == Eigen::Lower) {
111  for (size_type j = 0; j < adjA.cols(); j++)
112  for (size_type i = j; i < adjA.rows(); i++)
113  _variRefA[pos++]->adj_ += adjA(i, j);
114  } else if (TriView == Eigen::Upper) {
115  for (size_type j = 0; j < adjA.cols(); j++)
116  for (size_type i = 0; i < j+1; i++)
117  _variRefA[pos++]->adj_ += adjA(i, j);
118  }
119 
120  pos = 0;
121  for (size_type j = 0; j < adjB.cols(); j++)
122  for (size_type i = 0; i < adjB.rows(); i++)
123  _variRefB[pos++]->adj_ += adjB(i, j);
124  }
125  };
126 
127  template <int TriView, int R1, int C1, int R2, int C2>
128  class mdivide_left_tri_dv_vari : public vari {
129  public:
130  int M_; // A.rows() = A.cols() = B.rows()
131  int N_; // B.cols()
132  double* A_;
133  double* C_;
134  vari** _variRefB;
135  vari** _variRefC;
136 
137  mdivide_left_tri_dv_vari(const Eigen::Matrix<double, R1, C1> &A,
138  const Eigen::Matrix<var, R2, C2> &B)
139  : vari(0.0),
140  M_(A.rows()),
141  N_(B.cols()),
142  A_(reinterpret_cast<double*>
143  (stan::math::ChainableStack::memalloc_
144  .alloc(sizeof(double) * A.rows() * A.cols()))),
145  C_(reinterpret_cast<double*>
146  (stan::math::ChainableStack::memalloc_
147  .alloc(sizeof(double) * B.rows() * B.cols()))),
148  _variRefB(reinterpret_cast<vari**>
149  (stan::math::ChainableStack::memalloc_
150  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
151  _variRefC(reinterpret_cast<vari**>
152  (stan::math::ChainableStack::memalloc_
153  .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
154  using Eigen::Matrix;
155  using Eigen::Map;
156 
157  size_t pos = 0;
158  for (size_type j = 0; j < M_; j++) {
159  for (size_type i = 0; i < M_; i++) {
160  A_[pos++] = A(i, j);
161  }
162  }
163 
164  pos = 0;
165  for (size_type j = 0; j < N_; j++) {
166  for (size_type i = 0; i < M_; i++) {
167  _variRefB[pos] = B(i, j).vi_;
168  C_[pos++] = B(i, j).val();
169  }
170  }
171 
172  Matrix<double, R1, C2> C(M_, N_);
173  C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
174 
175  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
176  .template triangularView<TriView>().solve(C);
177 
178  pos = 0;
179  for (size_type j = 0; j < N_; j++) {
180  for (size_type i = 0; i < M_; i++) {
181  C_[pos] = C(i, j);
182  _variRefC[pos] = new vari(C_[pos], false);
183  pos++;
184  }
185  }
186  }
187 
188  virtual void chain() {
189  using Eigen::Matrix;
190  using Eigen::Map;
191  Matrix<double, R2, C2> adjB(M_, N_);
192  Matrix<double, R1, C2> adjC(M_, N_);
193 
194  size_t pos = 0;
195  for (size_type j = 0; j < adjC.cols(); j++)
196  for (size_type i = 0; i < adjC.rows(); i++)
197  adjC(i, j) = _variRefC[pos++]->adj_;
198 
199  adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
200  .template triangularView<TriView>().transpose().solve(adjC);
201 
202  pos = 0;
203  for (size_type j = 0; j < adjB.cols(); j++)
204  for (size_type i = 0; i < adjB.rows(); i++)
205  _variRefB[pos++]->adj_ += adjB(i, j);
206  }
207  };
208 
209  template <int TriView, int R1, int C1, int R2, int C2>
210  class mdivide_left_tri_vd_vari : public vari {
211  public:
212  int M_; // A.rows() = A.cols() = B.rows()
213  int N_; // B.cols()
214  double* A_;
215  double* C_;
216  vari** _variRefA;
217  vari** _variRefC;
218 
219  mdivide_left_tri_vd_vari(const Eigen::Matrix<var, R1, C1> &A,
220  const Eigen::Matrix<double, R2, C2> &B)
221  : vari(0.0),
222  M_(A.rows()),
223  N_(B.cols()),
224  A_(reinterpret_cast<double*>
225  (stan::math::ChainableStack::memalloc_
226  .alloc(sizeof(double) * A.rows() * A.cols()))),
227  C_(reinterpret_cast<double*>
228  (stan::math::ChainableStack::memalloc_
229  .alloc(sizeof(double) * B.rows() * B.cols()))),
230  _variRefA(reinterpret_cast<vari**>
231  (stan::math::ChainableStack::memalloc_
232  .alloc(sizeof(vari*) * A.rows() * (A.rows() + 1) / 2))),
233  _variRefC(reinterpret_cast<vari**>
234  (stan::math::ChainableStack::memalloc_
235  .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
236  using Eigen::Matrix;
237  using Eigen::Map;
238 
239  size_t pos = 0;
240  if (TriView == Eigen::Lower) {
241  for (size_type j = 0; j < M_; j++)
242  for (size_type i = j; i < M_; i++)
243  _variRefA[pos++] = A(i, j).vi_;
244  } else if (TriView == Eigen::Upper) {
245  for (size_type j = 0; j < M_; j++)
246  for (size_type i = 0; i < j+1; i++)
247  _variRefA[pos++] = A(i, j).vi_;
248  }
249 
250  pos = 0;
251  for (size_type j = 0; j < M_; j++) {
252  for (size_type i = 0; i < M_; i++) {
253  A_[pos++] = A(i, j).val();
254  }
255  }
256 
257  Matrix<double, R1, C2> C(M_, N_);
258  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
259  .template triangularView<TriView>().solve(B);
260 
261  pos = 0;
262  for (size_type j = 0; j < N_; j++) {
263  for (size_type i = 0; i < M_; i++) {
264  C_[pos] = C(i, j);
265  _variRefC[pos] = new vari(C_[pos], false);
266  pos++;
267  }
268  }
269  }
270 
271  virtual void chain() {
272  using Eigen::Matrix;
273  using Eigen::Map;
274  Matrix<double, R1, C1> adjA(M_, M_);
275  Matrix<double, R1, C2> adjC(M_, N_);
276 
277  size_t pos = 0;
278  for (size_type j = 0; j < adjC.cols(); j++)
279  for (size_type i = 0; i < adjC.rows(); i++)
280  adjC(i, j) = _variRefC[pos++]->adj_;
281 
282  adjA.noalias() = -Map<Matrix<double, R1, C1> >(A_, M_, M_)
283  .template triangularView<TriView>()
284  .transpose().solve(adjC * Map<Matrix<double, R1, C2> >(C_, M_, N_)
285  .transpose());
286 
287  pos = 0;
288  if (TriView == Eigen::Lower) {
289  for (size_type j = 0; j < adjA.cols(); j++)
290  for (size_type i = j; i < adjA.rows(); i++)
291  _variRefA[pos++]->adj_ += adjA(i, j);
292  } else if (TriView == Eigen::Upper) {
293  for (size_type j = 0; j < adjA.cols(); j++)
294  for (size_type i = 0; i < j+1; i++)
295  _variRefA[pos++]->adj_ += adjA(i, j);
296  }
297  }
298  };
299  }
300 
301  template <int TriView, int R1, int C1, int R2, int C2>
302  inline
303  Eigen::Matrix<var, R1, C2>
304  mdivide_left_tri(const Eigen::Matrix<var, R1, C1> &A,
305  const Eigen::Matrix<var, R2, C2> &b) {
306  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
307 
308  stan::math::check_square("mdivide_left_tri", "A", A);
309  stan::math::check_multiplicable("mdivide_left_tri",
310  "A", A,
311  "b", b);
312 
313  // NOTE: this is not a memory leak, this vari is used in the
314  // expression graph to evaluate the adjoint, but is not needed
315  // for the returned matrix. Memory will be cleaned up with the
316  // arena allocator.
317  mdivide_left_tri_vv_vari<TriView, R1, C1, R2, C2> *baseVari
318  = new mdivide_left_tri_vv_vari<TriView, R1, C1, R2, C2>(A, b);
319 
320  size_t pos = 0;
321  for (size_type j = 0; j < res.cols(); j++)
322  for (size_type i = 0; i < res.rows(); i++)
323  res(i, j).vi_ = baseVari->_variRefC[pos++];
324 
325  return res;
326  }
327  template <int TriView, int R1, int C1, int R2, int C2>
328  inline
329  Eigen::Matrix<var, R1, C2>
330  mdivide_left_tri(const Eigen::Matrix<double, R1, C1> &A,
331  const Eigen::Matrix<var, R2, C2> &b) {
332  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
333 
334  stan::math::check_square("mdivide_left_tri", "A", A);
335  stan::math::check_multiplicable("mdivide_left_tri",
336  "A", A,
337  "b", b);
338 
339  // NOTE: this is not a memory leak, this vari is used in the
340  // expression graph to evaluate the adjoint, but is not needed
341  // for the returned matrix. Memory will be cleaned up with the
342  // arena allocator.
343  mdivide_left_tri_dv_vari<TriView, R1, C1, R2, C2> *baseVari
344  = new mdivide_left_tri_dv_vari<TriView, R1, C1, R2, C2>(A, b);
345 
346  size_t pos = 0;
347  for (size_type j = 0; j < res.cols(); j++)
348  for (size_type i = 0; i < res.rows(); i++)
349  res(i, j).vi_ = baseVari->_variRefC[pos++];
350 
351  return res;
352  }
353  template <int TriView, int R1, int C1, int R2, int C2>
354  inline
355  Eigen::Matrix<var, R1, C2>
356  mdivide_left_tri(const Eigen::Matrix<var, R1, C1> &A,
357  const Eigen::Matrix<double, R2, C2> &b) {
358  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
359 
360  stan::math::check_square("mdivide_left_tri", "A", A);
361  stan::math::check_multiplicable("mdivide_left_tri",
362  "A", A,
363  "b", b);
364 
365  // NOTE: this is not a memory leak, this vari is used in the
366  // expression graph to evaluate the adjoint, but is not needed
367  // for the returned matrix. Memory will be cleaned up with the
368  // arena allocator.
369  mdivide_left_tri_vd_vari<TriView, R1, C1, R2, C2> *baseVari
370  = new mdivide_left_tri_vd_vari<TriView, R1, C1, R2, C2>(A, b);
371 
372  size_t pos = 0;
373  for (size_type j = 0; j < res.cols(); j++)
374  for (size_type i = 0; i < res.rows(); i++)
375  res(i, j).vi_ = baseVari->_variRefC[pos++];
376 
377  return res;
378  }
379 
380  }
381 }
382 #endif
int N_
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double e.
Definition: typedefs.hpp:13
double * C_
Eigen::Matrix< typename boost::math::tools::promote_args< T1, T2 >::type, R1, C2 > mdivide_left_tri(const Eigen::Matrix< T1, R1, C1 > &A, const Eigen::Matrix< T2, R2, C2 > &b)
Returns the solution of the system Ax=b when A is triangular.
size_t rows(const Eigen::Matrix< T, R, C > &m)
Definition: rows.hpp:12
AutodiffStackStorage< chainable, chainable_alloc > ChainableStack
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.
vari ** _variRefC
size_t cols(const Eigen::Matrix< T, R, C > &m)
Definition: cols.hpp:12
double * A_
vari ** _variRefA
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
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.
int M_
vari ** _variRefB

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