mdivide_left.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_HPP
2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_HPP
3 
7 #include <stan/math/rev/core.hpp>
10 #include <vector>
11 
12 namespace stan {
13 namespace math {
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_(reinterpret_cast<double *>(
33  ChainableStack::instance().memalloc_.alloc(sizeof(double) * A.rows()
34  * A.cols()))),
35  C_(reinterpret_cast<double *>(
36  ChainableStack::instance().memalloc_.alloc(sizeof(double) * B.rows()
37  * B.cols()))),
38  variRefA_(reinterpret_cast<vari **>(
39  ChainableStack::instance().memalloc_.alloc(sizeof(vari *) * A.rows()
40  * A.cols()))),
41  variRefB_(reinterpret_cast<vari **>(
42  ChainableStack::instance().memalloc_.alloc(sizeof(vari *) * B.rows()
43  * B.cols()))),
44  variRefC_(reinterpret_cast<vari **>(
45  ChainableStack::instance().memalloc_.alloc(sizeof(vari *) * B.rows()
46  * B.cols()))) {
47  using Eigen::Map;
48  using Eigen::Matrix;
49 
50  size_t pos = 0;
51  for (size_type j = 0; j < M_; j++) {
52  for (size_type i = 0; i < M_; i++) {
53  variRefA_[pos] = A(i, j).vi_;
54  A_[pos++] = A(i, j).val();
55  }
56  }
57 
58  pos = 0;
59  for (size_type j = 0; j < N_; j++) {
60  for (size_type i = 0; i < M_; i++) {
61  variRefB_[pos] = B(i, j).vi_;
62  C_[pos++] = B(i, j).val();
63  }
64  }
65 
66  Matrix<double, R1, C2> C(M_, N_);
67  C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
68 
69  C = Map<Matrix<double, R1, C1> >(A_, M_, M_).colPivHouseholderQr().solve(C);
70 
71  pos = 0;
72  for (size_type j = 0; j < N_; j++) {
73  for (size_type i = 0; i < M_; i++) {
74  C_[pos] = C(i, j);
75  variRefC_[pos] = new vari(C_[pos], false);
76  pos++;
77  }
78  }
79  }
80 
81  virtual void chain() {
82  using Eigen::Map;
83  using Eigen::Matrix;
84  Eigen::Matrix<double, R1, C1> adjA(M_, M_);
85  Eigen::Matrix<double, R2, C2> adjB(M_, N_);
86  Eigen::Matrix<double, R1, C2> adjC(M_, N_);
87 
88  size_t pos = 0;
89  for (size_type j = 0; j < adjC.cols(); j++)
90  for (size_type i = 0; i < adjC.rows(); i++)
91  adjC(i, j) = variRefC_[pos++]->adj_;
92 
93  adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
94  .transpose()
95  .colPivHouseholderQr()
96  .solve(adjC);
97  adjA.noalias()
98  = -adjB * Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose();
99 
100  pos = 0;
101  for (size_type j = 0; j < adjA.cols(); j++)
102  for (size_type i = 0; i < adjA.rows(); i++)
103  variRefA_[pos++]->adj_ += adjA(i, j);
104 
105  pos = 0;
106  for (size_type j = 0; j < adjB.cols(); j++)
107  for (size_type i = 0; i < adjB.rows(); i++)
108  variRefB_[pos++]->adj_ += adjB(i, j);
109  }
110 };
111 
112 template <int R1, int C1, int R2, int C2>
113 class mdivide_left_dv_vari : public vari {
114  public:
115  int M_; // A.rows() = A.cols() = B.rows()
116  int N_; // B.cols()
117  double *A_;
118  double *C_;
119  vari **variRefB_;
120  vari **variRefC_;
121 
122  mdivide_left_dv_vari(const Eigen::Matrix<double, R1, C1> &A,
123  const Eigen::Matrix<var, R2, C2> &B)
124  : vari(0.0),
125  M_(A.rows()),
126  N_(B.cols()),
127  A_(reinterpret_cast<double *>(
128  ChainableStack::instance().memalloc_.alloc(sizeof(double) * A.rows()
129  * A.cols()))),
130  C_(reinterpret_cast<double *>(
131  ChainableStack::instance().memalloc_.alloc(sizeof(double) * B.rows()
132  * B.cols()))),
133  variRefB_(reinterpret_cast<vari **>(
134  ChainableStack::instance().memalloc_.alloc(sizeof(vari *) * B.rows()
135  * B.cols()))),
136  variRefC_(reinterpret_cast<vari **>(
137  ChainableStack::instance().memalloc_.alloc(sizeof(vari *) * B.rows()
138  * B.cols()))) {
139  using Eigen::Map;
140  using Eigen::Matrix;
141 
142  size_t pos = 0;
143  for (size_type j = 0; j < M_; j++) {
144  for (size_type i = 0; i < M_; i++) {
145  A_[pos++] = A(i, j);
146  }
147  }
148 
149  pos = 0;
150  for (size_type j = 0; j < N_; j++) {
151  for (size_type i = 0; i < M_; i++) {
152  variRefB_[pos] = B(i, j).vi_;
153  C_[pos++] = B(i, j).val();
154  }
155  }
156 
157  Matrix<double, R1, C2> C(M_, N_);
158  C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
159 
160  C = Map<Matrix<double, R1, C1> >(A_, M_, M_).colPivHouseholderQr().solve(C);
161 
162  pos = 0;
163  for (size_type j = 0; j < N_; j++) {
164  for (size_type i = 0; i < M_; i++) {
165  C_[pos] = C(i, j);
166  variRefC_[pos] = new vari(C_[pos], false);
167  pos++;
168  }
169  }
170  }
171 
172  virtual void chain() {
173  using Eigen::Map;
174  using Eigen::Matrix;
175  Eigen::Matrix<double, R2, C2> adjB(M_, N_);
176  Eigen::Matrix<double, R1, C2> adjC(M_, N_);
177 
178  size_t pos = 0;
179  for (size_type j = 0; j < adjC.cols(); j++)
180  for (size_type i = 0; i < adjC.rows(); i++)
181  adjC(i, j) = variRefC_[pos++]->adj_;
182 
183  adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
184  .transpose()
185  .colPivHouseholderQr()
186  .solve(adjC);
187 
188  pos = 0;
189  for (size_type j = 0; j < adjB.cols(); j++)
190  for (size_type i = 0; i < adjB.rows(); i++)
191  variRefB_[pos++]->adj_ += adjB(i, j);
192  }
193 };
194 
195 template <int R1, int C1, int R2, int C2>
196 class mdivide_left_vd_vari : public vari {
197  public:
198  int M_; // A.rows() = A.cols() = B.rows()
199  int N_; // B.cols()
200  double *A_;
201  double *C_;
202  vari **variRefA_;
203  vari **variRefC_;
204 
205  mdivide_left_vd_vari(const Eigen::Matrix<var, R1, C1> &A,
206  const Eigen::Matrix<double, R2, C2> &B)
207  : vari(0.0),
208  M_(A.rows()),
209  N_(B.cols()),
210  A_(reinterpret_cast<double *>(
211  ChainableStack::instance().memalloc_.alloc(sizeof(double) * A.rows()
212  * A.cols()))),
213  C_(reinterpret_cast<double *>(
214  ChainableStack::instance().memalloc_.alloc(sizeof(double) * B.rows()
215  * B.cols()))),
216  variRefA_(reinterpret_cast<vari **>(
217  ChainableStack::instance().memalloc_.alloc(sizeof(vari *) * A.rows()
218  * A.cols()))),
219  variRefC_(reinterpret_cast<vari **>(
220  ChainableStack::instance().memalloc_.alloc(sizeof(vari *) * B.rows()
221  * B.cols()))) {
222  using Eigen::Map;
223  using Eigen::Matrix;
224 
225  size_t pos = 0;
226  for (size_type j = 0; j < M_; j++) {
227  for (size_type i = 0; i < M_; i++) {
228  variRefA_[pos] = A(i, j).vi_;
229  A_[pos++] = A(i, j).val();
230  }
231  }
232 
233  Matrix<double, R1, C2> C(M_, N_);
234  C = Map<Matrix<double, R1, C1> >(A_, M_, M_).colPivHouseholderQr().solve(B);
235 
236  pos = 0;
237  for (size_type j = 0; j < N_; j++) {
238  for (size_type i = 0; i < M_; i++) {
239  C_[pos] = C(i, j);
240  variRefC_[pos] = new vari(C_[pos], false);
241  pos++;
242  }
243  }
244  }
245 
246  virtual void chain() {
247  using Eigen::Map;
248  using Eigen::Matrix;
249  Eigen::Matrix<double, R1, C1> adjA(M_, M_);
250  Eigen::Matrix<double, R1, C2> adjC(M_, N_);
251 
252  size_t pos = 0;
253  for (size_type j = 0; j < adjC.cols(); j++)
254  for (size_type i = 0; i < adjC.rows(); i++)
255  adjC(i, j) = variRefC_[pos++]->adj_;
256 
257  // FIXME: add .noalias() to LHS
258  adjA = -Map<Matrix<double, R1, C1> >(A_, M_, M_)
259  .transpose()
260  .colPivHouseholderQr()
261  .solve(adjC
262  * Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose());
263 
264  pos = 0;
265  for (size_type j = 0; j < adjA.cols(); j++)
266  for (size_type i = 0; i < adjA.rows(); i++)
267  variRefA_[pos++]->adj_ += adjA(i, j);
268  }
269 };
270 } // namespace
271 
272 template <int R1, int C1, int R2, int C2>
273 inline Eigen::Matrix<var, R1, C2> mdivide_left(
274  const Eigen::Matrix<var, R1, C1> &A, const Eigen::Matrix<var, R2, C2> &b) {
275  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
276 
277  check_square("mdivide_left", "A", A);
278  check_multiplicable("mdivide_left", "A", A, "b", b);
279 
280  // NOTE: this is not a memory leak, this vari is used in the
281  // expression graph to evaluate the adjoint, but is not needed
282  // for the returned matrix. Memory will be cleaned up with the
283  // arena allocator.
284  mdivide_left_vv_vari<R1, C1, R2, C2> *baseVari
285  = new mdivide_left_vv_vari<R1, C1, R2, C2>(A, b);
286 
287  size_t pos = 0;
288  for (size_type j = 0; j < res.cols(); j++)
289  for (size_type i = 0; i < res.rows(); i++)
290  res(i, j).vi_ = baseVari->variRefC_[pos++];
291 
292  return res;
293 }
294 
295 template <int R1, int C1, int R2, int C2>
296 inline Eigen::Matrix<var, R1, C2> mdivide_left(
297  const Eigen::Matrix<var, R1, C1> &A,
298  const Eigen::Matrix<double, R2, C2> &b) {
299  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
300 
301  check_square("mdivide_left", "A", A);
302  check_multiplicable("mdivide_left", "A", A, "b", b);
303 
304  // NOTE: this is not a memory leak, this vari is used in the
305  // expression graph to evaluate the adjoint, but is not needed
306  // for the returned matrix. Memory will be cleaned up with the
307  // arena allocator.
308  mdivide_left_vd_vari<R1, C1, R2, C2> *baseVari
309  = new mdivide_left_vd_vari<R1, C1, R2, C2>(A, b);
310 
311  size_t pos = 0;
312  for (size_type j = 0; j < res.cols(); j++)
313  for (size_type i = 0; i < res.rows(); i++)
314  res(i, j).vi_ = baseVari->variRefC_[pos++];
315 
316  return res;
317 }
318 
319 template <int R1, int C1, int R2, int C2>
320 inline Eigen::Matrix<var, R1, C2> mdivide_left(
321  const Eigen::Matrix<double, R1, C1> &A,
322  const Eigen::Matrix<var, R2, C2> &b) {
323  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
324 
325  check_square("mdivide_left", "A", A);
326  check_multiplicable("mdivide_left", "A", A, "b", b);
327 
328  // NOTE: this is not a memory leak, this vari is used in the
329  // expression graph to evaluate the adjoint, but is not needed
330  // for the returned matrix. Memory will be cleaned up with the
331  // arena allocator.
332  mdivide_left_dv_vari<R1, C1, R2, C2> *baseVari
333  = new mdivide_left_dv_vari<R1, C1, R2, C2>(A, b);
334 
335  size_t pos = 0;
336  for (size_type j = 0; j < res.cols(); j++)
337  for (size_type i = 0; i < res.rows(); i++)
338  res(i, j).vi_ = baseVari->variRefC_[pos++];
339 
340  return res;
341 }
342 
343 } // namespace math
344 } // namespace stan
345 #endif
Eigen::Matrix< fvar< T >, R1, C2 > mdivide_left(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
int rows(const Eigen::Matrix< T, R, C > &m)
Definition: rows.hpp:20
double * C_
vari ** variRefB_
const double C
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
int M_
chain
Check that an output directory exists.
const double j
Definition: BetheBloch.cxx:29
AutodiffStackSingleton< vari, chainable_alloc > ChainableStack
int cols(const Eigen::Matrix< T, R, C > &m)
Definition: cols.hpp:20
double * A_
void check_multiplicable(const char *function, const char *name1, const T1 &y1, const char *name2, const T2 &y2)
void check_square(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
static const double A
Definition: Units.h:82
static AutodiffStackStorage & instance()
const hit & b
Definition: hits.cxx:21
vari ** variRefA_
vari ** variRefC_
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:10
int N_
void * alloc(size_t len)