diff --git a/stan/math/opencl/prim/eigenvectors_sym.hpp b/stan/math/opencl/prim/eigenvectors_sym.hpp
index bf61d8b3b68..73d65f0c27f 100644
--- a/stan/math/opencl/prim/eigenvectors_sym.hpp
+++ b/stan/math/opencl/prim/eigenvectors_sym.hpp
@@ -10,7 +10,7 @@ namespace math {
 
 inline matrix_cl<double> eigenvectors_sym(const matrix_cl<double>& m) {
   check_nonzero_size("eigenvectors_sym", "m", m);
-  check_symmetric("eigenvalues_sym", "m", m);
+  check_symmetric("eigenvectors_sym", "m", m);
 
   matrix_cl<double> eigenvalues, eigenvectors;
   symmetric_eigensolver(m, eigenvalues, eigenvectors);
diff --git a/stan/math/prim/err/check_symmetric.hpp b/stan/math/prim/err/check_symmetric.hpp
index 58b3f9f4e79..088023593ff 100644
--- a/stan/math/prim/err/check_symmetric.hpp
+++ b/stan/math/prim/err/check_symmetric.hpp
@@ -7,7 +7,7 @@
 #include <stan/math/prim/err/constraint_tolerance.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/abs.hpp>
-#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <sstream>
 #include <string>
@@ -32,9 +32,6 @@ template <typename EigMat, require_matrix_t<EigMat>* = nullptr>
 inline void check_symmetric(const char* function, const char* name,
                             const EigMat& y) {
   check_square(function, name, y);
-  using stan::math::fabs;
-  using std::fabs;
-
   Eigen::Index k = y.rows();
   if (k <= 1) {
     return;
@@ -42,7 +39,8 @@ inline void check_symmetric(const char* function, const char* name,
   const auto& y_ref = to_ref(y);
   for (Eigen::Index m = 0; m < k; ++m) {
     for (Eigen::Index n = m + 1; n < k; ++n) {
-      if (!(fabs(value_of(y_ref(m, n)) - value_of(y_ref(n, m)))
+      if (!(stan::math::abs(value_of_rec(y_ref(m, n))
+                            - value_of_rec(y_ref(n, m)))
             <= CONSTRAINT_TOLERANCE)) {
         [&]() STAN_COLD_PATH {
           std::ostringstream msg1;
diff --git a/stan/math/prim/fun/eigenvectors_sym.hpp b/stan/math/prim/fun/eigenvectors_sym.hpp
index 1aa9718f1e3..109009ca7e2 100644
--- a/stan/math/prim/fun/eigenvectors_sym.hpp
+++ b/stan/math/prim/fun/eigenvectors_sym.hpp
@@ -17,7 +17,7 @@ eigenvectors_sym(EigMat&& m) {
   }
   using PlainMat = plain_type_t<EigMat>;
   decltype(auto) m_ref = to_ref(std::forward<EigMat>(m));
-  check_symmetric("eigenvalues_sym", "m", m_ref);
+  check_symmetric("eigenvectors_sym", "m", m_ref);
 
   Eigen::SelfAdjointEigenSolver<PlainMat> solver(
       std::forward<decltype(m_ref)>(m_ref));
diff --git a/test/unit/math/mix/fun/eigendecompose_sym_test.cpp b/test/unit/math/mix/fun/eigendecompose_sym_test.cpp
index e4aba4f615b..d9b78d1a001 100644
--- a/test/unit/math/mix/fun/eigendecompose_sym_test.cpp
+++ b/test/unit/math/mix/fun/eigendecompose_sym_test.cpp
@@ -104,3 +104,30 @@ TEST(MathMixMatFun, eigendecomposeSym_varmat) {
   vec_tols.hessian_hessian_ = 1e-2;
   stan::test::expect_ad_matvar(vec_tols, g, a33);
 }
+
+TEST(MathMixMatFun, eigendecomposeSymComplex) {
+  // complex input with real entries, the only case where symmetric and
+  // Hermitian coincide; the imaginary parts are constant
+  auto g = [](const auto& y) {
+    auto a = ((y + y.transpose()) * 0.5).eval();
+    return std::get<0>(
+        stan::math::eigendecompose_sym(stan::math::to_complex(a, 0)));
+  };
+  auto f = [](const auto& y) {
+    auto a = ((y + y.transpose()) * 0.5).eval();
+    return std::get<1>(
+        stan::math::eigendecompose_sym(stan::math::to_complex(a, 0)));
+  };
+
+  stan::test::ad_tolerances values_tols;
+  values_tols.hessian_hessian_ = 5e-2;
+  values_tols.hessian_fvar_hessian_ = 5e-2;
+  stan::test::ad_tolerances vec_tols;
+  vec_tols.hessian_fvar_hessian_ = 1e-2;
+  vec_tols.hessian_hessian_ = 1e-2;
+
+  Eigen::MatrixXd a33(3, 3);
+  a33 << 1, 2, 3, 2, 5, 7.9, 3, 7.9, 1.08;
+  stan::test::expect_ad(values_tols, f, a33);
+  stan::test::expect_ad(vec_tols, g, a33);
+}
diff --git a/test/unit/math/mix/fun/eigenvalues_sym_test.cpp b/test/unit/math/mix/fun/eigenvalues_sym_test.cpp
index 4b8c3cde5be..cae15ae569d 100644
--- a/test/unit/math/mix/fun/eigenvalues_sym_test.cpp
+++ b/test/unit/math/mix/fun/eigenvalues_sym_test.cpp
@@ -69,3 +69,20 @@ TEST(MathMixMatFun, eigenvaluesSym_varmat) {
   stan::test::expect_ad(tols, f, a22);
   stan::test::expect_ad_matvar(tols, f, a22);
 }
+
+TEST(MathMixMatFun, eigenvaluesSymComplex) {
+  // complex input with real entries, the only case where symmetric and
+  // Hermitian coincide; the imaginary parts are constant
+  auto f = [](const auto& y) {
+    auto a = ((y + y.transpose()) * 0.5).eval();
+    return stan::math::eigenvalues_sym(stan::math::to_complex(a, 0));
+  };
+
+  stan::test::ad_tolerances tols;
+  tols.hessian_hessian_ = 5e-2;
+  tols.hessian_fvar_hessian_ = 5e-2;
+
+  Eigen::MatrixXd a33(3, 3);
+  a33 << 1, 2, 3, 2, 5, 7.9, 3, 7.9, 1.08;
+  stan::test::expect_ad(tols, f, a33);
+}
diff --git a/test/unit/math/mix/fun/eigenvectors_sym_test.cpp b/test/unit/math/mix/fun/eigenvectors_sym_test.cpp
index 199e34602da..b9d7ca709f8 100644
--- a/test/unit/math/mix/fun/eigenvectors_sym_test.cpp
+++ b/test/unit/math/mix/fun/eigenvectors_sym_test.cpp
@@ -67,3 +67,20 @@ TEST(MathMixMatFun, eigenvectorsSym_varmat) {
   stan::test::expect_ad(tols, f, a33);
   stan::test::expect_ad_matvar(tols, f, a33);
 }
+
+TEST(MathMixMatFun, eigenvectorsSymComplex) {
+  // complex input with real entries, the only case where symmetric and
+  // Hermitian coincide; the imaginary parts are constant
+  auto f = [](const auto& y) {
+    auto a = ((y + y.transpose()) * 0.5).eval();
+    return stan::math::eigenvectors_sym(stan::math::to_complex(a, 0));
+  };
+
+  stan::test::ad_tolerances tols;
+  tols.hessian_fvar_hessian_ = 1e-2;
+  tols.hessian_hessian_ = 1e-2;
+
+  Eigen::MatrixXd a33(3, 3);
+  a33 << 1, 2, 3, 2, 5, 7.9, 3, 7.9, 1.08;
+  stan::test::expect_ad(tols, f, a33);
+}