Stan Math Library  2.20.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 
10 #include <stan/math/fwd/core.hpp>
11 
12 namespace stan {
13 namespace math {
14 
15 template <typename T, int R1, int C1, int R2, int C2>
16 inline Eigen::Matrix<fvar<T>, R1, C1> mdivide_right_tri_low(
17  const Eigen::Matrix<fvar<T>, R1, C1> &A,
18  const Eigen::Matrix<fvar<T>, R2, C2> &b) {
19  check_square("mdivide_right_tri_low", "b", b);
20  check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
21 
22  Eigen::Matrix<T, R1, C2> A_mult_inv_b(A.rows(), b.cols());
23  Eigen::Matrix<T, R1, C2> deriv_A_mult_inv_b(A.rows(), b.cols());
24  Eigen::Matrix<T, R2, C2> deriv_b_mult_inv_b(b.rows(), b.cols());
25  Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
26  Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
27  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
28  Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
29  val_b.setZero();
30  deriv_b.setZero();
31 
32  for (size_type j = 0; j < A.cols(); j++) {
33  for (size_type i = 0; i < A.rows(); i++) {
34  val_A(i, j) = A(i, j).val_;
35  deriv_A(i, j) = A(i, j).d_;
36  }
37  }
38 
39  for (size_type j = 0; j < b.cols(); j++) {
40  for (size_type i = j; i < b.rows(); i++) {
41  val_b(i, j) = b(i, j).val_;
42  deriv_b(i, j) = b(i, j).d_;
43  }
44  }
45 
46  A_mult_inv_b = mdivide_right(val_A, val_b);
47  deriv_A_mult_inv_b = mdivide_right(deriv_A, val_b);
48  deriv_b_mult_inv_b = mdivide_right(deriv_b, val_b);
49 
50  Eigen::Matrix<T, R1, C2> deriv(A.rows(), b.cols());
51  deriv = deriv_A_mult_inv_b - multiply(A_mult_inv_b, deriv_b_mult_inv_b);
52 
53  return to_fvar(A_mult_inv_b, deriv);
54 }
55 
56 template <typename T, int R1, int C1, int R2, int C2>
57 inline Eigen::Matrix<fvar<T>, R1, C2> mdivide_right_tri_low(
58  const Eigen::Matrix<fvar<T>, R1, C1> &A,
59  const Eigen::Matrix<double, R2, C2> &b) {
60  check_square("mdivide_right_tri_low", "b", b);
61  check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
62 
63  Eigen::Matrix<T, R2, C2> deriv_b_mult_inv_b(b.rows(), b.cols());
64  Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
65  Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
66  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
67  val_b.setZero();
68 
69  for (int j = 0; j < A.cols(); j++) {
70  for (int i = 0; i < A.rows(); i++) {
71  val_A(i, j) = A(i, j).val_;
72  deriv_A(i, j) = A(i, j).d_;
73  }
74  }
75 
76  for (size_type j = 0; j < b.cols(); j++) {
77  for (size_type i = j; i < b.rows(); i++) {
78  val_b(i, j) = b(i, j);
79  }
80  }
81 
82  return to_fvar(mdivide_right(val_A, val_b), mdivide_right(deriv_A, val_b));
83 }
84 
85 template <typename T, int R1, int C1, int R2, int C2>
86 inline Eigen::Matrix<fvar<T>, R1, C2> mdivide_right_tri_low(
87  const Eigen::Matrix<double, R1, C1> &A,
88  const Eigen::Matrix<fvar<T>, R2, C2> &b) {
89  check_square("mdivide_right_tri_low", "b", b);
90  check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
91 
92  Eigen::Matrix<T, R1, C2> A_mult_inv_b(A.rows(), b.cols());
93  Eigen::Matrix<T, R2, C2> deriv_b_mult_inv_b(b.rows(), b.cols());
94  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
95  Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
96  val_b.setZero();
97  deriv_b.setZero();
98 
99  for (int j = 0; j < b.cols(); j++) {
100  for (int i = j; i < b.rows(); i++) {
101  val_b(i, j) = b(i, j).val_;
102  deriv_b(i, j) = b(i, j).d_;
103  }
104  }
105 
106  A_mult_inv_b = mdivide_right(A, val_b);
107  deriv_b_mult_inv_b = mdivide_right(deriv_b, val_b);
108 
109  Eigen::Matrix<T, R1, C2> deriv(A.rows(), b.cols());
110  deriv = -multiply(A_mult_inv_b, deriv_b_mult_inv_b);
111 
112  return to_fvar(A_mult_inv_b, deriv);
113 }
114 
115 } // namespace math
116 } // namespace stan
117 #endif
void check_square(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is square.
Eigen::Matrix< fvar< T >, R1, C1 > multiply(const Eigen::Matrix< fvar< T >, R1, C1 > &m, const fvar< T > &c)
Definition: multiply.hpp:14
std::vector< fvar< T > > to_fvar(const std::vector< T > &v)
Definition: to_fvar.hpp:12
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
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.
Eigen::Matrix< fvar< T >, R1, C2 > mdivide_right(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
This template class represents scalars used in forward-mode automatic differentiation, which consist of values and directional derivatives of the specified template type.
Definition: fvar.hpp:41
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–2018, Stan Development Team.