Automatic Differentiation
 
Loading...
Searching...
No Matches
mdivide_left_tri_low.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_FWD_FUN_MDIVIDE_LEFT_TRI_LOW_HPP
2#define STAN_MATH_FWD_FUN_MDIVIDE_LEFT_TRI_LOW_HPP
3
11
12namespace stan {
13namespace math {
14
15template <typename T1, typename T2,
16 require_all_eigen_vt<is_fvar, T1, T2>* = nullptr,
17 require_vt_same<T1, T2>* = nullptr>
18inline Eigen::Matrix<value_type_t<T1>, T1::RowsAtCompileTime,
19 T2::ColsAtCompileTime>
20mdivide_left_tri_low(const T1& A, const T2& b) {
21 using T = typename value_type_t<T1>::Scalar;
22 constexpr int S1 = T1::RowsAtCompileTime;
23 constexpr int C2 = T2::ColsAtCompileTime;
24
25 check_square("mdivide_left_tri_low", "A", A);
26 check_multiplicable("mdivide_left_tri_low", "A", A, "b", b);
27 if (A.size() == 0) {
28 return {0, b.cols()};
29 }
30
31 Eigen::Matrix<T, S1, S1> val_A(A.rows(), A.cols());
32 Eigen::Matrix<T, S1, S1> deriv_A(A.rows(), A.cols());
33 val_A.setZero();
34 deriv_A.setZero();
35
36 const Eigen::Ref<const plain_type_t<T2>>& b_ref = b;
37 const Eigen::Ref<const plain_type_t<T1>>& A_ref = A;
38 for (size_type j = 0; j < A.cols(); j++) {
39 for (size_type i = j; i < A.rows(); i++) {
40 val_A(i, j) = A_ref(i, j).val_;
41 deriv_A(i, j) = A_ref(i, j).d_;
42 }
43 }
44
45 Eigen::Matrix<T, S1, C2> inv_A_mult_b = mdivide_left(val_A, b_ref.val());
46
47 return to_fvar(inv_A_mult_b,
48 mdivide_left(val_A, b_ref.d())
49 - multiply(mdivide_left(val_A, deriv_A), inv_A_mult_b));
50}
51
52template <typename T1, typename T2, require_eigen_t<T1>* = nullptr,
53 require_vt_same<double, T1>* = nullptr,
54 require_eigen_vt<is_fvar, T2>* = nullptr>
55inline Eigen::Matrix<value_type_t<T2>, T1::RowsAtCompileTime,
56 T2::ColsAtCompileTime>
57mdivide_left_tri_low(const T1& A, const T2& b) {
58 constexpr int S1 = T1::RowsAtCompileTime;
59
60 check_square("mdivide_left_tri_low", "A", A);
61 check_multiplicable("mdivide_left_tri_low", "A", A, "b", b);
62 if (A.size() == 0) {
63 return {0, b.cols()};
64 }
65
66 Eigen::Matrix<double, S1, S1> val_A(A.rows(), A.cols());
67 val_A.setZero();
68
69 const Eigen::Ref<const plain_type_t<T2>>& b_ref = b;
70 const Eigen::Ref<const plain_type_t<T1>>& A_ref = A;
71 for (size_type j = 0; j < A.cols(); j++) {
72 for (size_type i = j; i < A.rows(); i++) {
73 val_A(i, j) = A_ref(i, j);
74 }
75 }
76
77 return to_fvar(mdivide_left(val_A, b_ref.val()),
78 mdivide_left(val_A, b_ref.d()));
79}
80
81template <typename T1, typename T2, require_eigen_vt<is_fvar, T1>* = nullptr,
82 require_eigen_t<T2>* = nullptr,
83 require_vt_same<double, T2>* = nullptr>
84inline Eigen::Matrix<value_type_t<T1>, T1::RowsAtCompileTime,
85 T2::ColsAtCompileTime>
86mdivide_left_tri_low(const T1& A, const T2& b) {
87 using T = typename value_type_t<T1>::Scalar;
88 constexpr int S1 = T1::RowsAtCompileTime;
89 constexpr int C2 = T2::ColsAtCompileTime;
90
91 check_square("mdivide_left_tri_low", "A", A);
92 check_multiplicable("mdivide_left_tri_low", "A", A, "b", b);
93 if (A.size() == 0) {
94 return {0, b.cols()};
95 }
96
97 Eigen::Matrix<T, S1, S1> val_A(A.rows(), A.cols());
98 Eigen::Matrix<T, S1, S1> deriv_A(A.rows(), A.cols());
99 val_A.setZero();
100 deriv_A.setZero();
101
102 const Eigen::Ref<const plain_type_t<T1>>& A_ref = A;
103 for (size_type j = 0; j < A.cols(); j++) {
104 for (size_type i = j; i < A.rows(); i++) {
105 val_A(i, j) = A_ref(i, j).val_;
106 deriv_A(i, j) = A_ref(i, j).d_;
107 }
108 }
109
110 Eigen::Matrix<T, S1, C2> inv_A_mult_b = mdivide_left(val_A, b);
111
112 return to_fvar(inv_A_mult_b,
113 -multiply(mdivide_left(val_A, deriv_A), inv_A_mult_b));
114}
115
116} // namespace math
117} // namespace stan
118#endif
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.
auto multiply(const Mat1 &m1, const Mat2 &m2)
Return the product of the specified matrices.
Definition multiply.hpp:18
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double elements.
Definition typedefs.hpp:11
fvar< T > to_fvar(const T &x)
Definition to_fvar.hpp:15
Eigen::Matrix< value_type_t< T1 >, T1::RowsAtCompileTime, T2::ColsAtCompileTime > mdivide_left(const T1 &A, const T2 &b)
Eigen::Matrix< value_type_t< T1 >, T1::RowsAtCompileTime, T2::ColsAtCompileTime > mdivide_left_tri_low(const T1 &A, const T2 &b)
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...