diff --git a/stan/math/prim/err/check_consistent_sizes.hpp b/stan/math/prim/err/check_consistent_sizes.hpp
index ec1cee072e6..d01636e33b0 100644
--- a/stan/math/prim/err/check_consistent_sizes.hpp
+++ b/stan/math/prim/err/check_consistent_sizes.hpp
@@ -1,8 +1,10 @@
 #ifndef STAN_MATH_PRIM_ERR_CHECK_CONSISTENT_SIZES_HPP
 #define STAN_MATH_PRIM_ERR_CHECK_CONSISTENT_SIZES_HPP
 
+#include <stan/math/prim/err/check_matching_sizes.hpp>
 #include <stan/math/prim/err/invalid_argument.hpp>
 #include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/meta/is_container.hpp>
 #include <stan/math/prim/meta/require_generics.hpp>
 #include <algorithm>
 #include <sstream>
@@ -66,6 +68,34 @@ inline void check_consistent_sizes(const char* function, const char* name1,
   }
 }
 
+/**
+ * Check that the unnamed container inputs have the same size. Inputs are
+ * labeled `arg1`, `arg2`, and so on in error messages.
+ *
+ * @tparam T type of the first input
+ * @tparam Types types of the remaining inputs
+ * @param function function name (for error messages)
+ * @param x first input
+ * @param xs remaining inputs
+ * @throw `invalid_argument` if sizes are inconsistent
+ */
+template <typename T, typename... Types,
+          require_all_container_t<T, Types...>* = nullptr>
+inline void check_consistent_sizes(const char* function, T&& x, Types&&... xs) {
+  std::size_t arg_idx = 2;
+  (
+      [&](const auto& y) {
+        if (x.size() != y.size()) {
+          [&]() STAN_COLD_PATH {
+            const std::string name = "arg" + std::to_string(arg_idx);
+            check_matching_sizes(function, "arg1", x, name.c_str(), y);
+          }();
+        }
+        ++arg_idx;
+      }(xs),
+      ...);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/err/check_matching_dims.hpp b/stan/math/prim/err/check_matching_dims.hpp
index 5308545834f..7856b9a7b97 100644
--- a/stan/math/prim/err/check_matching_dims.hpp
+++ b/stan/math/prim/err/check_matching_dims.hpp
@@ -155,6 +155,35 @@ inline void check_matching_dims(const char* function, const char* name1,
   check_matching_dims(function, name1, y1, name2, y2);
 }
 
+/**
+ * Check that the unnamed Eigen inputs have the same dimensions. Inputs are
+ * labeled `arg1`, `arg2`, and so on in error messages.
+ *
+ * @tparam T type of the first input
+ * @tparam Types types of the remaining inputs
+ * @param function function name (for error messages)
+ * @param x first input
+ * @param xs remaining inputs
+ * @throw `invalid_argument` if dimensions do not match
+ */
+template <typename T, typename... Types,
+          require_all_eigen_t<T, Types...>* = nullptr>
+inline void check_matching_dims(const char* function, const T& x,
+                                const Types&... xs) {
+  std::size_t arg_idx = 2;
+  (
+      [&](const auto& y) {
+        if (x.rows() != y.rows() || x.cols() != y.cols()) {
+          [&]() STAN_COLD_PATH {
+            const std::string name = "arg" + std::to_string(arg_idx);
+            check_matching_dims(function, "arg1", x, name.c_str(), y);
+          }();
+        }
+        ++arg_idx;
+      }(xs),
+      ...);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/functor.hpp b/stan/math/prim/functor.hpp
index 2f9de59c706..22a3934f99d 100644
--- a/stan/math/prim/functor.hpp
+++ b/stan/math/prim/functor.hpp
@@ -23,6 +23,7 @@
 #include <stan/math/prim/functor/ode_ckrk.hpp>
 #include <stan/math/prim/functor/ode_rk45.hpp>
 #include <stan/math/prim/functor/ode_store_sensitivities.hpp>
+#include <stan/math/prim/functor/map.hpp>
 #include <stan/math/prim/functor/map_if.hpp>
 #include <stan/math/prim/functor/map_rect.hpp>
 #include <stan/math/prim/functor/map_rect_combine.hpp>
diff --git a/stan/math/prim/functor/map.hpp b/stan/math/prim/functor/map.hpp
new file mode 100644
index 00000000000..c055693cf2f
--- /dev/null
+++ b/stan/math/prim/functor/map.hpp
@@ -0,0 +1,314 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_MAP_HPP
+#define STAN_MATH_PRIM_FUNCTOR_MAP_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/eval.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err/check_consistent_sizes.hpp>
+#include <stan/math/prim/err/check_matching_dims.hpp>
+#include <stan/math/prim/err/check_size_match.hpp>
+#include <stan/math/prim/functor/apply.hpp>
+
+#include <cstddef>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return a std::vector whose i-th element is `f(x[i], args...)`.
+ * Additional arguments are shared across calls.  If `x` is empty, the result
+ * is empty.  The element type of the result is deduced from `f`.
+ *
+ * @tparam F type of functor
+ * @tparam T type of input std::vector
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each element
+ * @param[in] x input std::vector
+ * @param[in] args shared arguments passed to each call
+ * @return std::vector of results
+ */
+template <typename F, typename T, typename... Args,
+          require_std_vector_t<T>* = nullptr>
+inline auto map(F&& f, T&& x, Args&&... args) {
+  using T_return = std::decay_t<decltype(f(x[0], args...))>;
+  std::vector<T_return> result;
+  result.reserve(x.size());
+
+  for (auto&& xi : x) {
+    result.push_back(f(xi, args...));
+  }
+
+  return result;
+}
+
+/**
+ * Return a std::vector whose i-th element is
+ * `f(get<0>(xs)[i], get<1>(xs)[i], ..., args...)`.
+ * The tuple `xs` must contain at least two equal-sized std::vector inputs.
+ * Additional trailing arguments are shared across calls.  If the inputs are
+ * empty, the result is empty.  The element type of the result is deduced from
+ * `f`.
+ *
+ * @tparam F type of functor
+ * @tparam Tuple type of tuple of std::vector inputs
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each tuple of elements
+ * @param[in] xs tuple of std::vector inputs
+ * @param[in] args shared arguments passed to each call
+ * @return std::vector of results
+ * @throw std::invalid_argument if the inputs have different sizes
+ */
+template <typename F, typename Tuple, typename... Args,
+          require_all_tuple_elements_t<is_std_vector, Tuple>* = nullptr>
+inline auto mapN(F&& f, Tuple&& xs, Args&&... args) {
+  static constexpr const char* function = "mapN";
+
+  return apply(
+      [&f, &args...](auto&&... xs_inner) {
+        check_consistent_sizes(function, xs_inner...);
+
+        const std::size_t n
+            = std::get<0>(std::forward_as_tuple(xs_inner...)).size();
+        using T_return = std::decay_t<decltype(f((xs_inner[0])..., args...))>;
+        std::vector<T_return> result;
+        result.reserve(n);
+
+        for (std::size_t i = 0; i < n; ++i) {
+          result.push_back(f((xs_inner[i])..., args...));
+        }
+        return result;
+      },
+      std::forward<Tuple>(xs));
+}
+
+/**
+ * Return a matrix whose i-th row is `f(m.row(i), args...)`.
+ * Additional arguments are shared across calls.  If `m` has zero rows, or if
+ * the first returned row is empty, returns a 0x0 matrix.  All returned rows
+ * must have the same number of columns.
+ *
+ * @tparam F type of functor
+ * @tparam T type of Eigen matrix
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each row
+ * @param[in] m input matrix
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if returned rows have inconsistent sizes
+ */
+template <typename F, typename T, typename... Args,
+          require_eigen_matrix_dynamic_t<T>* = nullptr>
+inline auto row_map(F&& f, T&& m, Args&&... args) {
+  static constexpr const char* function = "row_map";
+  decltype(auto) m_ref = to_ref(std::forward<T>(m));
+
+  using matrix_t = Eigen::Matrix<
+      scalar_type_t<plain_type_t<decltype(f(m_ref.row(0), args...))>>,
+      Eigen::Dynamic, Eigen::Dynamic>;
+
+  const Eigen::Index n_rows = m_ref.rows();
+  if (n_rows == 0) {
+    return matrix_t(0, 0);
+  }
+
+  auto row_i = eval(f(m_ref.row(0), args...));
+
+  if (row_i.size() == 0) {
+    return matrix_t(0, 0);
+  }
+  matrix_t result(n_rows, row_i.cols());
+  result.row(0) = row_i;
+  for (Eigen::Index i = 1; i < n_rows; ++i) {
+    row_i = f(m_ref.row(i), args...);
+    check_size_match(function, "columns of result", result.cols(),
+                     "columns of returned row", row_i.cols());
+    result.row(i) = row_i;
+  }
+  return result;
+}
+
+/**
+ * Return a matrix whose j-th column is `f(m.col(j), args...)`.
+ * Additional arguments are shared across calls.  If `m` has zero columns, or
+ * if the first returned column is empty, returns a 0x0 matrix.  All returned
+ * columns must have the same number of rows.
+ *
+ * @tparam F type of functor
+ * @tparam T type of Eigen matrix
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each column
+ * @param[in] m input matrix
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if returned columns have inconsistent sizes
+ */
+template <typename F, typename T, typename... Args,
+          require_eigen_matrix_dynamic_t<T>* = nullptr>
+inline auto col_map(F&& f, T&& m, Args&&... args) {
+  static constexpr const char* function = "col_map";
+  decltype(auto) m_ref = to_ref(std::forward<T>(m));
+
+  using matrix_t = Eigen::Matrix<
+      scalar_type_t<plain_type_t<decltype(f(m_ref.col(0), args...))>>,
+      Eigen::Dynamic, Eigen::Dynamic>;
+
+  const Eigen::Index n_cols = m_ref.cols();
+  if (n_cols == 0) {
+    return matrix_t(0, 0);
+  }
+
+  auto col_j = eval(f(m_ref.col(0), args...));
+
+  if (col_j.size() == 0) {
+    return matrix_t(0, 0);
+  }
+  matrix_t result(col_j.rows(), n_cols);
+  result.col(0) = col_j;
+  for (Eigen::Index j = 1; j < n_cols; ++j) {
+    col_j = f(m_ref.col(j), args...);
+    check_size_match(function, "rows of result", result.rows(),
+                     "rows of returned column", col_j.rows());
+    result.col(j) = col_j;
+  }
+  return result;
+}
+
+/**
+ * Return a matrix whose i-th row is
+ * `f(get<0>(ms).row(i), get<1>(ms).row(i), ..., args...)`.
+ * The tuple `ms` must contain at least two Eigen matrices with identical
+ * dimensions.  Additional trailing arguments are shared across calls.  If the
+ * inputs have zero rows, or if the first returned row is empty, returns a 0x0
+ * matrix.  All returned rows must have the same number of columns.
+ *
+ * @tparam F type of functor
+ * @tparam Tuple type of tuple of Eigen matrix inputs
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each tuple of rows
+ * @param[in] ms tuple of Eigen matrix inputs
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if the inputs have different dimensions or
+ * returned rows have inconsistent sizes
+ */
+template <typename F, typename Tuple, typename... Args,
+          require_tuple_t<Tuple>* = nullptr>
+inline auto row_mapN(F&& f, Tuple&& ms, Args&&... args) {
+  static constexpr const char* function = "row_mapN";
+
+  return apply(
+      [&f, &args...](auto&&... mats) {
+        static_assert(
+            (is_eigen_matrix_dynamic<std::decay_t<decltype(mats)>>::value
+             && ...),
+            "row_mapN tuple elements must be Eigen matrices.");
+        check_matching_dims(function, mats...);
+
+        auto m_refs = std::tuple{to_ref(std::forward<decltype(mats)>(mats))...};
+        const Eigen::Index n_rows = std::get<0>(m_refs).rows();
+
+        return apply(
+            [&f, &args..., &n_rows](auto&&... mrs) {
+              using matrix_t
+                  = Eigen::Matrix<scalar_type_t<plain_type_t<decltype(
+                                      f((mrs.row(0))..., args...))>>,
+                                  Eigen::Dynamic, Eigen::Dynamic>;
+
+              if (n_rows == 0) {
+                return matrix_t(0, 0);
+              }
+
+              auto row_i = eval(f((mrs.row(0))..., args...));
+
+              if (row_i.size() == 0) {
+                return matrix_t(0, 0);
+              }
+              matrix_t result(n_rows, row_i.cols());
+              result.row(0) = row_i;
+              for (Eigen::Index i = 1; i < n_rows; ++i) {
+                row_i = f((mrs.row(i))..., args...);
+                check_size_match(function, "columns of result", result.cols(),
+                                 "columns of returned row", row_i.cols());
+                result.row(i) = row_i;
+              }
+              return result;
+            },
+            m_refs);
+      },
+      std::forward<Tuple>(ms));
+}
+
+/**
+ * Return a matrix whose j-th column is
+ * `f(get<0>(ms).col(j), get<1>(ms).col(j), ..., args...)`.
+ * The tuple `ms` must contain at least two Eigen matrices with identical
+ * dimensions.  Additional trailing arguments are shared across calls.  If the
+ * inputs have zero columns, or if the first returned column is empty, returns
+ * a 0x0 matrix.  All returned columns must have the same number of rows.
+ *
+ * @tparam F type of functor
+ * @tparam Tuple type of tuple of Eigen matrix inputs
+ * @tparam Args types of shared arguments
+ * @param[in] f functor applied to each tuple of columns
+ * @param[in] ms tuple of Eigen matrix inputs
+ * @param[in] args shared arguments passed to each call
+ * @return matrix of results
+ * @throw std::invalid_argument if the inputs have different dimensions or
+ * returned columns have inconsistent sizes
+ */
+template <typename F, typename Tuple, typename... Args,
+          require_tuple_t<Tuple>* = nullptr>
+inline auto col_mapN(F&& f, Tuple&& ms, Args&&... args) {
+  static constexpr const char* function = "col_mapN";
+
+  return apply(
+      [&f, &args...](auto&&... mats) {
+        static_assert(
+            (is_eigen_matrix_dynamic<std::decay_t<decltype(mats)>>::value
+             && ...),
+            "col_mapN tuple elements must be Eigen matrices.");
+        check_matching_dims(function, mats...);
+
+        auto m_refs = std::tuple{to_ref(std::forward<decltype(mats)>(mats))...};
+        const Eigen::Index n_cols = std::get<0>(m_refs).cols();
+
+        return apply(
+            [&f, &args..., &n_cols](auto&&... mrs) {
+              using matrix_t
+                  = Eigen::Matrix<scalar_type_t<plain_type_t<decltype(
+                                      f((mrs.col(0))..., args...))>>,
+                                  Eigen::Dynamic, Eigen::Dynamic>;
+
+              if (n_cols == 0) {
+                return matrix_t(0, 0);
+              }
+
+              auto col_j = eval(f((mrs.col(0))..., args...));
+
+              if (col_j.size() == 0) {
+                return matrix_t(0, 0);
+              }
+              matrix_t result(col_j.rows(), n_cols);
+              result.col(0) = col_j;
+              for (Eigen::Index j = 1; j < n_cols; ++j) {
+                col_j = f((mrs.col(j))..., args...);
+                check_size_match(function, "rows of result", result.rows(),
+                                 "rows of returned column", col_j.rows());
+                result.col(j) = col_j;
+              }
+              return result;
+            },
+            m_refs);
+      },
+      std::forward<Tuple>(ms));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/meta/is_tuple.hpp b/stan/math/prim/meta/is_tuple.hpp
index e0c5e2949fc..1ef4cca1e22 100644
--- a/stan/math/prim/meta/is_tuple.hpp
+++ b/stan/math/prim/meta/is_tuple.hpp
@@ -3,6 +3,7 @@
 
 #include <cstddef>
 #include <stan/math/prim/meta/require_helpers.hpp>
+#include <stan/math/prim/meta/conjunction.hpp>
 #include <tuple>
 #include <type_traits>
 
@@ -15,6 +16,13 @@ struct is_tuple_impl : std::false_type {};
 
 template <typename... Types>
 struct is_tuple_impl<std::tuple<Types...>> : std::true_type {};
+
+template <template <class...> class TypeCheck, typename T>
+struct all_tuple_elements : std::false_type {};
+
+template <template <class...> class TypeCheck, typename... Types>
+struct all_tuple_elements<TypeCheck, std::tuple<Types...>>
+    : math::conjunction<TypeCheck<std::decay_t<Types>>...> {};
 }  // namespace internal
 
 template <typename T>
@@ -32,6 +40,13 @@ constexpr bool is_tuple_v = is_tuple<T>::value;
 template <typename T>
 using require_tuple_t = require_t<is_tuple<std::decay_t<T>>>;
 
+/*! \brief Require all elements of a tuple satisfy `TypeCheck` */
+/*! @tparam TypeCheck the type trait applied to each tuple element */
+/*! @tparam T the tuple type to check */
+template <template <class...> class TypeCheck, typename T>
+using require_all_tuple_elements_t
+    = require_t<internal::all_tuple_elements<TypeCheck, std::decay_t<T>>>;
+
 /*! \brief Require type does not satisfy @ref is_tuple */
 /*! @tparam T the type to check */
 template <typename T>
diff --git a/test/unit/math/prim/err/check_consistent_sizes_test.cpp b/test/unit/math/prim/err/check_consistent_sizes_test.cpp
new file mode 100644
index 00000000000..62fcc9eaa46
--- /dev/null
+++ b/test/unit/math/prim/err/check_consistent_sizes_test.cpp
@@ -0,0 +1,32 @@
+#include <stan/math/prim/err/check_consistent_sizes.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+#include <vector>
+
+TEST(ErrorHandling, checkConsistentSizesUnnamedArguments) {
+  using stan::math::check_consistent_sizes;
+
+  std::vector<double> x1(3);
+  std::vector<int> x2(3);
+  std::vector<double> x3(2);
+
+  EXPECT_NO_THROW(check_consistent_sizes("check_consistent_sizes", x1, x2));
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_consistent_sizes("check_consistent_sizes", x1, x2, x3),
+      std::invalid_argument, "arg1", 1);
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_consistent_sizes("check_consistent_sizes", x1, x2, x3),
+      std::invalid_argument, "arg3", 1);
+  EXPECT_NO_THROW(check_consistent_sizes("check_consistent_sizes", x1));
+}
+
+TEST(ErrorHandling, checkConsistentSizesNamedArgumentsStillResolve) {
+  using stan::math::check_consistent_sizes;
+
+  std::vector<double> x1(3);
+  std::vector<double> x2(2);
+
+  EXPECT_THROW_MSG(
+      check_consistent_sizes("check_consistent_sizes", "x1", x1, "x2", x2),
+      std::invalid_argument, "x2");
+}
diff --git a/test/unit/math/prim/err/check_matching_dims_test.cpp b/test/unit/math/prim/err/check_matching_dims_test.cpp
index 4643355b52f..f0b62a70493 100644
--- a/test/unit/math/prim/err/check_matching_dims_test.cpp
+++ b/test/unit/math/prim/err/check_matching_dims_test.cpp
@@ -27,6 +27,23 @@ TEST(ErrorHandlingMatrix, checkMatchingDimsMatrix) {
       std::invalid_argument);
 }
 
