Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
mdivide_left_spd.hpp
Go to the documentation of this file.
1 #ifndef __STAN__AGRAD__REV__MATRIX__MDIVIDE_LEFT_SPD_HPP__
2 #define __STAN__AGRAD__REV__MATRIX__MDIVIDE_LEFT_SPD_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 R1,int C1,int R2,int C2>
17  class mdivide_left_spd_alloc : public chainable_alloc {
18  public:
19  virtual ~mdivide_left_spd_alloc() {}
20 
21  Eigen::LLT< Eigen::Matrix<double,R1,C1> > _llt;
22  Eigen::Matrix<double,R2,C2> _C;
23  };
24 
25  template <int R1,int C1,int R2,int C2>
26  class mdivide_left_spd_vv_vari : public vari {
27  public:
28  int _M; // A.rows() = A.cols() = B.rows()
29  int _N; // B.cols()
30  vari** _variRefA;
31  vari** _variRefB;
32  vari** _variRefC;
33  mdivide_left_spd_alloc<R1,C1,R2,C2> *_alloc;
34 
35  mdivide_left_spd_vv_vari(const Eigen::Matrix<var,R1,C1> &A,
36  const Eigen::Matrix<var,R2,C2> &B)
37  : vari(0.0),
38  _M(A.rows()),
39  _N(B.cols()),
40  _variRefA((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*)
41  * A.rows() * A.cols())),
42  _variRefB((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*)
43  * B.rows() * B.cols())),
44  _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*)
45  * B.rows() * B.cols())),
46  _alloc(new mdivide_left_spd_alloc<R1,C1,R2,C2>())
47  {
48  using Eigen::Matrix;
49  using Eigen::Map;
50 
51  Matrix<double,R1,C1> Ad(A.rows(),A.cols());
52 
53  size_t pos = 0;
54  for (size_type j = 0; j < _M; j++) {
55  for (size_type i = 0; i < _M; i++) {
56  _variRefA[pos] = A(i,j).vi_;
57  Ad(i,j) = A(i,j).val();
58  pos++;
59  }
60  }
61 
62  pos = 0;
63  _alloc->_C.resize(_M,_N);
64  for (size_type j = 0; j < _N; j++) {
65  for (size_type i = 0; i < _M; i++) {
66  _variRefB[pos] = B(i,j).vi_;
67  _alloc->_C(i,j) = B(i,j).val();
68  pos++;
69  }
70  }
71 
72  _alloc->_llt = Ad.llt();
73  _alloc->_llt.solveInPlace(_alloc->_C);
74 
75  pos = 0;
76  for (size_type j = 0; j < _N; j++) {
77  for (size_type i = 0; i < _M; i++) {
78  _variRefC[pos] = new vari(_alloc->_C(i,j),false);
79  pos++;
80  }
81  }
82  }
83 
84  virtual void chain() {
85  using Eigen::Matrix;
86  using Eigen::Map;
87  Eigen::Matrix<double,R1,C1> adjA(_M,_M);
88  Eigen::Matrix<double,R2,C2> adjB(_M,_N);
89 
90  size_t pos = 0;
91  for (size_type j = 0; j < _N; j++)
92  for (size_type i = 0; i < _M; i++)
93  adjB(i,j) = _variRefC[pos++]->adj_;
94 
95  _alloc->_llt.solveInPlace(adjB);
96  adjA.noalias() = -adjB * _alloc->_C.transpose();
97 
98  pos = 0;
99  for (size_type j = 0; j < _M; j++)
100  for (size_type i = 0; i < _M; i++)
101  _variRefA[pos++]->adj_ += adjA(i,j);
102 
103  pos = 0;
104  for (size_type j = 0; j < _N; j++)
105  for (size_type i = 0; i < _M; i++)
106  _variRefB[pos++]->adj_ += adjB(i,j);
107  }
108  };
109 
110  template <int R1,int C1,int R2,int C2>
111  class mdivide_left_spd_dv_vari : public vari {
112  public:
113  int _M; // A.rows() = A.cols() = B.rows()
114  int _N; // B.cols()
115  vari** _variRefB;
116  vari** _variRefC;
117  mdivide_left_spd_alloc<R1,C1,R2,C2> *_alloc;
118 
119  mdivide_left_spd_dv_vari(const Eigen::Matrix<double,R1,C1> &A,
120  const Eigen::Matrix<var,R2,C2> &B)
121  : vari(0.0),
122  _M(A.rows()),
123  _N(B.cols()),
124  _variRefB((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*)
125  * B.rows() * B.cols())),
126  _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*)
127  * B.rows() * B.cols())),
128  _alloc(new mdivide_left_spd_alloc<R1,C1,R2,C2>())
129  {
130  using Eigen::Matrix;
131  using Eigen::Map;
132 
133  size_t pos = 0;
134  _alloc->_C.resize(_M,_N);
135  for (size_type j = 0; j < _N; j++) {
136  for (size_type i = 0; i < _M; i++) {
137  _variRefB[pos] = B(i,j).vi_;
138  _alloc->_C(i,j) = B(i,j).val();
139  pos++;
140  }
141  }
142 
143  _alloc->_llt = A.llt();
144  _alloc->_llt.solveInPlace(_alloc->_C);
145 
146  pos = 0;
147  for (size_type j = 0; j < _N; j++) {
148  for (size_type i = 0; i < _M; i++) {
149  _variRefC[pos] = new vari(_alloc->_C(i,j),false);
150  pos++;
151  }
152  }
153  }
154 
155  virtual void chain() {
156  using Eigen::Matrix;
157  using Eigen::Map;
158  Eigen::Matrix<double,R2,C2> adjB(_M,_N);
159 
160  size_t pos = 0;
161  for (size_type j = 0; j < adjB.cols(); j++)
162  for (size_type i = 0; i < adjB.rows(); i++)
163  adjB(i,j) = _variRefC[pos++]->adj_;
164 
165  _alloc->_llt.solveInPlace(adjB);
166 
167  pos = 0;
168  for (size_type j = 0; j < adjB.cols(); j++)
169  for (size_type i = 0; i < adjB.rows(); i++)
170  _variRefB[pos++]->adj_ += adjB(i,j);
171  }
172  };
173 
174  template <int R1,int C1,int R2,int C2>
175  class mdivide_left_spd_vd_vari : public vari {
176  public:
177  int _M; // A.rows() = A.cols() = B.rows()
178  int _N; // B.cols()
179  vari** _variRefA;
180  vari** _variRefC;
181  mdivide_left_spd_alloc<R1,C1,R2,C2> *_alloc;
182 
183  mdivide_left_spd_vd_vari(const Eigen::Matrix<var,R1,C1> &A,
184  const Eigen::Matrix<double,R2,C2> &B)
185  : vari(0.0),
186  _M(A.rows()),
187  _N(B.cols()),
188  _variRefA((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*)
189  * A.rows() * A.cols())),
190  _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*)
191  * B.rows() * B.cols())),
192  _alloc(new mdivide_left_spd_alloc<R1,C1,R2,C2>())
193  {
194  using Eigen::Matrix;
195  using Eigen::Map;
196 
197  Matrix<double,R1,C1> Ad(A.rows(),A.cols());
198 
199  size_t pos = 0;
200  for (size_type j = 0; j < _M; j++) {
201  for (size_type i = 0; i < _M; i++) {
202  _variRefA[pos] = A(i,j).vi_;
203  Ad(i,j) = A(i,j).val();
204  pos++;
205  }
206  }
207 
208  _alloc->_llt = Ad.llt();
209  _alloc->_C = _alloc->_llt.solve(B);
210 
211  pos = 0;
212  for (size_type j = 0; j < _N; j++) {
213  for (size_type i = 0; i < _M; i++) {
214  _variRefC[pos] = new vari(_alloc->_C(i,j),false);
215  pos++;
216  }
217  }
218  }
219 
220  virtual void chain() {
221  using Eigen::Matrix;
222  using Eigen::Map;
223  Eigen::Matrix<double,R1,C1> adjA(_M,_M);
224  Eigen::Matrix<double,R1,C2> adjC(_M,_N);
225 
226  size_t pos = 0;
227  for (size_type j = 0; j < adjC.cols(); j++)
228  for (size_type i = 0; i < adjC.rows(); i++)
229  adjC(i,j) = _variRefC[pos++]->adj_;
230 
231  adjA = -_alloc->_llt.solve(adjC*_alloc->_C.transpose());
232 
233  pos = 0;
234  for (size_type j = 0; j < adjA.cols(); j++)
235  for (size_type i = 0; i < adjA.rows(); i++)
236  _variRefA[pos++]->adj_ += adjA(i,j);
237  }
238  };
239  }
240 
241  template <int R1,int C1,int R2,int C2>
242  inline
243  Eigen::Matrix<var,R1,C2>
244  mdivide_left_spd(const Eigen::Matrix<var,R1,C1> &A,
245  const Eigen::Matrix<var,R2,C2> &b) {
246  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
247 
248  stan::math::validate_square(A,"mdivide_left_spd");
249  stan::math::validate_multiplicable(A,b,"mdivide_left_spd");
250 
251  // NOTE: this is not a memory leak, this vari is used in the
252  // expression graph to evaluate the adjoint, but is not needed
253  // for the returned matrix. Memory will be cleaned up with the arena allocator.
254  mdivide_left_spd_vv_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_spd_vv_vari<R1,C1,R2,C2>(A,b);
255 
256  size_t pos = 0;
257  for (size_type j = 0; j < res.cols(); j++)
258  for (size_type i = 0; i < res.rows(); i++)
259  res(i,j).vi_ = baseVari->_variRefC[pos++];
260 
261  return res;
262  }
263 
264  template <int R1,int C1,int R2,int C2>
265  inline
266  Eigen::Matrix<var,R1,C2>
267  mdivide_left_spd(const Eigen::Matrix<var,R1,C1> &A,
268  const Eigen::Matrix<double,R2,C2> &b) {
269  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
270 
271  stan::math::validate_square(A,"mdivide_left_spd");
272  stan::math::validate_multiplicable(A,b,"mdivide_left_spd");
273 
274  // NOTE: this is not a memory leak, this vari is used in the
275  // expression graph to evaluate the adjoint, but is not needed
276  // for the returned matrix. Memory will be cleaned up with the arena allocator.
277  mdivide_left_spd_vd_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_spd_vd_vari<R1,C1,R2,C2>(A,b);
278 
279  size_t pos = 0;
280  for (size_type j = 0; j < res.cols(); j++)
281  for (size_type i = 0; i < res.rows(); i++)
282  res(i,j).vi_ = baseVari->_variRefC[pos++];
283 
284  return res;
285  }
286 
287  template <int R1,int C1,int R2,int C2>
288  inline
289  Eigen::Matrix<var,R1,C2>
290  mdivide_left_spd(const Eigen::Matrix<double,R1,C1> &A,
291  const Eigen::Matrix<var,R2,C2> &b) {
292  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
293 
294  stan::math::validate_square(A,"mdivide_left_spd");
295  stan::math::validate_multiplicable(A,b,"mdivide_left_spd");
296 
297  // NOTE: this is not a memory leak, this vari is used in the
298  // expression graph to evaluate the adjoint, but is not needed
299  // for the returned matrix. Memory will be cleaned up with the arena allocator.
300  mdivide_left_spd_dv_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_spd_dv_vari<R1,C1,R2,C2>(A,b);
301 
302  size_t pos = 0;
303  for (size_type j = 0; j < res.cols(); j++)
304  for (size_type i = 0; i < res.rows(); i++)
305  res(i,j).vi_ = baseVari->_variRefC[pos++];
306 
307  return res;
308  }
309 
310  }
311 }
312 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.