1 #ifndef STAN_MATH_REV_MAT_FUN_TRACE_INV_QUAD_FORM_LDLT_HPP 2 #define STAN_MATH_REV_MAT_FUN_TRACE_INV_QUAD_FORM_LDLT_HPP 10 #include <type_traits> 16 template <
typename T2,
int R2,
int C2,
typename T3,
int R3,
int C3>
19 inline void initializeB(
const Eigen::Matrix<var, R3, C3> &B,
bool haveD) {
20 Eigen::Matrix<double, R3, C3> Bd(B.rows(), B.cols());
21 variB_.resize(B.rows(), B.cols());
22 for (
int j = 0; j < B.cols(); j++) {
23 for (
int i = 0; i < B.rows(); i++) {
24 variB_(i, j) = B(i, j).vi_;
25 Bd(i, j) = B(i, j).val();
30 C_.noalias() = Bd.transpose() *
AinvB_;
34 inline void initializeB(
const Eigen::Matrix<double, R3, C3> &B,
bool haveD) {
37 C_.noalias() = B.transpose() *
AinvB_;
42 template <
int R1,
int C1>
43 inline void initializeD(
const Eigen::Matrix<var, R1, C1> &D) {
44 D_.resize(D.rows(), D.cols());
45 variD_.resize(D.rows(), D.cols());
46 for (
int j = 0; j < D.cols(); j++) {
47 for (
int i = 0; i < D.rows(); i++) {
48 variD_(i, j) = D(i, j).vi_;
49 D_(i, j) = D(i, j).val();
53 template <
int R1,
int C1>
54 inline void initializeD(
const Eigen::Matrix<double, R1, C1> &D) {
59 template <
typename T1,
int R1,
int C1>
62 const Eigen::Matrix<T3, R3, C3> &B)
71 const Eigen::Matrix<T3, R3, C3> &B)
78 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
D_;
79 Eigen::Matrix<vari *, Eigen::Dynamic, Eigen::Dynamic>
variD_;
80 Eigen::Matrix<vari *, R3, C3>
variB_;
81 Eigen::Matrix<double, R3, C3>
AinvB_;
82 Eigen::Matrix<double, C3, C3>
C_;
86 template <
typename T2,
int R2,
int C2,
typename T3,
int R3,
int C3>
99 Eigen::Matrix<double, R2, C2> aA;
104 * (impl->
AinvB_ * impl->
D_.transpose() * impl->
AinvB_.transpose());
106 aA.noalias() = -adj * (impl->
AinvB_ * impl->
AinvB_.transpose());
108 for (
int j = 0; j < aA.cols(); j++)
109 for (
int i = 0; i < aA.rows(); i++)
110 impl->
ldlt_.alloc_->variA_(i, j)->adj_ += aA(i, j);
115 Eigen::Matrix<double, R3, C3> aB;
118 aB.noalias() = adj * impl->
AinvB_ * (impl->
D_ + impl->
D_.transpose());
120 aB.noalias() = 2 * adj * impl->
AinvB_;
122 for (
int j = 0; j < aB.cols(); j++)
123 for (
int i = 0; i < aB.rows(); i++)
124 impl->
variB_(i, j)->adj_ += aB(i, j);
141 if (impl_->Dtype_ == 1) {
142 for (
int j = 0; j < impl_->variD_.cols(); j++)
143 for (
int i = 0; i < impl_->variD_.rows(); i++)
144 impl_->variD_(i, j)->adj_ += adj_ * impl_->C_(i, j);
158 template <
typename T2,
int R2,
int C2,
typename T3,
int R3,
int C3>
163 const Eigen::Matrix<T3, R3, C3> &B) {
Defines a public enum named value which is defined to be false as the primitive scalar types cannot b...
The variable implementation base class.
const Eigen::internal::solve_retval< ldlt_t, Rhs > solve(const Eigen::MatrixBase< Rhs > &b) const
Independent (input) and dependent (output) variables for gradients.
std::enable_if< !stan::is_var< T1 >::value &&!stan::is_var< T2 >::value, typename boost::math::tools::promote_args< T1, T2 >::type >::type trace_inv_quad_form_ldlt(const LDLT_factor< T1, R2, C2 > &A, const Eigen::Matrix< T2, R3, C3 > &B)
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.
T trace(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Returns the trace of the specified matrix.
A chainable_alloc is an object which is constructed and destructed normally but the memory lifespan i...