ad embedded expression tests to expect_ad testing framework

stan/math/fwd/fun/mdivide_left_ldlt.hpp

@@ -32,16 +32,17 @@ mdivide_left_ldlt(LDLT_factor<T>& A, const EigMat& b) {
   check_multiplicable("mdivide_left_ldlt", "A", A.matrix(), "b", b);
 
   const auto& b_ref = to_ref(b);
-  Eigen::Matrix<EigMatValueScalar, R2, C2> b_val(b.rows(), b.cols());
-  Eigen::Matrix<EigMatValueScalar, R2, C2> b_der(b.rows(), b.cols());
-  for (int j = 0; j < b.cols(); j++) {
-    for (int i = 0; i < b.rows(); i++) {
+  Eigen::Matrix<EigMatValueScalar, R2, C2> b_val(b_ref.rows(), b_ref.cols());
+  Eigen::Matrix<EigMatValueScalar, R2, C2> b_der(b_ref.rows(), b_ref.cols());
+  for (int j = 0; j < b_ref.cols(); j++) {
+    for (int i = 0; i < b_ref.rows(); i++) {
       b_val.coeffRef(i, j) = b_ref.coeff(i, j).val_;
       b_der.coeffRef(i, j) = b_ref.coeff(i, j).d_;
     }
   }
-
-  return to_fvar(mdivide_left_ldlt(A, b_val), mdivide_left_ldlt(A, b_der));
+  auto&& A_ref = to_ref(A);
+  return to_fvar(mdivide_left_ldlt(A_ref, b_val),
+                 mdivide_left_ldlt(A_ref, b_der));
 }
 
 }  // namespace math

stan/math/prim/constraint/cov_matrix_constrain_lkj.hpp

@@ -32,9 +32,9 @@ namespace math {
  */
 template <typename T, require_eigen_vector_t<T>* = nullptr>
 inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-cov_matrix_constrain_lkj(const T& x, size_t k) {
+cov_matrix_constrain_lkj(T&& x, size_t k) {
   size_t k_choose_2 = (k * (k - 1)) / 2;
-  const auto& x_ref = to_ref(x);
+  auto&& x_ref = to_ref(std::forward<T>(x));
   return read_cov_matrix(corr_constrain(x_ref.head(k_choose_2)),
                          positive_constrain(x_ref.tail(k)));
 }
@@ -56,9 +56,9 @@ cov_matrix_constrain_lkj(const T& x, size_t k) {
  */
 template <typename T, require_eigen_vector_t<T>* = nullptr>
 inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-cov_matrix_constrain_lkj(const T& x, size_t k, return_type_t<T>& lp) {
+cov_matrix_constrain_lkj(T&& x, size_t k, return_type_t<T>& lp) {
   size_t k_choose_2 = (k * (k - 1)) / 2;
-  const auto& x_ref = x;
+  auto&& x_ref = to_ref(std::forward<T>(x));
   return read_cov_matrix(corr_constrain(x_ref.head(k_choose_2)),
                          positive_constrain(x_ref.tail(k)), lp);
 }

stan/math/prim/constraint/lb_constrain.hpp

@@ -152,8 +152,9 @@ inline auto lb_constrain(const T& x, const L& lb, return_type_t<T, L>& lp) {
 template <typename T, typename L, require_not_std_vector_t<L>* = nullptr>
 inline auto lb_constrain(const std::vector<T>& x, const L& lb) {
   std::vector<plain_type_t<decltype(lb_constrain(x[0], lb))>> ret(x.size());
+  auto&& lb_ref = to_ref(lb);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lb_constrain(x[i], lb);
+    ret[i] = lb_constrain(x[i], lb_ref);
   }
   return ret;
 }
@@ -173,8 +174,9 @@ template <typename T, typename L, require_not_std_vector_t<L>* = nullptr>
 inline auto lb_constrain(const std::vector<T>& x, const L& lb,
                          return_type_t<T, L>& lp) {
   std::vector<plain_type_t<decltype(lb_constrain(x[0], lb))>> ret(x.size());
+  auto&& lb_ref = to_ref(lb);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lb_constrain(x[i], lb, lp);
+    ret[i] = lb_constrain(x[i], lb_ref, lp);
   }
   return ret;
 }

stan/math/prim/constraint/lub_constrain.hpp

@@ -119,8 +119,8 @@ template <typename T, typename L, typename U, require_eigen_t<T>* = nullptr,
           require_all_stan_scalar_t<L, U>* = nullptr,
           require_not_var_t<return_type_t<T, L, U>>* = nullptr>
 inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
-  return eval(
-      x.unaryExpr([ub, lb](auto&& xx) { return lub_constrain(xx, lb, ub); }));
+  return eval(to_ref(x).unaryExpr(
+      [ub, lb](auto&& xx) { return lub_constrain(xx, lb, ub); }));
 }
 
 /**
@@ -131,7 +131,7 @@ template <typename T, typename L, typename U, require_eigen_t<T>* = nullptr,
           require_not_var_t<return_type_t<T, L, U>>* = nullptr>
 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
                           return_type_t<T, L, U>& lp) {
-  return eval(x.unaryExpr(
+  return eval(to_ref(x).unaryExpr(
       [lb, ub, &lp](auto&& xx) { return lub_constrain(xx, lb, ub, lp); }));
 }
 
@@ -145,8 +145,9 @@ template <typename T, typename L, typename U,
           require_not_var_t<return_type_t<T, L, U>>* = nullptr>
 inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
   check_matching_dims("lub_constrain", "x", x, "lb", lb);
-  return eval(x.binaryExpr(
-      lb, [ub](auto&& x, auto&& lb) { return lub_constrain(x, lb, ub); }));
+  return eval(to_ref(x).binaryExpr(to_ref(lb), [ub](auto&& x, auto&& lb) {
+    return lub_constrain(x, lb, ub);
+  }));
 }
 
 /**
@@ -160,7 +161,7 @@ template <typename T, typename L, typename U,
 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
                           return_type_t<T, L, U>& lp) {
   check_matching_dims("lub_constrain", "x", x, "lb", lb);
-  return eval(x.binaryExpr(lb, [ub, &lp](auto&& x, auto&& lb) {
+  return eval(to_ref(x).binaryExpr(to_ref(lb), [ub, &lp](auto&& x, auto&& lb) {
     return lub_constrain(x, lb, ub, lp);
   }));
 }
@@ -175,8 +176,9 @@ template <typename T, typename L, typename U,
           require_not_var_t<return_type_t<T, L, U>>* = nullptr>
 inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
   check_matching_dims("lub_constrain", "x", x, "ub", ub);
-  return eval(x.binaryExpr(
-      ub, [lb](auto&& x, auto&& ub) { return lub_constrain(x, lb, ub); }));
+  return eval(to_ref(x).binaryExpr(to_ref(ub), [lb](auto&& x, auto&& ub) {
+    return lub_constrain(x, lb, ub);
+  }));
 }
 
 /**
@@ -190,7 +192,7 @@ template <typename T, typename L, typename U,
 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
                           return_type_t<T, L, U>& lp) {
   check_matching_dims("lub_constrain", "x", x, "ub", ub);
-  return eval(x.binaryExpr(ub, [lb, &lp](auto&& x, auto&& ub) {
+  return eval(to_ref(x).binaryExpr(to_ref(ub), [lb, &lp](auto&& x, auto&& ub) {
     return lub_constrain(x, lb, ub, lp);
   }));
 }
@@ -248,8 +250,10 @@ template <typename T, typename L, typename U,
 inline auto lub_constrain(const std::vector<T>& x, const L& lb, const U& ub) {
   std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub))>> ret(
       x.size());
+  auto&& lb_ref = to_ref(lb);
+  auto&& ub_ref = to_ref(ub);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lub_constrain(x[i], lb, ub);
+    ret[i] = lub_constrain(x[i], lb_ref, ub_ref);
   }
   return ret;
 }
@@ -263,8 +267,10 @@ inline auto lub_constrain(const std::vector<T>& x, const L& lb, const U& ub,
                           return_type_t<T, L, U>& lp) {
   std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub))>> ret(
       x.size());
+  auto&& lb_ref = to_ref(lb);
+  auto&& ub_ref = to_ref(ub);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lub_constrain(x[i], lb, ub, lp);
+    ret[i] = lub_constrain(x[i], lb_ref, ub_ref, lp);
   }
   return ret;
 }
@@ -279,8 +285,9 @@ inline auto lub_constrain(const std::vector<T>& x, const L& lb,
   check_matching_dims("lub_constrain", "x", x, "ub", ub);
   std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub[0]))>> ret(
       x.size());
+  auto&& lb_ref = to_ref(lb);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lub_constrain(x[i], lb, ub[i]);
+    ret[i] = lub_constrain(x[i], lb_ref, ub[i]);
   }
   return ret;
 }
@@ -296,8 +303,9 @@ inline auto lub_constrain(const std::vector<T>& x, const L& lb,
   check_matching_dims("lub_constrain", "x", x, "ub", ub);
   std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub[0]))>> ret(
       x.size());
+  auto&& lb_ref = to_ref(lb);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lub_constrain(x[i], lb, ub[i], lp);
+    ret[i] = lub_constrain(x[i], lb_ref, ub[i], lp);
   }
   return ret;
 }
@@ -312,8 +320,9 @@ inline auto lub_constrain(const std::vector<T>& x, const std::vector<L>& lb,
   check_matching_dims("lub_constrain", "x", x, "lb", lb);
   std::vector<plain_type_t<decltype(lub_constrain(x[0], lb[0], ub))>> ret(
       x.size());
+  auto&& ub_ref = to_ref(ub);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lub_constrain(x[i], lb[i], ub);
+    ret[i] = lub_constrain(x[i], lb[i], ub_ref);
   }
   return ret;
 }
@@ -328,8 +337,9 @@ inline auto lub_constrain(const std::vector<T>& x, const std::vector<L>& lb,
   check_matching_dims("lub_constrain", "x", x, "lb", lb);
   std::vector<plain_type_t<decltype(lub_constrain(x[0], lb[0], ub))>> ret(
       x.size());
+  auto&& ub_ref = to_ref(ub);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lub_constrain(x[i], lb[i], ub, lp);
+    ret[i] = lub_constrain(x[i], lb[i], ub_ref, lp);
   }
   return ret;
 }

stan/math/prim/constraint/ub_constrain.hpp

@@ -162,8 +162,9 @@ inline auto ub_constrain(const T& x, const U& ub,
 template <typename T, typename U, require_not_std_vector_t<U>* = nullptr>
 inline auto ub_constrain(const std::vector<T>& x, const U& ub) {
   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub))>> ret(x.size());
+  auto&& ub_ref = to_ref(ub);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = ub_constrain(x[i], ub);
+    ret[i] = ub_constrain(x[i], ub_ref);
   }
   return ret;
 }
@@ -183,8 +184,9 @@ template <typename T, typename U, require_not_std_vector_t<U>* = nullptr>
 inline auto ub_constrain(const std::vector<T>& x, const U& ub,
                          return_type_t<T, U>& lp) {
   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub))>> ret(x.size());
+  auto&& ub_ref = to_ref(ub);
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = ub_constrain(x[i], ub, lp);
+    ret[i] = ub_constrain(x[i], ub_ref, lp);
   }
   return ret;
 }

stan/math/prim/fun/diag_post_multiply.hpp

@@ -25,7 +25,7 @@ template <typename T1, typename T2, require_eigen_t<T1>* = nullptr,
 auto diag_post_multiply(const T1& m1, const T2& m2) {
   check_size_match("diag_post_multiply", "m2.size()", m2.size(), "m1.cols()",
                    m1.cols());
-  return m1 * m2.asDiagonal();
+  return to_ref(to_ref(m1) * to_ref(m2).asDiagonal());
 }
 
 }  // namespace math

stan/math/prim/fun/diag_pre_multiply.hpp

@@ -26,7 +26,7 @@ auto diag_pre_multiply(const T1& m1, const T2& m2) {
   check_size_match("diag_pre_multiply", "m1.size()", m1.size(), "m2.rows()",
                    m2.rows());
 
-  return m1.asDiagonal() * m2;
+  return to_ref(to_ref(m1).asDiagonal() * to_ref(m2));
 }
 
 }  // namespace math

stan/math/prim/fun/fma.hpp

@@ -44,7 +44,8 @@ inline auto fma(T1&& x, T2&& y, T3&& z) {
       [](auto&& x, auto&& y, auto&& z) {
         return ((as_array_or_scalar(x) * as_array_or_scalar(y))
                 + as_array_or_scalar(z))
-            .matrix();
+            .matrix()
+            .eval();
       },
       std::forward<T1>(x), std::forward<T2>(y), std::forward<T3>(z));
 }

stan/math/prim/fun/matrix_exp.hpp

@@ -24,14 +24,14 @@ namespace math {
 template <typename T, typename = require_eigen_t<T>>
 inline plain_type_t<T> matrix_exp(const T& A_in) {
   using std::exp;
-  const auto& A = A_in.eval();
+  const auto& A = to_ref(A_in);
   check_square("matrix_exp", "input matrix", A);
   if (T::RowsAtCompileTime == 1 && T::ColsAtCompileTime == 1) {
     plain_type_t<T> res;
     res << exp(A(0));
     return res;
   }
-  if (A_in.size() == 0) {
+  if (A.size() == 0) {
     return {};
   }
   return (A.cols() == 2

stan/math/prim/fun/matrix_exp_2x2.hpp

@@ -28,9 +28,9 @@ matrix_exp_2x2(const EigMat& A) {
   using std::exp;
   using std::sinh;
   using std::sqrt;
-
+  auto&& A_ref = to_ref(A);
   using T = value_type_t<EigMat>;
-  T a = A(0, 0), b = A(0, 1), c = A(1, 0), d = A(1, 1), delta;
+  T a = A_ref(0, 0), b = A_ref(0, 1), c = A_ref(1, 0), d = A_ref(1, 1), delta;
   delta = sqrt(square(a - d) + 4 * b * c);
 
   Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> B(2, 2);
@@ -49,7 +49,7 @@ matrix_exp_2x2(const EigMat& A) {
 
   // use pade approximation if cosh & sinh ops overflow to NaN
   if (B.hasNaN()) {
-    return matrix_exp_pade(A);
+    return matrix_exp_pade(A_ref);
   } else {
     return B / delta;
   }

stan/math/prim/fun/matrix_exp_multiply.hpp

@@ -22,13 +22,15 @@ template <typename EigMat1, typename EigMat2,
           require_all_st_same<double, EigMat1, EigMat2>* = nullptr>
 inline Eigen::Matrix<double, Eigen::Dynamic, EigMat2::ColsAtCompileTime>
 matrix_exp_multiply(const EigMat1& A, const EigMat2& B) {
-  check_square("matrix_exp_multiply", "input matrix", A);
-  check_multiplicable("matrix_exp_multiply", "A", A, "B", B);
-  if (A.size() == 0) {
-    return {0, B.cols()};
+  auto&& A_ref = to_ref(A);
+  auto&& B_ref = to_ref(B);
+  check_square("matrix_exp_multiply", "input matrix", A_ref);
+  check_multiplicable("matrix_exp_multiply", "A", A_ref, "B", B_ref);
+  if (A_ref.size() == 0) {
+    return {0, B_ref.cols()};
   }
 
-  return matrix_exp_action_handler().action(A, B);
+  return matrix_exp_action_handler().action(A_ref, B_ref);
 }
 
 }  // namespace math

stan/math/prim/fun/matrix_exp_pade.hpp

@@ -21,18 +21,20 @@ template <typename EigMat, require_eigen_t<EigMat>* = nullptr>
 Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime,
               EigMat::ColsAtCompileTime>
 matrix_exp_pade(const EigMat& arg) {
+  auto&& arg_ref = to_ref(arg);
   using MatrixType
       = Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime,
                       EigMat::ColsAtCompileTime>;
-  check_square("matrix_exp_pade", "arg", arg);
-  if (arg.size() == 0) {
+  check_square("matrix_exp_pade", "arg", arg_ref);
+  if (arg_ref.size() == 0) {
     return {};
   }
 
   MatrixType U, V;
   int squarings;
 
-  Eigen::matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings, arg(0, 0));
+  Eigen::matrix_exp_computeUV<MatrixType>::run(arg_ref, U, V, squarings,
+                                               arg_ref(0, 0));
   // Pade approximant is
   // (U+V) / (-U+V)
   MatrixType numer = U + V;

stan/math/prim/fun/pseudo_eigenvalues.hpp

@@ -7,12 +7,13 @@
 namespace stan {
 namespace math {
 
-template <typename T>
-Eigen::Matrix<T, -1, -1> pseudo_eigenvalues(const Eigen::Matrix<T, -1, -1>& m) {
+template <typename EigMat, require_eigen_matrix_base_t<EigMat>* = nullptr>
+inline Eigen::Matrix<scalar_type_t<EigMat>, -1, -1> pseudo_eigenvalues(
+    EigMat&& m) {
   check_nonzero_size("pseudo_eigenvalues", "m", m);
   check_square("pseudo_eigenvalues", "m", m);
-
-  Eigen::EigenSolver<Eigen::Matrix<T, -1, -1>> solver(m);
+  Eigen::EigenSolver<Eigen::Matrix<scalar_type_t<EigMat>, -1, -1>> solver(
+      to_ref(std::forward<EigMat>(m)));
   return solver.pseudoEigenvalueMatrix();
 }
 

stan/math/prim/fun/pseudo_eigenvectors.hpp

@@ -7,13 +7,13 @@
 namespace stan {
 namespace math {
 
-template <typename T>
-Eigen::Matrix<T, -1, -1> pseudo_eigenvectors(
-    const Eigen::Matrix<T, -1, -1>& m) {
+template <typename EigMat, require_eigen_matrix_base_t<EigMat>* = nullptr>
+inline Eigen::Matrix<scalar_type_t<EigMat>, -1, -1> pseudo_eigenvectors(
+    EigMat&& m) {
   check_nonzero_size("pseudo_eigenvectors", "m", m);
   check_square("pseudo_eigenvectors", "m", m);
-
-  Eigen::EigenSolver<Eigen::Matrix<T, -1, -1>> solver(m);
+  Eigen::EigenSolver<Eigen::Matrix<scalar_type_t<EigMat>, -1, -1>> solver(
+      to_ref(std::forward<EigMat>(m)));
   return solver.pseudoEigenvectors();
 }