Stan Math Library  2.12.0
reverse mode automatic differentiation
mdivide_left_tri_low.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_FWD_MAT_FUN_MDIVIDE_LEFT_TRI_LOW_HPP
2 #define STAN_MATH_FWD_MAT_FUN_MDIVIDE_LEFT_TRI_LOW_HPP
3 
12 #include <stan/math/fwd/core.hpp>
14 #include <vector>
15 
16 namespace stan {
17  namespace math {
18 
19  template<typename T, int R1, int C1, int R2, int C2>
20  inline
21  Eigen::Matrix<fvar<T>, R1, C1>
22  mdivide_left_tri_low(const Eigen::Matrix<fvar<T>, R1, C1>& A,
23  const Eigen::Matrix<fvar<T>, R2, C2>& b) {
24  check_square("mdivide_left_tri_low", "A", A);
25  check_multiplicable("mdivide_left_tri_low", "A", A, "b", b);
26 
27  Eigen::Matrix<T, R1, C2> inv_A_mult_b(A.rows(), b.cols());
28  Eigen::Matrix<T, R1, C2> inv_A_mult_deriv_b(A.rows(), b.cols());
29  Eigen::Matrix<T, R1, C1> inv_A_mult_deriv_A(A.rows(), A.cols());
30  Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
31  Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
32  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
33  Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
34  val_A.setZero();
35  deriv_A.setZero();
36 
37  for (size_type j = 0; j < A.cols(); j++) {
38  for (size_type i = j; i < A.rows(); i++) {
39  val_A(i, j) = A(i, j).val_;
40  deriv_A(i, j) = A(i, j).d_;
41  }
42  }
43 
44  for (size_type j = 0; j < b.cols(); j++) {
45  for (size_type i = 0; i < b.rows(); i++) {
46  val_b(i, j) = b(i, j).val_;
47  deriv_b(i, j) = b(i, j).d_;
48  }
49  }
50 
51  inv_A_mult_b = mdivide_left(val_A, val_b);
52  inv_A_mult_deriv_b = mdivide_left(val_A, deriv_b);
53  inv_A_mult_deriv_A = mdivide_left(val_A, deriv_A);
54 
55  Eigen::Matrix<T, R1, C2> deriv(A.rows(), b.cols());
56  deriv = inv_A_mult_deriv_b - multiply(inv_A_mult_deriv_A, inv_A_mult_b);
57 
58  return to_fvar(inv_A_mult_b, deriv);
59  }
60 
61  template<typename T, int R1, int C1, int R2, int C2>
62  inline
63  Eigen::Matrix<fvar<T>, R1, C1>
64  mdivide_left_tri_low(const Eigen::Matrix<double, R1, C1>& A,
65  const Eigen::Matrix<fvar<T>, R2, C2>& b) {
66  check_square("mdivide_left_tri_low", "A", A);
67  check_multiplicable("mdivide_left_tri_low", "A", A, "b", b);
68 
69  Eigen::Matrix<T, R1, C2> inv_A_mult_b(A.rows(), b.cols());
70  Eigen::Matrix<T, R1, C2> inv_A_mult_deriv_b(A.rows(), b.cols());
71  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
72  Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
73  Eigen::Matrix<double, R1, C1> val_A(A.rows(), A.cols());
74  val_A.setZero();
75 
76  for (size_type j = 0; j < A.cols(); j++) {
77  for (size_type i = j; i < A.rows(); i++) {
78  val_A(i, j) = A(i, j);
79  }
80  }
81 
82  for (size_type j = 0; j < b.cols(); j++) {
83  for (size_type i = 0; i < b.rows(); i++) {
84  val_b(i, j) = b(i, j).val_;
85  deriv_b(i, j) = b(i, j).d_;
86  }
87  }
88 
89  inv_A_mult_b = mdivide_left(val_A, val_b);
90  inv_A_mult_deriv_b = mdivide_left(val_A, deriv_b);
91 
92  Eigen::Matrix<T, R1, C2> deriv(A.rows(), b.cols());
93  deriv = inv_A_mult_deriv_b;
94 
95  return to_fvar(inv_A_mult_b, deriv);
96  }
97 
98  template<typename T, int R1, int C1, int R2, int C2>
99  inline
100  Eigen::Matrix<fvar<T>, R1, C1>
101  mdivide_left_tri_low(const Eigen::Matrix<fvar<T>, R1, C1>& A,
102  const Eigen::Matrix<double, R2, C2>& b) {
103  check_square("mdivide_left_tri_low", "A", A);
104  check_multiplicable("mdivide_left_tri_low", "A", A, "b", b);
105 
106  Eigen::Matrix<T, R1, C2> inv_A_mult_b(A.rows(), b.cols());
107  Eigen::Matrix<T, R1, C1> inv_A_mult_deriv_A(A.rows(), A.cols());
108  Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
109  Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
110  val_A.setZero();
111  deriv_A.setZero();
112 
113  for (size_type j = 0; j < A.cols(); j++) {
114  for (size_type i = j; i < A.rows(); i++) {
115  val_A(i, j) = A(i, j).val_;
116  deriv_A(i, j) = A(i, j).d_;
117  }
118  }
119 
120  inv_A_mult_b = mdivide_left(val_A, b);
121  inv_A_mult_deriv_A = mdivide_left(val_A, deriv_A);
122 
123  Eigen::Matrix<T, R1, C2> deriv(A.rows(), b.cols());
124  deriv = -multiply(inv_A_mult_deriv_A, inv_A_mult_b);
125 
126  return to_fvar(inv_A_mult_b, deriv);
127  }
128 
129  }
130 }
131 #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)
Eigen::Matrix< fvar< T >, R1, C1 > multiply(const Eigen::Matrix< fvar< T >, R1, C1 > &m, const fvar< T > &c)
Definition: multiply.hpp:20
std::vector< fvar< T > > to_fvar(const std::vector< T > &v)
Definition: to_fvar.hpp:14
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
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.
Eigen::Matrix< fvar< T >, R1, C1 > mdivide_left_tri_low(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
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.

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