diff --git a/stan/math/prim/prob/multi_normal_lpdf.hpp b/stan/math/prim/prob/multi_normal_lpdf.hpp
index 38f3f41da82..ec04aa25a40 100644
--- a/stan/math/prim/prob/multi_normal_lpdf.hpp
+++ b/stan/math/prim/prob/multi_normal_lpdf.hpp
@@ -15,7 +15,8 @@
 namespace stan {
 namespace math {
 
-template <bool propto, typename T_y, typename T_loc, typename T_covar>
+template <bool propto, typename T_y, typename T_loc, typename T_covar,
+          require_all_vector_t<T_y, T_loc>* = nullptr>
 return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(const T_y& y,
                                                      const T_loc& mu,
                                                      const T_covar& Sigma) {
@@ -102,6 +103,52 @@ return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(const T_y& y,
   return lp;
 }
 
+template <bool propto, typename T_y, typename T_loc, typename T_covar,
+          require_all_dense_dynamic_t<T_y, T_loc, T_covar>* = nullptr>
+return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(const T_y& y,
+                                                     const T_loc& mu,
+                                                     const T_covar& Sigma) {
+  using T_covar_elem = typename scalar_type<T_covar>::type;
+  static const char* function = "multi_normal_lpdf";
+  check_positive(function, "Covariance matrix rows", Sigma.rows());
+  auto&& Sigma_ref = to_ref(Sigma);
+  const Eigen::Index number_of_y = y.rows();
+  const Eigen::Index number_of_mu = mu.rows();
+  if (unlikely(number_of_y == 0 || number_of_mu == 0)) {
+    return 0.0;
+  }
+  check_size_match(function, "Rows of random variable", y.rows(),
+                  "rows of location parameter", mu.rows());
+  check_size_match(function, "Columns of random variable", y.cols(),
+                  "columns of location parameter", mu.cols());
+  check_size_match(function, "Columns of random variable", y.cols(),
+                  "rows of covariance parameter", Sigma.rows());
+  return_type_t<T_y, T_loc, T_covar> lp(0.0);
+  auto&& y_ref = to_ref(y);
+  auto&& mu_ref = to_ref(mu);
+  const Eigen::Index size_y = y_ref.cols();
+  const Eigen::Index size_mu = mu_ref.cols();
+  check_finite(function, "Location parameter", mu_ref);
+  check_not_nan(function, "Random variable", y_ref);
+  check_symmetric(function, "Covariance matrix", Sigma_ref);
+  auto ldlt_Sigma = make_ldlt_factor(Sigma_ref);
+  check_ldlt_factor(function, "LDLT_Factor of covariance parameter",
+                   ldlt_Sigma);
+
+  if (include_summand<propto>::value) {
+   lp += NEG_LOG_SQRT_TWO_PI * y.rows() * y.cols();
+  }
+
+  if (include_summand<propto, T_covar_elem>::value) {
+   lp -= 0.5 * log_determinant_ldlt(ldlt_Sigma) * y.rows();
+  }
+
+  if (include_summand<propto, T_y, T_loc, T_covar_elem>::value) {
+   lp -= 0.5 * trace_inv_quad_form_ldlt(ldlt_Sigma, (y_ref - mu_ref).transpose());
+  }
+  return lp;
+}
+
 template <typename T_y, typename T_loc, typename T_covar>
 inline return_type_t<T_y, T_loc, T_covar> multi_normal_lpdf(
     const T_y& y, const T_loc& mu, const T_covar& Sigma) {
diff --git a/stan/math/prim/prob/multi_normal_rng.hpp b/stan/math/prim/prob/multi_normal_rng.hpp
index 8f342306164..e1394192306 100644
--- a/stan/math/prim/prob/multi_normal_rng.hpp
+++ b/stan/math/prim/prob/multi_normal_rng.hpp
@@ -29,7 +29,7 @@ namespace math {
  * std::invalid_argument if the length of (each) mu is not equal to the number
  * of rows and columns in S
  */
-template <typename T_loc, class RNG>
+template <typename T_loc, class RNG, require_vector_t<T_loc>* = nullptr>
 inline typename StdVectorBuilder<true, Eigen::VectorXd, T_loc>::type
 multi_normal_rng(const T_loc& mu,
                  const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& S,
@@ -80,6 +80,23 @@ multi_normal_rng(const T_loc& mu,
   return output.data();
 }
 
+template <typename T_loc, class RNG, require_not_vector_t<T_loc>* = nullptr>
+inline typename Eigen::MatrixXd multi_normal_rng(
+    const T_loc& mu,
+    const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& S, RNG& rng) {
+  using Eigen::MatrixXd;
+
+  const size_t N = mu.rows();
+  const size_t M = mu.cols();
+
+  MatrixXd output(N, M);
+  for (int n = 0; n < N; ++n) {
+    auto current_mu = mu.row(n);
+    output.row(n) = multi_normal_rng(current_mu, S, rng);
+  }
+  return output;
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/test/unit/math/mix/prob/multi_normal_test.cpp b/test/unit/math/mix/prob/multi_normal_test.cpp
index bb1fc347d70..ddbd6c37f90 100644
--- a/test/unit/math/mix/prob/multi_normal_test.cpp
+++ b/test/unit/math/mix/prob/multi_normal_test.cpp
@@ -47,6 +47,22 @@ TEST(ProbDistributionsMultiNormal, matvar) {
   stan::test::expect_ad_matvar(f, y1, mu1, Sigma00);
 }
 
+TEST(ProbDistributionsMultiNormal, matrix) {
+  auto f = [](const auto& y, const auto& mu, const auto& sigma) {
+    auto sigma_sym = stan::math::multiply(0.5, sigma + sigma.transpose());
+    return stan::math::multi_normal_lpdf(y, mu, sigma_sym);
+  };
+
+  Eigen::MatrixXd y(2, 3);
+  y << 2.0, -2.0, 11.0, 2.0, -2.0, 11.0;
+  Eigen::MatrixXd mu(2, 3);
+  mu << 1.0, -1.0, 3.0, 1.0, -1.0, 3.0;
+  Eigen::MatrixXd Sigma(3, 3);
+  Sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 0.0, 0.0, 0.0, 5.0;
+  stan::test::expect_ad(f, y, mu, Sigma);
+  stan::test::expect_ad_matvar(f, y, mu, Sigma);
+}
+
 TEST(ProbDistributionsMultiNormal, fvar_var) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
diff --git a/test/unit/math/prim/prob/multi_normal_test.cpp b/test/unit/math/prim/prob/multi_normal_test.cpp
index 4cf0f4bcf3d..7e450d28fde 100644
--- a/test/unit/math/prim/prob/multi_normal_test.cpp
+++ b/test/unit/math/prim/prob/multi_normal_test.cpp
@@ -21,6 +21,29 @@ TEST(ProbDistributionsMultiNormal, NotVectorized) {
   EXPECT_NO_THROW(stan::math::multi_normal_rng(mu, Sigma, rng));
 }
 
+TEST(ProbDistributionsMultiNormal, MatrixNotVectorized) {
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+
+  boost::random::mt19937 rng;
+  Matrix<double, Dynamic, Dynamic> y(2, 3);
+  y << 2.0, -2.0, 11.0, 2.0, -2.0, 11.0;
+  Matrix<double, Dynamic, Dynamic> mu(2, 3);
+  mu << 1.0, -1.0, 3.0, 1.0, -1.0, 3.0;
+  Matrix<double, Dynamic, Dynamic> Sigma(3, 3);
+  Sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 0.0, 0.0, 0.0, 5.0;
+  std::vector<Eigen::VectorXd> test_y;
+  std::vector<Eigen::VectorXd> test_mu;
+  test_y.push_back(y.row(0).transpose());
+  test_y.push_back(y.row(1).transpose());
+  test_mu.push_back(mu.row(0).transpose());
+  test_mu.push_back(mu.row(1).transpose());
+  auto test_eval = stan::math::multi_normal_lpdf(y, mu, Sigma);
+  auto test_orig = stan::math::multi_normal_lpdf(test_y, test_mu, Sigma);
+  EXPECT_FLOAT_EQ(test_orig, test_eval);
+  EXPECT_NO_THROW(stan::math::multi_normal_rng(mu, Sigma, rng));
+}
+
 TEST(ProbDistributionsMultiNormal, Vectorized) {
   using Eigen::Dynamic;
   using Eigen::Matrix;