1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_HPP
2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_HPP
16 template <
int R1,
int C1,
int R2,
int C2>
17 class mdivide_left_vv_vari :
public vari {
27 mdivide_left_vv_vari(
const Eigen::Matrix<var, R1, C1> &A,
28 const Eigen::Matrix<var, R2, C2> &B)
32 A_(reinterpret_cast<double*>
34 .alloc(sizeof(double) * A.
rows() * A.
cols()))),
35 C_(reinterpret_cast<double*>
37 .alloc(sizeof(double) * B.
rows() * B.
cols()))),
38 _variRefA(reinterpret_cast<vari**>
40 .alloc(sizeof(vari*) * A.
rows() * A.
cols()))),
41 _variRefB(reinterpret_cast<vari**>
43 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
44 _variRefC(reinterpret_cast<vari**>
46 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))) {
53 _variRefA[pos] = A(i, j).vi_;
54 A_[pos++] = A(i, j).val();
61 _variRefB[pos] = B(i, j).vi_;
62 C_[pos++] = B(i, j).val();
66 Matrix<double, R1, C2> C(M_, N_);
67 C = Map<Matrix<double, R1, C2> >(
C_,
M_,
N_);
69 C = Map<Matrix<double, R1, C1> >(
A_,
M_,
M_)
70 .colPivHouseholderQr().solve(C);
76 _variRefC[pos] =
new vari(C_[pos],
false);
82 virtual void chain() {
85 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
86 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
87 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
90 for (
size_type j = 0; j < adjC.cols(); j++)
91 for (
size_type i = 0; i < adjC.rows(); i++)
92 adjC(i, j) = _variRefC[pos++]->adj_;
95 adjB = Map<Matrix<double, R1, C1> >(
A_,
M_,
M_)
96 .
transpose().colPivHouseholderQr().solve(adjC);
97 adjA.noalias() = -adjB
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);
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);
112 template <
int R1,
int C1,
int R2,
int C2>
113 class mdivide_left_dv_vari :
public vari {
122 mdivide_left_dv_vari(
const Eigen::Matrix<double, R1, C1> &A,
123 const Eigen::Matrix<var, R2, C2> &B)
127 A_(reinterpret_cast<double*>
129 .alloc(sizeof(double) * A.
rows() * A.
cols()))),
130 C_(reinterpret_cast<double*>
132 .alloc(sizeof(double) * B.
rows() * B.
cols()))),
133 _variRefB(reinterpret_cast<vari**>
135 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))),
136 _variRefC(reinterpret_cast<vari**>
138 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))) {
152 _variRefB[pos] = B(i, j).vi_;
153 C_[pos++] = B(i, j).val();
157 Matrix<double, R1, C2> C(M_, N_);
158 C = Map<Matrix<double, R1, C2> >(
C_,
M_,
N_);
160 C = Map<Matrix<double, R1, C1> >(
A_,
M_,
M_)
161 .colPivHouseholderQr().solve(C);
167 _variRefC[pos] =
new vari(C_[pos],
false);
173 virtual void chain() {
176 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
177 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
180 for (
size_type j = 0; j < adjC.cols(); j++)
181 for (
size_type i = 0; i < adjC.rows(); i++)
182 adjC(i, j) = _variRefC[pos++]->adj_;
184 adjB = Map<Matrix<double, R1, C1> >(
A_,
M_,
M_)
185 .
transpose().colPivHouseholderQr().solve(adjC);
188 for (
size_type j = 0; j < adjB.cols(); j++)
189 for (
size_type i = 0; i < adjB.rows(); i++)
190 _variRefB[pos++]->adj_ += adjB(i, j);
194 template <
int R1,
int C1,
int R2,
int C2>
195 class mdivide_left_vd_vari :
public vari {
204 mdivide_left_vd_vari(
const Eigen::Matrix<var, R1, C1> &A,
205 const Eigen::Matrix<double, R2, C2> &B)
209 A_(reinterpret_cast<double*>
211 .alloc(sizeof(double) * A.
rows() * A.
cols()))),
212 C_(reinterpret_cast<double*>
214 .alloc(sizeof(double) * B.
rows() * B.
cols()))),
215 _variRefA(reinterpret_cast<vari**>
217 .alloc(sizeof(vari*) * A.
rows() * A.
cols()))),
218 _variRefC(reinterpret_cast<vari**>
220 .alloc(sizeof(vari*) * B.
rows() * B.
cols()))) {
227 _variRefA[pos] = A(i, j).vi_;
228 A_[pos++] = A(i, j).val();
232 Matrix<double, R1, C2> C(M_, N_);
233 C = Map<Matrix<double, R1, C1> >(
A_,
M_,
M_)
234 .colPivHouseholderQr().solve(B);
240 _variRefC[pos] =
new vari(C_[pos],
false);
246 virtual void chain() {
249 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
250 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
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_;
258 adjA = -Map<Matrix<double, R1, C1> >(
A_,
M_,
M_)
260 .colPivHouseholderQr()
261 .solve(adjC*Map<Matrix<double, R1, C2> >(C_, M_, N_).
transpose());
264 for (
size_type j = 0; j < adjA.cols(); j++)
265 for (
size_type i = 0; i < adjA.rows(); i++)
266 _variRefA[pos++]->adj_ += adjA(i, j);
271 template <
int R1,
int C1,
int R2,
int C2>
273 Eigen::Matrix<var, R1, C2>
275 const Eigen::Matrix<var, R2, C2> &b) {
276 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
287 mdivide_left_vv_vari<R1, C1, R2, C2> *baseVari
288 =
new mdivide_left_vv_vari<R1, C1, R2, C2>(A, b);
291 for (
size_type j = 0; j < res.cols(); j++)
292 for (
size_type i = 0; i < res.rows(); i++)
293 res(i, j).vi_ = baseVari->_variRefC[pos++];
298 template <
int R1,
int C1,
int R2,
int C2>
300 Eigen::Matrix<var, R1, C2>
302 const Eigen::Matrix<double, R2, C2> &b) {
303 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
314 mdivide_left_vd_vari<R1, C1, R2, C2> *baseVari
315 =
new mdivide_left_vd_vari<R1, C1, R2, C2>(A, b);
318 for (
size_type j = 0; j < res.cols(); j++)
319 for (
size_type i = 0; i < res.rows(); i++)
320 res(i, j).vi_ = baseVari->_variRefC[pos++];
325 template <
int R1,
int C1,
int R2,
int C2>
327 Eigen::Matrix<var, R1, C2>
329 const Eigen::Matrix<var, R2, C2> &b) {
330 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
341 mdivide_left_dv_vari<R1, C1, R2, C2> *baseVari
342 =
new mdivide_left_dv_vari<R1, C1, R2, C2>(A, b);
345 for (
size_type j = 0; j < res.cols(); j++)
346 for (
size_type i = 0; i < res.rows(); i++)
347 res(i, j).vi_ = baseVari->_variRefC[pos++];
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)
Return the number of rows in the specified matrix, vector, or row vector.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double e.
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.
int cols(const Eigen::Matrix< T, R, C > &m)
Return the number of columns in the specified matrix, vector, or row vector.
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
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.
AutodiffStackStorage< vari, chainable_alloc > ChainableStack