| Stan Math Library 2.12.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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

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