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)20 inline 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