chol2inv.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_CHOL2INV_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_CHOL2INV_HPP
3 
11 
12 namespace stan {
13 namespace math {
14 
15 /**
16  * Returns the inverse of the matrix whose Cholesky factor is L
17  * @tparam T The scalar type of the matrix
18  * @param L Matrix that is a Cholesky factor.
19  * @return The matrix inverse of L * L'
20  * @throw std::domain_error If the input matrix is not square or
21  * lower triangular
22  */
23 template <typename T>
24 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> chol2inv(
25  const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& L) {
26  check_square("chol2inv", "L", L);
27  check_lower_triangular("chol2inv", "L", L);
28  int K = L.rows();
29  typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
30  if (K == 0)
31  return L;
32  if (K == 1) {
33  matrix_t X(1, 1);
34  X.coeffRef(0) = inv_square(L.coeff(0));
35  return X;
36  }
37  matrix_t L_inv = mdivide_left_tri_low(L, matrix_t::Identity(K, K).eval());
38  matrix_t X(K, K);
39  for (int k = 0; k < K; ++k) {
40  X.coeffRef(k, k) = dot_self(L_inv.col(k).tail(K - k).eval());
41  for (int j = k + 1; j < K; ++j) {
42  int Kmj = K - j;
43  X.coeffRef(k, j) = X.coeffRef(j, k) = dot_product(
44  L_inv.col(k).tail(Kmj).eval(), L_inv.col(j).tail(Kmj).eval());
45  }
46  }
47  return X;
48 }
49 
50 } // namespace math
51 } // namespace stan
52 #endif
void check_lower_triangular(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > chol2inv(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &L)
Definition: chol2inv.hpp:24
fvar< T > dot_self(const Eigen::Matrix< fvar< T >, R, C > &v)
Definition: dot_self.hpp:15
Double_t K
static constexpr double L
const double j
Definition: BetheBloch.cxx:29
fvar< T > dot_product(const Eigen::Matrix< fvar< T >, R1, C1 > &v1, const Eigen::Matrix< fvar< T >, R2, C2 > &v2)
Definition: dot_product.hpp:16
fvar< T > inv_square(const fvar< T > &x)
Definition: inv_square.hpp:12
void check_square(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Eigen::Matrix< fvar< T >, R1, C1 > mdivide_left_tri_low(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
Float_t X
Definition: plot.C:38