Automatic Differentiation
 
Loading...
Searching...
No Matches
mdivide_right_tri_low.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_FWD_FUN_MDIVIDE_RIGHT_TRI_LOW_HPP
2#define STAN_MATH_FWD_FUN_MDIVIDE_RIGHT_TRI_LOW_HPP
3
10
11namespace stan {
12namespace math {
13
14template <typename EigMat1, typename EigMat2,
15 require_all_eigen_vt<is_fvar, EigMat1, EigMat2>* = nullptr,
16 require_vt_same<EigMat1, EigMat2>* = nullptr>
17inline Eigen::Matrix<value_type_t<EigMat1>, EigMat1::RowsAtCompileTime,
18 EigMat2::ColsAtCompileTime>
19mdivide_right_tri_low(const EigMat1& A, const EigMat2& b) {
20 using T = typename value_type_t<EigMat1>::Scalar;
21 constexpr int R1 = EigMat1::RowsAtCompileTime;
22 constexpr int C1 = EigMat1::ColsAtCompileTime;
23 constexpr int R2 = EigMat2::RowsAtCompileTime;
24 constexpr int C2 = EigMat2::ColsAtCompileTime;
25
26 check_square("mdivide_right_tri_low", "b", b);
27 check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
28 if (b.size() == 0) {
29 return {A.rows(), 0};
30 }
31
32 Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
33 Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
34 Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
35 Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
36 val_b.setZero();
37 deriv_b.setZero();
38
39 const Eigen::Ref<const plain_type_t<EigMat1>>& A_ref = A;
40 for (int j = 0; j < A.cols(); j++) {
41 for (int i = 0; i < A.rows(); i++) {
42 val_A.coeffRef(i, j) = A_ref.coeff(i, j).val_;
43 deriv_A.coeffRef(i, j) = A_ref.coeff(i, j).d_;
44 }
45 }
46
47 const Eigen::Ref<const plain_type_t<EigMat2>>& b_ref = b;
48 for (int j = 0; j < b.cols(); j++) {
49 for (int i = j; i < b.rows(); i++) {
50 val_b.coeffRef(i, j) = b_ref.coeff(i, j).val_;
51 deriv_b.coeffRef(i, j) = b_ref.coeff(i, j).d_;
52 }
53 }
54
55 Eigen::Matrix<T, R1, C2> A_mult_inv_b = mdivide_right(val_A, val_b);
56 return to_fvar(A_mult_inv_b,
57 mdivide_right(deriv_A, val_b)
58 - A_mult_inv_b * mdivide_right(deriv_b, val_b));
59}
60
61template <typename EigMat1, typename EigMat2,
64inline Eigen::Matrix<value_type_t<EigMat1>, EigMat1::RowsAtCompileTime,
65 EigMat2::ColsAtCompileTime>
66mdivide_right_tri_low(const EigMat1& A, const EigMat2& b) {
67 using T = typename value_type_t<EigMat1>::Scalar;
68 constexpr int R1 = EigMat1::RowsAtCompileTime;
69 constexpr int C1 = EigMat1::ColsAtCompileTime;
70
71 check_square("mdivide_right_tri_low", "b", b);
72 check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
73 if (b.size() == 0) {
74 return {A.rows(), 0};
75 }
76
77 Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
78 Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
79
80 const Eigen::Ref<const plain_type_t<EigMat1>>& A_ref = A;
81 for (int j = 0; j < A.cols(); j++) {
82 for (int i = 0; i < A.rows(); i++) {
83 val_A.coeffRef(i, j) = A_ref.coeff(i, j).val_;
84 deriv_A.coeffRef(i, j) = A_ref.coeff(i, j).d_;
85 }
86 }
87
88 plain_type_t<EigMat2> val_b = b.template triangularView<Eigen::Lower>();
89
90 return to_fvar(mdivide_right(val_A, val_b), mdivide_right(deriv_A, val_b));
91}
92
93template <typename EigMat1, typename EigMat2,
94 require_eigen_vt<std::is_arithmetic, EigMat1>* = nullptr,
95 require_eigen_vt<is_fvar, EigMat2>* = nullptr>
96inline Eigen::Matrix<value_type_t<EigMat2>, EigMat1::RowsAtCompileTime,
97 EigMat2::ColsAtCompileTime>
98mdivide_right_tri_low(const EigMat1& A, const EigMat2& b) {
99 using T = typename value_type_t<EigMat2>::Scalar;
100 constexpr int R1 = EigMat1::RowsAtCompileTime;
101 constexpr int R2 = EigMat2::RowsAtCompileTime;
102 constexpr int C2 = EigMat2::ColsAtCompileTime;
103 check_square("mdivide_right_tri_low", "b", b);
104 check_multiplicable("mdivide_right_tri_low", "A", A, "b", b);
105 if (b.size() == 0) {
106 return {A.rows(), 0};
107 }
108
109 Eigen::Matrix<T, R1, C2> A_mult_inv_b(A.rows(), b.cols());
110 Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
111 Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
112 val_b.setZero();
113 deriv_b.setZero();
114
115 for (int j = 0; j < b.cols(); j++) {
116 for (int i = j; i < b.rows(); i++) {
117 val_b(i, j) = b(i, j).val_;
118 deriv_b(i, j) = b(i, j).d_;
119 }
120 }
121
122 A_mult_inv_b = mdivide_right(A, val_b);
123
124 return to_fvar(A_mult_inv_b,
125 -multiply(A_mult_inv_b, mdivide_right(deriv_b, val_b)));
126}
127
128} // namespace math
129} // namespace stan
130#endif
require_t< container_type_check_base< is_eigen, value_type_t, TypeCheck, Check... > > require_eigen_vt
Require type satisfies is_eigen.
Definition is_eigen.hpp:97
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
void check_square(const char *function, const char *name, const T_y &y)
Check if the specified matrix is square.
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< value_type_t< EigMat1 >, EigMat1::RowsAtCompileTime, EigMat2::ColsAtCompileTime > mdivide_right_tri_low(const EigMat1 &A, const EigMat2 &b)
auto multiply(const Mat1 &m1, const Mat2 &m2)
Return the product of the specified matrices.
Definition multiply.hpp:18
Eigen::Matrix< value_type_t< EigMat1 >, EigMat1::RowsAtCompileTime, EigMat2::ColsAtCompileTime > mdivide_right(const EigMat1 &A, const EigMat2 &b)
fvar< T > to_fvar(const T &x)
Definition to_fvar.hpp:15
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...