Stan Math Library  2.12.0
reverse mode automatic differentiation
mdivide_right_tri_low.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_FWD_MAT_FUN_MDIVIDE_RIGHT_TRI_LOW_HPP
2 #define STAN_MATH_FWD_MAT_FUN_MDIVIDE_RIGHT_TRI_LOW_HPP
3 
13 #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_right_tri_low(const Eigen::Matrix<fvar<T>, R1, C1>& A,
23  const Eigen::Matrix<fvar<T>, R2, C2>& b) {
24  check_square("mdivide_right_tri_low", "b", b);
25  check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
26 
27  Eigen::Matrix<T, R1, C2> A_mult_inv_b(A.rows(), b.cols());
28  Eigen::Matrix<T, R1, C2> deriv_A_mult_inv_b(A.rows(), b.cols());
29  Eigen::Matrix<T, R2, C2> deriv_b_mult_inv_b(b.rows(), b.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_b.setZero();
35  deriv_b.setZero();
36 
37  for (size_type j = 0; j < A.cols(); j++) {
38  for (size_type i = 0; 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 = j; 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  A_mult_inv_b = mdivide_right(val_A, val_b);
52  deriv_A_mult_inv_b = mdivide_right(deriv_A, val_b);
53  deriv_b_mult_inv_b = mdivide_right(deriv_b, val_b);
54 
55  Eigen::Matrix<T, R1, C2> deriv(A.rows(), b.cols());
56  deriv = deriv_A_mult_inv_b - multiply(A_mult_inv_b, deriv_b_mult_inv_b);
57 
58  return to_fvar(A_mult_inv_b, deriv);
59  }
60 
61  template <typename T, int R1, int C1, int R2, int C2>
62  inline
63  Eigen::Matrix<fvar<T>, R1, C2>
64  mdivide_right_tri_low(const Eigen::Matrix<fvar<T>, R1, C1> &A,
65  const Eigen::Matrix<double, R2, C2> &b) {
66  check_square("mdivide_right_tri_low", "b", b);
67  check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
68 
69  Eigen::Matrix<T, R2, C2> deriv_b_mult_inv_b(b.rows(), b.cols());
70  Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
71  Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
72  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
73  val_b.setZero();
74 
75  for (int j = 0; j < A.cols(); j++) {
76  for (int i = 0; i < A.rows(); i++) {
77  val_A(i, j) = A(i, j).val_;
78  deriv_A(i, j) = A(i, j).d_;
79  }
80  }
81 
82  for (size_type j = 0; j < b.cols(); j++) {
83  for (size_type i = j; i < b.rows(); i++) {
84  val_b(i, j) = b(i, j);
85  }
86  }
87 
88  return to_fvar(mdivide_right(val_A, val_b),
89  mdivide_right(deriv_A, val_b));
90  }
91 
92  template <typename T, int R1, int C1, int R2, int C2>
93  inline
94  Eigen::Matrix<fvar<T>, R1, C2>
95  mdivide_right_tri_low(const Eigen::Matrix<double, R1, C1> &A,
96  const Eigen::Matrix<fvar<T>, R2, C2> &b) {
97  check_square("mdivide_right_tri_low", "b", b);
98  check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
99 
100  Eigen::Matrix<T, R1, C2>
101  A_mult_inv_b(A.rows(), b.cols());
102  Eigen::Matrix<T, R2, C2> deriv_b_mult_inv_b(b.rows(), b.cols());
103  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
104  Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
105  val_b.setZero();
106  deriv_b.setZero();
107 
108  for (int j = 0; j < b.cols(); j++) {
109  for (int i = j; i < b.rows(); i++) {
110  val_b(i, j) = b(i, j).val_;
111  deriv_b(i, j) = b(i, j).d_;
112  }
113  }
114 
115  A_mult_inv_b = mdivide_right(A, val_b);
116  deriv_b_mult_inv_b = mdivide_right(deriv_b, val_b);
117 
118  Eigen::Matrix<T, R1, C2>
119  deriv(A.rows(), b.cols());
120  deriv = -multiply(A_mult_inv_b, deriv_b_mult_inv_b);
121 
122  return to_fvar(A_mult_inv_b, deriv);
123  }
124 
125  }
126 }
127 #endif
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.
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.
Eigen::Matrix< fvar< T >, R1, C2 > mdivide_right(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
Eigen::Matrix< fvar< T >, R1, C1 > mdivide_right_tri_low(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)

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