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