Stan Math Library  2.20.0
reverse mode automatic differentiation
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 
4 #include <stan/math/rev/meta.hpp>
8 #include <stan/math/rev/core.hpp>
10 
11 namespace stan {
12 namespace math {
13 
14 namespace internal {
15 template <int TriView, int R1, int C1, int R2, int C2>
17  public:
18  int M_; // A.rows() = A.cols() = B.rows()
19  int N_; // B.cols()
20  double *A_;
21  double *C_;
25 
26  mdivide_left_tri_vv_vari(const Eigen::Matrix<var, R1, C1> &A,
27  const Eigen::Matrix<var, R2, C2> &B)
28  : vari(0.0),
29  M_(A.rows()),
30  N_(B.cols()),
31  A_(reinterpret_cast<double *>(
32  ChainableStack::instance_->memalloc_.alloc(sizeof(double) * A.rows()
33  * A.cols()))),
34  C_(reinterpret_cast<double *>(
35  ChainableStack::instance_->memalloc_.alloc(sizeof(double) * B.rows()
36  * B.cols()))),
37  variRefA_(reinterpret_cast<vari **>(
38  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * A.rows()
39  * (A.rows() + 1) / 2))),
40  variRefB_(reinterpret_cast<vari **>(
41  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
42  * B.cols()))),
43  variRefC_(reinterpret_cast<vari **>(
44  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
45  * B.cols()))) {
46  using Eigen::Map;
47  using Eigen::Matrix;
48 
49  size_t pos = 0;
50  if (TriView == Eigen::Lower) {
51  for (size_type j = 0; j < M_; j++)
52  for (size_type i = j; i < M_; i++)
53  variRefA_[pos++] = A(i, j).vi_;
54  } else if (TriView == Eigen::Upper) {
55  for (size_type j = 0; j < M_; j++)
56  for (size_type i = 0; i < j + 1; i++)
57  variRefA_[pos++] = A(i, j).vi_;
58  }
59 
60  pos = 0;
61  for (size_type j = 0; j < M_; j++) {
62  for (size_type i = 0; i < M_; i++) {
63  A_[pos++] = A(i, j).val();
64  }
65  }
66 
67  pos = 0;
68  for (size_type j = 0; j < N_; j++) {
69  for (size_type i = 0; i < M_; i++) {
70  variRefB_[pos] = B(i, j).vi_;
71  C_[pos++] = B(i, j).val();
72  }
73  }
74 
75  Matrix<double, R1, C2> C(M_, N_);
76  C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
77 
78  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
79  .template triangularView<TriView>()
80  .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::Map;
94  using Eigen::Matrix;
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>()
106  .transpose()
107  .solve(adjC);
108  adjA.noalias()
109  = -adjB * Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose();
110 
111  pos = 0;
112  if (TriView == Eigen::Lower) {
113  for (size_type j = 0; j < adjA.cols(); j++)
114  for (size_type i = j; i < adjA.rows(); i++)
115  variRefA_[pos++]->adj_ += adjA(i, j);
116  } else if (TriView == Eigen::Upper) {
117  for (size_type j = 0; j < adjA.cols(); j++)
118  for (size_type i = 0; i < j + 1; i++)
119  variRefA_[pos++]->adj_ += adjA(i, j);
120  }
121 
122  pos = 0;
123  for (size_type j = 0; j < adjB.cols(); j++)
124  for (size_type i = 0; i < adjB.rows(); i++)
125  variRefB_[pos++]->adj_ += adjB(i, j);
126  }
127 };
128 
129 template <int TriView, int R1, int C1, int R2, int C2>
131  public:
132  int M_; // A.rows() = A.cols() = B.rows()
133  int N_; // B.cols()
134  double *A_;
135  double *C_;
138 
139  mdivide_left_tri_dv_vari(const Eigen::Matrix<double, R1, C1> &A,
140  const Eigen::Matrix<var, R2, C2> &B)
141  : vari(0.0),
142  M_(A.rows()),
143  N_(B.cols()),
144  A_(reinterpret_cast<double *>(
145  ChainableStack::instance_->memalloc_.alloc(sizeof(double) * A.rows()
146  * A.cols()))),
147  C_(reinterpret_cast<double *>(
148  ChainableStack::instance_->memalloc_.alloc(sizeof(double) * B.rows()
149  * B.cols()))),
150  variRefB_(reinterpret_cast<vari **>(
151  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
152  * B.cols()))),
153  variRefC_(reinterpret_cast<vari **>(
154  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
155  * B.cols()))) {
156  using Eigen::Map;
157  using Eigen::Matrix;
158 
159  size_t pos = 0;
160  for (size_type j = 0; j < M_; j++) {
161  for (size_type i = 0; i < M_; i++) {
162  A_[pos++] = A(i, j);
163  }
164  }
165 
166  pos = 0;
167  for (size_type j = 0; j < N_; j++) {
168  for (size_type i = 0; i < M_; i++) {
169  variRefB_[pos] = B(i, j).vi_;
170  C_[pos++] = B(i, j).val();
171  }
172  }
173 
174  Matrix<double, R1, C2> C(M_, N_);
175  C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
176 
177  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
178  .template triangularView<TriView>()
179  .solve(C);
180 
181  pos = 0;
182  for (size_type j = 0; j < N_; j++) {
183  for (size_type i = 0; i < M_; i++) {
184  C_[pos] = C(i, j);
185  variRefC_[pos] = new vari(C_[pos], false);
186  pos++;
187  }
188  }
189  }
190 
191  virtual void chain() {
192  using Eigen::Map;
193  using Eigen::Matrix;
194  Matrix<double, R2, C2> adjB(M_, N_);
195  Matrix<double, R1, C2> adjC(M_, N_);
196 
197  size_t pos = 0;
198  for (size_type j = 0; j < adjC.cols(); j++)
199  for (size_type i = 0; i < adjC.rows(); i++)
200  adjC(i, j) = variRefC_[pos++]->adj_;
201 
202  adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
203  .template triangularView<TriView>()
204  .transpose()
205  .solve(adjC);
206 
207  pos = 0;
208  for (size_type j = 0; j < adjB.cols(); j++)
209  for (size_type i = 0; i < adjB.rows(); i++)
210  variRefB_[pos++]->adj_ += adjB(i, j);
211  }
212 };
213 
214 template <int TriView, int R1, int C1, int R2, int C2>
216  public:
217  int M_; // A.rows() = A.cols() = B.rows()
218  int N_; // B.cols()
219  double *A_;
220  double *C_;
223 
224  mdivide_left_tri_vd_vari(const Eigen::Matrix<var, R1, C1> &A,
225  const Eigen::Matrix<double, R2, C2> &B)
226  : vari(0.0),
227  M_(A.rows()),
228  N_(B.cols()),
229  A_(reinterpret_cast<double *>(
230  ChainableStack::instance_->memalloc_.alloc(sizeof(double) * A.rows()
231  * A.cols()))),
232  C_(reinterpret_cast<double *>(
233  ChainableStack::instance_->memalloc_.alloc(sizeof(double) * B.rows()
234  * B.cols()))),
235  variRefA_(reinterpret_cast<vari **>(
236  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * A.rows()
237  * (A.rows() + 1) / 2))),
238  variRefC_(reinterpret_cast<vari **>(
239  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
240  * B.cols()))) {
241  using Eigen::Map;
242  using Eigen::Matrix;
243 
244  size_t pos = 0;
245  if (TriView == Eigen::Lower) {
246  for (size_type j = 0; j < M_; j++)
247  for (size_type i = j; i < M_; i++)
248  variRefA_[pos++] = A(i, j).vi_;
249  } else if (TriView == Eigen::Upper) {
250  for (size_type j = 0; j < M_; j++)
251  for (size_type i = 0; i < j + 1; i++)
252  variRefA_[pos++] = A(i, j).vi_;
253  }
254 
255  pos = 0;
256  for (size_type j = 0; j < M_; j++) {
257  for (size_type i = 0; i < M_; i++) {
258  A_[pos++] = A(i, j).val();
259  }
260  }
261 
262  Matrix<double, R1, C2> C(M_, N_);
263  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
264  .template triangularView<TriView>()
265  .solve(B);
266 
267  pos = 0;
268  for (size_type j = 0; j < N_; j++) {
269  for (size_type i = 0; i < M_; i++) {
270  C_[pos] = C(i, j);
271  variRefC_[pos] = new vari(C_[pos], false);
272  pos++;
273  }
274  }
275  }
276 
277  virtual void chain() {
278  using Eigen::Map;
279  using Eigen::Matrix;
280  Matrix<double, R1, C1> adjA(M_, M_);
281  Matrix<double, R1, C2> adjC(M_, N_);
282 
283  size_t pos = 0;
284  for (size_type j = 0; j < adjC.cols(); j++)
285  for (size_type i = 0; i < adjC.rows(); i++)
286  adjC(i, j) = variRefC_[pos++]->adj_;
287 
288  adjA.noalias()
289  = -Map<Matrix<double, R1, C1> >(A_, M_, M_)
290  .template triangularView<TriView>()
291  .transpose()
292  .solve(adjC
293  * Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose());
294 
295  pos = 0;
296  if (TriView == Eigen::Lower) {
297  for (size_type j = 0; j < adjA.cols(); j++)
298  for (size_type i = j; i < adjA.rows(); i++)
299  variRefA_[pos++]->adj_ += adjA(i, j);
300  } else if (TriView == Eigen::Upper) {
301  for (size_type j = 0; j < adjA.cols(); j++)
302  for (size_type i = 0; i < j + 1; i++)
303  variRefA_[pos++]->adj_ += adjA(i, j);
304  }
305  }
306 };
307 } // namespace internal
308 
309 template <int TriView, int R1, int C1, int R2, int C2>
310 inline Eigen::Matrix<var, R1, C2> mdivide_left_tri(
311  const Eigen::Matrix<var, R1, C1> &A, const Eigen::Matrix<var, R2, C2> &b) {
312  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
313 
314  check_square("mdivide_left_tri", "A", A);
315  check_multiplicable("mdivide_left_tri", "A", A, "b", b);
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
320  // arena allocator.
323 
324  size_t pos = 0;
325  for (size_type j = 0; j < res.cols(); j++)
326  for (size_type i = 0; i < res.rows(); i++)
327  res(i, j).vi_ = baseVari->variRefC_[pos++];
328 
329  return res;
330 }
331 template <int TriView, int R1, int C1, int R2, int C2>
332 inline Eigen::Matrix<var, R1, C2> mdivide_left_tri(
333  const Eigen::Matrix<double, R1, C1> &A,
334  const Eigen::Matrix<var, R2, C2> &b) {
335  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
336 
337  check_square("mdivide_left_tri", "A", A);
338  check_multiplicable("mdivide_left_tri", "A", A, "b", b);
339 
340  // NOTE: this is not a memory leak, this vari is used in the
341  // expression graph to evaluate the adjoint, but is not needed
342  // for the returned matrix. Memory will be cleaned up with the
343  // arena allocator.
346 
347  size_t pos = 0;
348  for (size_type j = 0; j < res.cols(); j++)
349  for (size_type i = 0; i < res.rows(); i++)
350  res(i, j).vi_ = baseVari->variRefC_[pos++];
351 
352  return res;
353 }
354 template <int TriView, int R1, int C1, int R2, int C2>
355 inline Eigen::Matrix<var, R1, C2> mdivide_left_tri(
356  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  check_square("mdivide_left_tri", "A", A);
361  check_multiplicable("mdivide_left_tri", "A", A, "b", b);
362 
363  // NOTE: this is not a memory leak, this vari is used in the
364  // expression graph to evaluate the adjoint, but is not needed
365  // for the returned matrix. Memory will be cleaned up with the
366  // arena allocator.
369 
370  size_t pos = 0;
371  for (size_type j = 0; j < res.cols(); j++)
372  for (size_type i = 0; i < res.rows(); i++)
373  res(i, j).vi_ = baseVari->variRefC_[pos++];
374 
375  return res;
376 }
377 
378 } // namespace math
379 } // namespace stan
380 #endif
int rows(const Eigen::Matrix< T, R, C > &m)
Return the number of rows in the specified matrix, vector, or row vector.
Definition: rows.hpp:20
mdivide_left_tri_vd_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< double, R2, C2 > &B)
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
void check_square(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is square.
The variable implementation base class.
Definition: vari.hpp:30
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:11
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.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
int cols(const Eigen::Matrix< T, R, C > &m)
Return the number of columns in the specified matrix, vector, or row vector.
Definition: cols.hpp:20
void check_multiplicable(const char *function, const char *name1, const T1 &y1, const char *name2, const T2 &y2)
Check if the matrices can be multiplied.
matrix_cl transpose(const matrix_cl &src)
Takes the transpose of the matrix on the OpenCL device.
Definition: transpose.hpp:20
mdivide_left_tri_vv_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)
vari(double x)
Construct a variable implementation from a value.
Definition: vari.hpp:58
mdivide_left_tri_dv_vari(const Eigen::Matrix< double, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
This struct always provides access to the autodiff stack using the singleton pattern.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...

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