+TEST(ErrorHandlingMatrix, checkMatchingDimsUnnamedArguments) {
+  using stan::math::check_matching_dims;
+
+  Eigen::MatrixXd x1(2, 3);
+  Eigen::MatrixXd x2(2, 3);
+  Eigen::MatrixXd x3(3, 2);
+
+  EXPECT_NO_THROW(check_matching_dims("checkMatchingDims", x1));
+  EXPECT_NO_THROW(check_matching_dims("checkMatchingDims", x1, x2));
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_matching_dims("checkMatchingDims", x1, x2, x3),
+      std::invalid_argument, "arg1", 1);
+  EXPECT_THROW_MSG_WITH_COUNT(
+      check_matching_dims("checkMatchingDims", x1, x2, x3),
+      std::invalid_argument, "arg3", 1);
+}
+
 TEST(ErrorHandlingMatrix, checkMatchingDimsArray) {
   std::vector<double> y(3);
   std::vector<double> x(3);
diff --git a/test/unit/math/prim/functor/map_test.cpp b/test/unit/math/prim/functor/map_test.cpp
new file mode 100644
index 00000000000..670c891812a
--- /dev/null
+++ b/test/unit/math/prim/functor/map_test.cpp
@@ -0,0 +1,536 @@
+#include <stan/math/prim.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <memory>
+#include <stdexcept>
+#include <tuple>
+#include <vector>
+
+TEST(MathFunctor, map_square) {
+  std::vector<double> x{1.0, 2.0, 3.0};
+  auto y = stan::math::map([](double v) { return v * v; }, x);
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 1.0);
+  EXPECT_FLOAT_EQ(y[1], 4.0);
+  EXPECT_FLOAT_EQ(y[2], 9.0);
+}
+
+TEST(MathFunctor, map_type_change) {
+  std::vector<int> x{1, 2, 3};
+  auto y = stan::math::map([](int v) { return v + 0.5; }, x);
+  EXPECT_FLOAT_EQ(y[0], 1.5);
+  EXPECT_FLOAT_EQ(y[1], 2.5);
+  EXPECT_FLOAT_EQ(y[2], 3.5);
+}
+
+TEST(MathFunctor, map_empty) {
+  std::vector<double> x;
+  auto y = stan::math::map([](double v) { return v; }, x);
+  EXPECT_EQ(y.size(), 0u);
+}
+
+TEST(MathFunctor, map_const_vector) {
+  const std::vector<double> x{1.0, 2.0, 3.0};
+  auto y = stan::math::map([](double v) { return 2.0 * v; }, x);
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 2.0);
+  EXPECT_FLOAT_EQ(y[1], 4.0);
+  EXPECT_FLOAT_EQ(y[2], 6.0);
+}
+
+TEST(MathFunctor, map_variadic_args) {
+  std::vector<double> x{1.0, 2.0, 3.0};
+  double offset = 10.0;
+  double scale = 2.0;
+  auto y = stan::math::map(
+      [](double v, double mu, double sigma) { return mu + sigma * v; }, x,
+      offset, scale);
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 12.0);
+  EXPECT_FLOAT_EQ(y[1], 14.0);
+  EXPECT_FLOAT_EQ(y[2], 16.0);
+}
+
+TEST(MathFunctor, map_shared_arg_reused) {
+  std::vector<double> x{1.0, 2.0, 3.0};
+  auto offset = std::make_unique<double>(10.0);
+  auto y = stan::math::map(
+      [](double v, const std::unique_ptr<double>& mu) { return v + *mu; }, x,
+      std::move(offset));
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 11.0);
+  EXPECT_FLOAT_EQ(y[1], 12.0);
+  EXPECT_FLOAT_EQ(y[2], 13.0);
+}
+
+TEST(MathFunctor, mapN_add) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0, 20.0};
+  auto y = stan::math::mapN([](double u, double v) { return u + v; },
+                            std::forward_as_tuple(a, b));
+  ASSERT_EQ(y.size(), 2u);
+  EXPECT_FLOAT_EQ(y[0], 11.0);
+  EXPECT_FLOAT_EQ(y[1], 22.0);
+}
+
+TEST(MathFunctor, mapN_mixed_scalar_types) {
+  std::vector<int> a{1, 2, 3};
+  std::vector<double> b{0.5, 1.5, 2.5};
+  auto y = stan::math::mapN([](int u, double v) { return u + v; },
+                            std::forward_as_tuple(a, b));
+  ASSERT_EQ(y.size(), 3u);
+  EXPECT_FLOAT_EQ(y[0], 1.5);
+  EXPECT_FLOAT_EQ(y[1], 3.5);
+  EXPECT_FLOAT_EQ(y[2], 5.5);
+}
+
+TEST(MathFunctor, mapN_three_vectors) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0, 20.0};
+  std::vector<double> c{100.0, 200.0};
+  auto y
+      = stan::math::mapN([](double u, double v, double w) { return u + v + w; },
+                         std::forward_as_tuple(a, b, c));
+  ASSERT_EQ(y.size(), 2u);
+  EXPECT_FLOAT_EQ(y[0], 111.0);
+  EXPECT_FLOAT_EQ(y[1], 222.0);
+}
+
+TEST(MathFunctor, mapN_shared_arg_reused) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0, 20.0};
+  auto offset = std::make_unique<double>(100.0);
+  auto y = stan::math::mapN(
+      [](double u, double v, const std::unique_ptr<double>& mu) {
+        return u + v + *mu;
+      },
+      std::forward_as_tuple(a, b), std::move(offset));
+  EXPECT_STD_VECTOR_FLOAT_EQ(y, std::vector<double>({111.0, 122.0}));
+}
+
+TEST(MathFunctor, mapN_empty) {
+  std::vector<double> a;
+  std::vector<double> b;
+  auto y = stan::math::mapN([](double u, double v) { return u + v; },
+                            std::forward_as_tuple(a, b));
+  EXPECT_EQ(y.size(), 0u);
+}
+
+TEST(MathFunctor, mapN_size_mismatch) {
+  std::vector<double> a{1.0, 2.0};
+  std::vector<double> b{10.0};
+  EXPECT_THROW(stan::math::mapN([](double u, double v) { return u + v; },
+                                std::forward_as_tuple(a, b)),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, row_map_standardize) {
+  Eigen::MatrixXd x(2, 3);
+  x << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0;
+  Eigen::RowVectorXd mu(3);
+  mu << 1.0, 2.0, 3.0;
+  Eigen::RowVectorXd sigma(3);
+  sigma << 1.0, 2.0, 3.0;
+
+  auto y = stan::math::row_map(
+      [](const Eigen::RowVectorXd& row, const Eigen::RowVectorXd& m,
+         const Eigen::RowVectorXd& s) {
+        return (row.array() - m.array()) / s.array();
+      },
+      x, mu, sigma);
+
+  Eigen::MatrixXd expected(2, 3);
+  expected << 0.0, 0.0, 0.0, 3.0, 1.5, 1.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_map_empty) {
+  Eigen::MatrixXd x(0, 3);
+  auto y = stan::math::row_map(
+      [](const Eigen::RowVectorXd& row) { return row; }, x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, row_map_empty_result) {
+  Eigen::MatrixXd x = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::row_map(
+      [&calls](const Eigen::RowVectorXd&) {
+        ++calls;
+        return Eigen::RowVectorXd(0);
+      },
+      x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, row_map_inconsistent_row_size) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  int call = 0;
+  EXPECT_THROW(stan::math::row_map(
+                   [&call](const Eigen::RowVectorXd&) {
+                     ++call;
+                     if (call == 1) {
+                       return Eigen::RowVectorXd::Ones(2);
+                     }
+                     return Eigen::RowVectorXd::Ones(3);
+                   },
+                   x),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, row_map_shared_arg_reused) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  auto offset = std::make_unique<double>(10.0);
+  auto y = stan::math::row_map(
+      [](const Eigen::RowVectorXd& row, const std::unique_ptr<double>& mu) {
+        return row.array() + *mu;
+      },
+      x, std::move(offset));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 12.0, 13.0, 14.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_map_sum_rows) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+
+  auto y = stan::math::col_map(
+      [](const Eigen::VectorXd& col) {
+        return Eigen::VectorXd::Constant(1, col.sum());
+      },
+      x);
+
+  Eigen::MatrixXd expected(1, 2);
+  expected << 4.0, 6.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_map_empty) {
+  Eigen::MatrixXd x(3, 0);
+  auto y
+      = stan::math::col_map([](const Eigen::VectorXd& col) { return col; }, x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, col_map_empty_result) {
+  Eigen::MatrixXd x = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::col_map(
+      [&calls](const Eigen::VectorXd&) {
+        ++calls;
+        return Eigen::VectorXd(0);
+      },
+      x);
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, col_map_inconsistent_col_size) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  int call = 0;
+  EXPECT_THROW(stan::math::col_map(
+                   [&call](const Eigen::VectorXd&) {
+                     ++call;
+                     if (call == 1) {
+                       return Eigen::VectorXd::Ones(2);
+                     }
+                     return Eigen::VectorXd::Ones(3);
+                   },
+                   x),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, col_map_shared_arg_reused) {
+  Eigen::MatrixXd x(2, 2);
+  x << 1.0, 2.0, 3.0, 4.0;
+  auto scale = std::make_unique<double>(2.0);
+  auto y = stan::math::col_map(
+      [](const Eigen::VectorXd& col, const std::unique_ptr<double>& s) {
+        return (*s) * col;
+      },
+      x, std::move(scale));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 2.0, 4.0, 6.0, 8.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_add) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+
+  auto y
+      = stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                                const Eigen::RowVectorXd& v) { return u + v; },
+                             std::forward_as_tuple(a, b));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_shared_arg_reused) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  auto offset = std::make_unique<double>(100.0);
+  auto y = stan::math::row_mapN(
+      [](const Eigen::RowVectorXd& u, const Eigen::RowVectorXd& v,
+         const std::unique_ptr<double>& mu) {
+        return *mu + u.array() + v.array();
+      },
+      std::forward_as_tuple(a, b), std::move(offset));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 111.0, 122.0, 133.0, 144.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_dim_mismatch) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 3);
+  b << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0;
+
+  EXPECT_THROW(
+      stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                              const Eigen::RowVectorXd& v) { return u + v; },
+                           std::forward_as_tuple(a, b)),
+      std::invalid_argument);
+}
+
+TEST(MathFunctor, row_mapN_empty) {
+  Eigen::MatrixXd a(0, 2);
+  Eigen::MatrixXd b(0, 2);
+  auto y
+      = stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                                const Eigen::RowVectorXd& v) { return u + v; },
+                             std::forward_as_tuple(a, b));
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, row_mapN_empty_result) {
+  Eigen::MatrixXd a = Eigen::MatrixXd::Ones(2, 2);
+  Eigen::MatrixXd b = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::row_mapN(
+      [&calls](const Eigen::RowVectorXd&, const Eigen::RowVectorXd&) {
+        ++calls;
+        return Eigen::RowVectorXd(0);
+      },
+      std::forward_as_tuple(a, b));
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, row_mapN_mixed_scalar_types) {
+  Eigen::MatrixXi a(2, 2);
+  a << 1, 2, 3, 4;
+  Eigen::MatrixXd b(2, 2);
+  b << 0.5, 1.5, 2.5, 3.5;
+  auto y = stan::math::row_mapN(
+      [](const auto& u, const auto& v) {
+        return u.template cast<double>() + v;
+      },
+      std::forward_as_tuple(a, b));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 1.5, 3.5, 5.5, 7.5;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_inconsistent_row_size) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  int call = 0;
+  EXPECT_THROW(
+      stan::math::row_mapN(
+          [&call](const Eigen::RowVectorXd&, const Eigen::RowVectorXd&) {
+            ++call;
+            if (call == 1) {
+              return Eigen::RowVectorXd::Ones(2);
+            }
+            return Eigen::RowVectorXd::Ones(3);
+          },
+          std::forward_as_tuple(a, b)),
+      std::invalid_argument);
+}
+
+TEST(MathFunctor, col_mapN_add) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v) { return u + v; },
+      std::forward_as_tuple(a, b));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_shared_arg_reused) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  auto offset = std::make_unique<double>(100.0);
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v,
+         const std::unique_ptr<double>& mu) {
+        return *mu + u.array() + v.array();
+      },
+      std::forward_as_tuple(a, b), std::move(offset));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 111.0, 122.0, 133.0, 144.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_dim_mismatch) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(3, 2);
+  b << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0;
+
+  EXPECT_THROW(
+      stan::math::col_mapN([](const Eigen::VectorXd& u,
+                              const Eigen::VectorXd& v) { return u + v; },
+                           std::forward_as_tuple(a, b)),
+      std::invalid_argument);
+}
+
+TEST(MathFunctor, col_mapN_empty) {
+  Eigen::MatrixXd a(2, 0);
+  Eigen::MatrixXd b(2, 0);
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v) { return u + v; },
+      std::forward_as_tuple(a, b));
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+}
+
+TEST(MathFunctor, col_mapN_empty_result) {
+  Eigen::MatrixXd a = Eigen::MatrixXd::Ones(2, 2);
+  Eigen::MatrixXd b = Eigen::MatrixXd::Ones(2, 2);
+  int calls = 0;
+  auto y = stan::math::col_mapN(
+      [&calls](const Eigen::VectorXd&, const Eigen::VectorXd&) {
+        ++calls;
+        return Eigen::VectorXd(0);
+      },
+      std::forward_as_tuple(a, b));
+  EXPECT_EQ(y.rows(), 0);
+  EXPECT_EQ(y.cols(), 0);
+  EXPECT_EQ(calls, 1);
+}
+
+TEST(MathFunctor, col_mapN_mixed_scalar_types) {
+  Eigen::MatrixXi a(2, 2);
+  a << 1, 2, 3, 4;
+  Eigen::MatrixXd b(2, 2);
+  b << 0.5, 1.5, 2.5, 3.5;
+  auto y = stan::math::col_mapN(
+      [](const auto& u, const auto& v) {
+        return u.template cast<double>() + v;
+      },
+      std::forward_as_tuple(a, b));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 1.5, 3.5, 5.5, 7.5;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_inconsistent_col_size) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 10.0, 20.0, 30.0, 40.0;
+  int call = 0;
+  EXPECT_THROW(stan::math::col_mapN(
+                   [&call](const Eigen::VectorXd&, const Eigen::VectorXd&) {
+                     ++call;
+                     if (call == 1) {
+                       return Eigen::VectorXd::Ones(2);
+                     }
+                     return Eigen::VectorXd::Ones(3);
+                   },
+                   std::forward_as_tuple(a, b)),
+               std::invalid_argument);
+}
+
+TEST(MathFunctor, row_map_block_expression) {
+  Eigen::MatrixXd x(3, 3);
+  x << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0;
+  auto y = stan::math::row_map([](const auto& row) { return 2.0 * row; },
+                               x.block(0, 0, 2, 3));
+
+  Eigen::MatrixXd expected(2, 3);
+  expected << 2.0, 4.0, 6.0, 8.0, 10.0, 12.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_map_block_expression) {
+  Eigen::MatrixXd x(3, 3);
+  x << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0;
+  auto y = stan::math::col_map([](const auto& col) { return 2.0 * col; },
+                               x.block(0, 0, 3, 2));
+
+  Eigen::MatrixXd expected(3, 2);
+  expected << 2.0, 4.0, 8.0, 10.0, 14.0, 16.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, row_mapN_expression_inputs) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 0.0, 0.0, 1.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd c(2, 2);
+  c << 10.0, 20.0, 30.0, 40.0;
+
+  auto y
+      = stan::math::row_mapN([](const Eigen::RowVectorXd& u,
+                                const Eigen::RowVectorXd& v) { return u + v; },
+                             std::forward_as_tuple(a * b, c));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
+
+TEST(MathFunctor, col_mapN_expression_inputs) {
+  Eigen::MatrixXd a(2, 2);
+  a << 1.0, 0.0, 0.0, 1.0;
+  Eigen::MatrixXd b(2, 2);
+  b << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd c(2, 2);
+  c << 10.0, 20.0, 30.0, 40.0;
+
+  auto y = stan::math::col_mapN(
+      [](const Eigen::VectorXd& u, const Eigen::VectorXd& v) { return u + v; },
+      std::forward_as_tuple(a * b, c));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 11.0, 22.0, 33.0, 44.0;
+  EXPECT_MATRIX_FLOAT_EQ(y, expected);
+}
diff --git a/test/unit/math/prim/meta/is_tuple_test.cpp b/test/unit/math/prim/meta/is_tuple_test.cpp
index dfc85e51a1e..b21e22f21d4 100644
--- a/test/unit/math/prim/meta/is_tuple_test.cpp
+++ b/test/unit/math/prim/meta/is_tuple_test.cpp
@@ -10,3 +10,19 @@ TEST(MathMetaPrim, is_tuple) {
   EXPECT_FALSE((stan::math::is_tuple<double>::value));
   EXPECT_FALSE((stan::math::is_tuple<std::vector<double>>::value));
 }
+
+TEST(MathMetaPrim, all_tuple_elements) {
+  using stan::is_std_vector;
+  using stan::math::internal::all_tuple_elements;
+  EXPECT_TRUE((all_tuple_elements<
+               is_std_vector,
+               std::tuple<std::vector<int>, std::vector<double>>>::value));
+  EXPECT_FALSE(
+      (all_tuple_elements<is_std_vector,
+                          std::tuple<std::vector<int>, double>>::value));
+  EXPECT_TRUE((all_tuple_elements<
+               is_std_vector,
+               std::tuple<std::vector<double>&, std::vector<double>&>>::value));
+  EXPECT_FALSE((all_tuple_elements<is_std_vector, std::vector<double>>::value));
+  EXPECT_FALSE((all_tuple_elements<is_std_vector, double>::value));
+}