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