1 #ifndef STAN_MATH_PRIM_FUN_MULTIPLY_LOWER_TRI_SELF_TRANSPOSE_HPP 2 #define STAN_MATH_PRIM_FUN_MULTIPLY_LOWER_TRI_SELF_TRANSPOSE_HPP 3 4 #include <stan/math/prim/fun/typedefs.hpp> 5 #include <stan/math/prim/fun/square.hpp> 6 7 namespace stan { 8 namespace math { 9 10 /** 11 * Returns the result of multiplying the lower triangular 12 * portion of the input matrix by its own transpose. 13 * 14 * @param L Matrix to multiply. 15 * @return The lower triangular values in L times their own 16 * transpose. 17 */ 18 template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr, 19 require_not_st_autodiff<EigMat>* = nullptr> multiply_lower_tri_self_transpose(const EigMat & L)20inline matrix_d multiply_lower_tri_self_transpose(const EigMat& L) { 21 int K = L.rows(); 22 if (K == 0) { 23 return L; 24 } 25 if (K == 1) { 26 matrix_d result(1, 1); 27 result.coeffRef(0) = square(L.coeff(0, 0)); 28 return result; 29 } 30 int J = L.cols(); 31 matrix_d LLt(K, K); 32 matrix_d Lt = L.transpose(); 33 for (int m = 0; m < K; ++m) { 34 int k = (J < m + 1) ? J : m + 1; 35 LLt(m, m) = Lt.col(m).head(k).squaredNorm(); 36 for (int n = (m + 1); n < K; ++n) { 37 LLt(n, m) = LLt(m, n) = Lt.col(m).head(k).dot(Lt.col(n).head(k)); 38 } 39 } 40 return LLt; 41 } 42 43 } // namespace math 44 } // namespace stan 45 46 #endif 47