diff --git a/Jenkinsfile b/Jenkinsfile
index cb3217f4025..b039a14c314 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -150,7 +150,7 @@ pipeline {
         stage('Verify changes') {
             agent { label 'linux' }
             steps {
-                script {         
+                script {
 
                     retry(3) { checkout scm }
                     sh 'git clean -xffd'
@@ -181,7 +181,7 @@ pipeline {
                 deleteDir()
                 unstash 'MathSetup'
                 sh "echo CXX=${MPICXX} >> make/local"
-                sh "echo CXX_TYPE=gcc >> make/local"                        
+                sh "echo CXX_TYPE=gcc >> make/local"
                 sh "echo STAN_MPI=true >> make/local"
                 runTests("test/unit")
             }
@@ -242,6 +242,7 @@ pipeline {
                         unstash 'MathSetup'
                         sh "echo CXX=${env.CXX} -Werror > make/local"
                         sh "echo CPPFLAGS+=-DSTAN_THREADS >> make/local"
+                        sh "export STAN_NUM_THREADS=4"
                         runTests("test/unit -f thread")
                         sh "find . -name *_test.xml | xargs rm"
                         runTests("test/unit -f map_rect")
@@ -272,7 +273,7 @@ pipeline {
             }
         }
         stage('Additional merge tests') {
-            when { 
+            when {
                 allOf {
                     anyOf {
                         branch 'develop'
@@ -309,12 +310,12 @@ pipeline {
             }
         }
         stage('Upstream tests') {
-            when { 
+            when {
                 allOf {
-                    expression { 
-                        env.BRANCH_NAME ==~ /PR-\d+/ 
+                    expression {
+                        env.BRANCH_NAME ==~ /PR-\d+/
                     }
-                    expression { 
+                    expression {
                         !skipRemainingStages
                     }
                 }
diff --git a/stan/math/fwd/functor.hpp b/stan/math/fwd/functor.hpp
index c65b467edf4..cb6b200f3a6 100644
--- a/stan/math/fwd/functor.hpp
+++ b/stan/math/fwd/functor.hpp
@@ -4,5 +4,6 @@
 #include <stan/math/fwd/functor/gradient.hpp>
 #include <stan/math/fwd/functor/hessian.hpp>
 #include <stan/math/fwd/functor/jacobian.hpp>
+#include <stan/math/fwd/functor/reduce_sum.hpp>
 
 #endif
diff --git a/stan/math/fwd/functor/reduce_sum.hpp b/stan/math/fwd/functor/reduce_sum.hpp
new file mode 100644
index 00000000000..4ebf14fa245
--- /dev/null
+++ b/stan/math/fwd/functor/reduce_sum.hpp
@@ -0,0 +1,104 @@
+#ifndef STAN_MATH_FWD_FUNCTOR_REDUCE_SUM_HPP
+#define STAN_MATH_FWD_FUNCTOR_REDUCE_SUM_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+
+#include <tbb/task_arena.h>
+#include <tbb/parallel_reduce.h>
+#include <tbb/blocked_range.h>
+
+#include <algorithm>
+#include <tuple>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+template <typename ReduceFunction, typename Enable, typename ReturnType,
+          typename Vec, typename... Args>
+struct reduce_sum_impl;
+
+/**
+ * reduce_sum_impl implementation for any autodiff type.
+ *
+ * @tparam ReduceFunction Type of reducer function
+ * @tparam ReturnType An arithmetic type
+ * @tparam Vec Type of sliced argument
+ * @tparam Args Types of shared arguments
+ */
+template <typename ReduceFunction, typename Enable, typename ReturnType,
+          typename Vec, typename... Args>
+struct reduce_sum_impl {
+  /**
+   * Call an instance of the function `ReduceFunction` on every element
+   *   of an input sequence and sum these terms.
+   *
+   * This specialization is not parallelized and works for any autodiff types.
+   *
+   * An instance, f, of `ReduceFunction` should have the signature:
+   *   T f(int start, int end, Vec&& vmapped_subset, std::ostream* msgs,
+   * Args&&... args)
+   *
+   * `ReduceFunction` must be default constructible without any arguments
+   *
+   * Each call to `ReduceFunction` is responsible for computing the
+   *   start through end terms (inclusive) of the overall sum. All args are
+   * passed from this function through to the `ReduceFunction` instances.
+   *   However, only the start through end  (inclusive) elements of the vmapped
+   * argument are passed to the `ReduceFunction` instances (as the
+   * `vmapped_subset` argument).
+   *
+   * If auto partitioning is true, do the calculation with one
+   *  ReduceFunction call. If false, break work into pieces strictly smaller
+   *  than grainsize.
+   *
+   * grainsize must be greater than or equal to 1
+   *
+   * @param vmapped Vector containing one element per term of sum
+   * @param auto_partitioning Work partitioning style (ignored)
+   * @param grainsize Suggested grainsize for tbb
+   * @param[in, out] msgs The print stream for warning messages
+   * @param args Shared arguments used in every sum term
+   * @return Summation of all terms
+   */
+  inline return_type_t<Vec, Args...> operator()(Vec&& vmapped,
+                                                bool auto_partitioning,
+                                                int grainsize,
+                                                std::ostream* msgs,
+                                                Args&&... args) const {
+    if (vmapped.empty()) {
+      return 0.0;
+    }
+
+    if (auto_partitioning) {
+      return ReduceFunction()(0, vmapped.size() - 1, std::forward<Vec>(vmapped),
+                              msgs, std::forward<Args>(args)...);
+    } else {
+      return_type_t<Vec, Args...> sum = 0.0;
+      for (size_t i = 0; i < (vmapped.size() + grainsize - 1) / grainsize;
+           ++i) {
+        size_t start = i * grainsize;
+        size_t end = std::min((i + 1) * grainsize, vmapped.size()) - 1;
+
+        std::decay_t<Vec> sub_slice;
+        sub_slice.reserve(end - start + 1);
+        for (size_t i = start; i <= end; ++i) {
+          sub_slice.emplace_back(vmapped[i]);
+        }
+
+        sum += ReduceFunction()(start, end, std::forward<Vec>(sub_slice), msgs,
+                                std::forward<Args>(args)...);
+      }
+      return sum;
+    }
+  }
+};
+
+}  // namespace internal
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/functor.hpp b/stan/math/prim/functor.hpp
index bfa72de41af..f33209dcac3 100644
--- a/stan/math/prim/functor.hpp
+++ b/stan/math/prim/functor.hpp
@@ -18,5 +18,6 @@
 #include <stan/math/prim/functor/mpi_cluster.hpp>
 #include <stan/math/prim/functor/mpi_command.hpp>
 #include <stan/math/prim/functor/mpi_distributed_apply.hpp>
-
+#include <stan/math/prim/functor/reduce_sum.hpp>
+#include <stan/math/prim/functor/reduce_sum_static.hpp>
 #endif
diff --git a/stan/math/prim/functor/reduce_sum.hpp b/stan/math/prim/functor/reduce_sum.hpp
new file mode 100644
index 00000000000..3acd5ddd2d0
--- /dev/null
+++ b/stan/math/prim/functor/reduce_sum.hpp
@@ -0,0 +1,220 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_REDUCE_SUM_HPP
+#define STAN_MATH_PRIM_FUNCTOR_REDUCE_SUM_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+
+#include <tbb/task_arena.h>
+#include <tbb/parallel_reduce.h>
+#include <tbb/blocked_range.h>
+
+#include <algorithm>
+#include <tuple>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+
+template <typename ReduceFunction, typename Enable, typename ReturnType,
+          typename Vec, typename... Args>
+struct reduce_sum_impl;
+
+/**
+ * Specialization of reduce_sum_impl for arithmetic types
+ *
+ * @tparam ReduceFunction Type of reducer function
+ * @tparam ReturnType An arithmetic type
+ * @tparam Vec Type of sliced argument
+ * @tparam Args Types of shared arguments
+ */
+template <typename ReduceFunction, typename ReturnType, typename Vec,
+          typename... Args>
+struct reduce_sum_impl<ReduceFunction, require_arithmetic_t<ReturnType>,
+                       ReturnType, Vec, Args...> {
+  /**
+   * This struct is used by the TBB to accumulate partial
+   *  sums over consecutive ranges of the input. To distribute the workload,
+   *  the TBB can split larger partial sums into smaller ones in which
+   *  case the splitting copy constructor is used. It is designed to
+   *  meet the Imperative form requirements of `tbb::parallel_reduce`.
+   *
+   * @note see link [here](https://tinyurl.com/vp7xw2t) for requirements.
+   */
+  struct recursive_reducer {
+    Vec vmapped_;
+    std::ostream* msgs_;
+    std::tuple<Args...> args_tuple_;
+    return_type_t<Vec, Args...> sum_{0.0};
+
+    recursive_reducer(Vec&& vmapped, std::ostream* msgs, Args&&... args)
+        : vmapped_(std::forward<Vec>(vmapped)),
+          msgs_(msgs),
+          args_tuple_(std::forward<Args>(args)...) {}
+
+    /**
+     * This is the copy operator as required for tbb::parallel_reduce
+     *   Imperative form. This requires sum_ be reset to zero since
+     *   the newly created reducer is used to accumulate an independent
+     *   partial sum.
+     */
+    recursive_reducer(recursive_reducer& other, tbb::split)
+        : vmapped_(other.vmapped_),
+          msgs_(other.msgs_),
+          args_tuple_(other.args_tuple_) {}
+
+    /**
+     * Compute the value and of `ReduceFunction` over the range defined by r
+     *   and accumulate those in member variable sum_. This function may
+     *   be called multiple times per object instantiation (so the sum_
+     *   must be accumulated, not just assigned).
+     *
+     * @param r Range over which to compute `ReduceFunction`
+     */
+    inline void operator()(const tbb::blocked_range<size_t>& r) {
+      if (r.empty()) {
+        return;
+      }
+
+      std::decay_t<Vec> sub_slice;
+      sub_slice.reserve(r.size());
+      for (size_t i = r.begin(); i < r.end(); ++i) {
+        sub_slice.emplace_back(vmapped_[i]);
+      }
+
+      sum_ += apply(
+          [&](auto&&... args) {
+            return ReduceFunction()(r.begin(), r.end() - 1, sub_slice, msgs_,
+                                    args...);
+          },
+          args_tuple_);
+    }
+
+    /**
+     * Join reducers. Accumuluate the value (sum_) of the other reducer.
+     *
+     * @param rhs Another partial sum
+     */
+    inline void join(const recursive_reducer& child) { sum_ += child.sum_; }
+  };
+
+  /**
+   * Call an instance of the function `ReduceFunction` on every element
+   *   of an input sequence and sum these terms.
+   *
+   * This specialization is parallelized using tbb and works only for
+   *   arithmetic types.
+   *
+   * ReduceFunction must define an operator() with the same signature as:
+   *   double f(int start, int end, Vec&& vmapped_subset, std::ostream* msgs,
+   * Args&&... args)
+   *
+   * `ReduceFunction` must be default constructible without any arguments
+   *
+   * Each call to `ReduceFunction` is responsible for computing the
+   *   start through end (inclusive) terms of the overall sum. All args are
+   * passed from this function through to the `ReduceFunction` instances.
+   *   However, only the start through end (inclusive) elements of the vmapped
+   * argument are passed to the `ReduceFunction` instances (as the
+   * `vmapped_subset` argument).
+   *
+   * This function distributes computation of the desired sum
+   *   over multiple threads by coordinating calls to `ReduceFunction`
+   * instances.
+   *
+   * If auto partitioning is true, break work into pieces automatically,
+   *  taking grainsize as a recommended work size. The partitioning is
+   *  not deterministic nor is the order guaranteed in which partial
+   *  sums are accumulated. Due to floating point imprecisions this will likely
+   *  lead to slight differences in the accumulated results between
+   *  multiple runs. If false, break work deterministically into pieces smaller
+   *  than or equal to grainsize and accumulate all the partial sums
+   *  in the same order. This still may not achieve bitwise reproducibility.
+   *
+   * grainsize must be greater than or equal to 1
+   *
+   * @param vmapped Vector containing one element per term of sum
+   * @param auto_partitioning Work partitioning style
+   * @param grainsize Suggested grainsize for tbb
+   * @param[in, out] msgs The print stream for warning messages
+   * @param args Shared arguments used in every sum term
+   * @return Summation of all terms
+   */
+  inline ReturnType operator()(Vec&& vmapped, bool auto_partitioning,
+                               int grainsize, std::ostream* msgs,
+                               Args&&... args) const {
+    const std::size_t num_terms = vmapped.size();
+    if (vmapped.empty()) {
+      return 0.0;
+    }
+    recursive_reducer worker(std::forward<Vec>(vmapped), msgs,
+                             std::forward<Args>(args)...);
+
+    if (auto_partitioning) {
+      tbb::parallel_reduce(
+          tbb::blocked_range<std::size_t>(0, num_terms, grainsize), worker);
+    } else {
+      tbb::simple_partitioner partitioner;
+      tbb::parallel_deterministic_reduce(
+          tbb::blocked_range<std::size_t>(0, num_terms, grainsize), worker,
+          partitioner);
+    }
+
+    return worker.sum_;
+  }
+};
+
+}  // namespace internal
+
+/**
+ * Call an instance of the function `ReduceFunction` on every element
+ *   of an input sequence and sum these terms.
+ *
+ * This defers to reduce_sum_impl for the appropriate implementation
+ *
+ * ReduceFunction must define an operator() with the same signature as:
+ *   T f(int start, int end, Vec&& vmapped_subset, std::ostream* msgs, Args&&...
+ * args)
+ *
+ * `ReduceFunction` must be default constructible without any arguments
+ *
+ * grainsize must be greater than or equal to 1
+ *
+ * @tparam ReduceFunction Type of reducer function
+ * @tparam ReturnType An arithmetic type
+ * @tparam Vec Type of sliced argument
+ * @tparam Args Types of shared arguments
+ * @param vmapped Vector containing one element per term of sum
+ * @param grainsize Suggested grainsize for tbb
+ * @param[in, out] msgs The print stream for warning messages
+ * @param args Shared arguments used in every sum term
+ * @return Sum of terms
+ */
+template <typename ReduceFunction, typename Vec,
+          typename = require_vector_like_t<Vec>, typename... Args>
+inline auto reduce_sum(Vec&& vmapped, int grainsize, std::ostream* msgs,
+                       Args&&... args) {
+  using return_type = return_type_t<Vec, Args...>;
+
+  check_positive("reduce_sum", "grainsize", grainsize);
+
+#ifdef STAN_THREADS
+  return internal::reduce_sum_impl<ReduceFunction, void, return_type, Vec,
+                                   Args...>()(std::forward<Vec>(vmapped), true,
+                                              grainsize, msgs,
+                                              std::forward<Args>(args)...);
+#else
+  if (vmapped.empty()) {
+    return return_type(0.0);
+  }
+
+  return ReduceFunction()(0, vmapped.size() - 1, std::forward<Vec>(vmapped),
+                          msgs, std::forward<Args>(args)...);
+#endif
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/functor/reduce_sum_static.hpp b/stan/math/prim/functor/reduce_sum_static.hpp
new file mode 100644
index 00000000000..0411a998b86
--- /dev/null
+++ b/stan/math/prim/functor/reduce_sum_static.hpp
@@ -0,0 +1,69 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_REDUCE_SUM_STATIC_HPP
+#define STAN_MATH_PRIM_FUNCTOR_REDUCE_SUM_STATIC_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/functor/reduce_sum.hpp>
+#include <tbb/task_arena.h>
+#include <tbb/parallel_reduce.h>
+#include <tbb/blocked_range.h>
+
+#include <tuple>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Call an instance of the function `ReduceFunction` on every element
+ *   of an input sequence and sum these terms.
+ *
+ * This defers to reduce_sum_impl for the appropriate implementation
+ *
+ * ReduceFunction must define an operator() with the same signature as:
+ *   T f(int start, int end, Vec&& vmapped_subset, std::ostream* msgs, Args&&...
+ * args)
+ *
+ * `ReduceFunction` must be default constructible without any arguments
+ *
+ * grainsize must be greater than or equal to 1
+ *
+ * If STAN_THREADS is not defined, do all the work with one ReduceFunction call.
+ *
+ * @tparam ReduceFunction Type of reducer function
+ * @tparam ReturnType An arithmetic type
+ * @tparam Vec Type of sliced argument
+ * @tparam Args Types of shared arguments
+ * @param vmapped Sliced arguments used only in some sum terms
+ * @param grainsize Suggested grainsize for tbb
+ * @param[in, out] msgs The print stream for warning messages
+ * @param args Shared arguments used in every sum term
+ * @return Sum of terms
+ */
+template <typename ReduceFunction, typename Vec,
+          typename = require_vector_like_t<Vec>, typename... Args>
+auto reduce_sum_static(Vec&& vmapped, int grainsize, std::ostream* msgs,
+                       Args&&... args) {
+  using return_type = return_type_t<Vec, Args...>;
+
+  check_positive("reduce_sum", "grainsize", grainsize);
+
+#ifdef STAN_THREADS
+  return internal::reduce_sum_impl<ReduceFunction, void, return_type, Vec,
+                                   Args...>()(std::forward<Vec>(vmapped), false,
+                                              grainsize, msgs,
+                                              std::forward<Args>(args)...);
+#else
+  if (vmapped.empty()) {
+    return return_type(0);
+  }
+
+  return ReduceFunction()(0, vmapped.size() - 1, std::forward<Vec>(vmapped),
+                          msgs, std::forward<Args>(args)...);
+#endif
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/functor.hpp b/stan/math/rev/functor.hpp
index 0af65fe88e5..13566d8b8bc 100644
--- a/stan/math/rev/functor.hpp
+++ b/stan/math/rev/functor.hpp
@@ -20,5 +20,6 @@
 #include <stan/math/rev/functor/kinsol_solve.hpp>
 #include <stan/math/rev/functor/map_rect_concurrent.hpp>
 #include <stan/math/rev/functor/map_rect_reduce.hpp>
+#include <stan/math/rev/functor/reduce_sum.hpp>
 
 #endif
diff --git a/stan/math/rev/functor/reduce_sum.hpp b/stan/math/rev/functor/reduce_sum.hpp
new file mode 100644
index 00000000000..4f0d3300fd2
--- /dev/null
+++ b/stan/math/rev/functor/reduce_sum.hpp
@@ -0,0 +1,251 @@
+#ifndef STAN_MATH_REV_FUNCTOR_REDUCE_SUM_HPP
+#define STAN_MATH_REV_FUNCTOR_REDUCE_SUM_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/functor.hpp>
+#include <stan/math/rev/core.hpp>
+#include <tbb/task_arena.h>
+#include <tbb/parallel_reduce.h>
+#include <tbb/blocked_range.h>
+
+#include <tuple>
+#include <vector>
+
+namespace stan {
+namespace math {
+namespace internal {
+
+/**
+ * Var specialization of reduce_sum_impl
+ *
+ * @tparam ReduceFunction Type of reducer function
+ * @tparam ReturnType Must be var
+ * @tparam Vec Type of sliced argument
+ * @tparam Args Types of shared arguments
+ */
+template <typename ReduceFunction, typename ReturnType, typename Vec,
+          typename... Args>
+struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
+                       Vec, Args...> {
+  /**
+   * This struct is used by the TBB to accumulate partial
+   *  sums over consecutive ranges of the input. To distribute the workload,
+   *  the TBB can split larger partial sums into smaller ones in which
+   *  case the splitting copy constructor is used. It is designed to
+   *  meet the Imperative form requirements of `tbb::parallel_reduce`.
+   *
+   * @note see link [here](https://tinyurl.com/vp7xw2t) for requirements.
+   */
+  struct recursive_reducer {
+    const size_t num_vars_per_term_;
+    const size_t num_vars_shared_terms_;  // Number of vars in shared arguments
+    double* sliced_partials_;  // Points to adjoints of the partial calculations
+    Vec vmapped_;
+    std::ostream* msgs_;
+    std::tuple<Args...> args_tuple_;
+    double sum_{0.0};
+    Eigen::VectorXd args_adjoints_{0};
+
+    template <typename VecT, typename... ArgsT>
+    recursive_reducer(size_t num_vars_per_term, size_t num_vars_shared_terms,
+                      double* sliced_partials, VecT&& vmapped,
+                      std::ostream* msgs, ArgsT&&... args)
+        : num_vars_per_term_(num_vars_per_term),
+          num_vars_shared_terms_(num_vars_shared_terms),
+          sliced_partials_(sliced_partials),
+          vmapped_(std::forward<VecT>(vmapped)),
+          msgs_(msgs),
+          args_tuple_(std::forward<ArgsT>(args)...) {}
+
+    /*
+     * This is the copy operator as required for tbb::parallel_reduce
+     *   Imperative form. This requires sum_ and arg_adjoints_ be reset
+     *   to zero since the newly created reducer is used to accumulate
+     *   an independent partial sum.
+     */
+    recursive_reducer(recursive_reducer& other, tbb::split)
+        : num_vars_per_term_(other.num_vars_per_term_),
+          num_vars_shared_terms_(other.num_vars_shared_terms_),
+          sliced_partials_(other.sliced_partials_),
+          vmapped_(other.vmapped_),
+          msgs_(other.msgs_),
+          args_tuple_(other.args_tuple_) {}
+
+    /**
+     * Compute, using nested autodiff, the value and Jacobian of
+     *  `ReduceFunction` called over the range defined by r and accumulate those
+     *  in member variable sum_ (for the value) and args_adjoints_ (for the
+     *  Jacobian). The nested autodiff uses deep copies of the involved operands
+     *  ensuring that no side effects are implied to the adjoints of the input
+     *  operands which reside potentially on a autodiff tape stored in a
+     *  different thread other than the current thread of execution. This
+     * function may be called multiple times per object instantiation (so the
+     * sum_ and args_adjoints_ must be accumulated, not just assigned).
+     *
+     * @param r Range over which to compute reduce_sum
+     */
+    inline void operator()(const tbb::blocked_range<size_t>& r) {
+      if (r.empty()) {
+        return;
+      }
+
+      if (args_adjoints_.size() == 0) {
+        args_adjoints_ = Eigen::VectorXd::Zero(num_vars_shared_terms_);
+      }
+
+      // Initialize nested autodiff stack
+      const nested_rev_autodiff begin_nest;
+
+      // Create nested autodiff copies of sliced argument that do not point
+      //   back to main autodiff stack
+      std::decay_t<Vec> local_sub_slice;
+      local_sub_slice.reserve(r.size());
+      for (size_t i = r.begin(); i < r.end(); ++i) {
+        local_sub_slice.emplace_back(deep_copy_vars(vmapped_[i]));
+      }
+
+      // Create nested autodiff copies of all shared arguments that do not point
+      //   back to main autodiff stack
+      auto args_tuple_local_copy = apply(
+          [&](auto&&... args) {
+            return std::tuple<decltype(deep_copy_vars(args))...>(
+                deep_copy_vars(args)...);
+          },
+          args_tuple_);
+
+      // Perform calculation
+      var sub_sum_v = apply(
+          [&](auto&&... args) {
+            return ReduceFunction()(r.begin(), r.end() - 1, local_sub_slice,
+                                    msgs_, args...);
+          },
+          args_tuple_local_copy);
+
+      // Compute Jacobian
+      sub_sum_v.grad();
+
+      // Accumulate value of reduce_sum
+      sum_ += sub_sum_v.val();
+
+      // Accumulate adjoints of sliced_arguments
+      accumulate_adjoints(sliced_partials_ + r.begin() * num_vars_per_term_,
+                          std::move(local_sub_slice));
+
+      // Accumulate adjoints of shared_arguments
+      apply(
+          [&](auto&&... args) {
+            accumulate_adjoints(args_adjoints_.data(),
+                                std::forward<decltype(args)>(args)...);
+          },
+          std::move(args_tuple_local_copy));
+    }
+
+    /**
+     * Join reducers. Accumuluate the value (sum_) and Jacobian (arg_adoints_)
+     *   of the other reducer.
+     *
+     * @param rhs Another partial sum
+     */
+    inline void join(const recursive_reducer& rhs) {
+      sum_ += rhs.sum_;
+      if (args_adjoints_.size() != 0 && rhs.args_adjoints_.size() != 0) {
+        args_adjoints_ += rhs.args_adjoints_;
+      } else if (args_adjoints_.size() == 0 && rhs.args_adjoints_.size() != 0) {
+        args_adjoints_ = rhs.args_adjoints_;
+      }
+    }
+  };
+
+  /**
+   * Call an instance of the function `ReduceFunction` on every element
+   *   of an input sequence and sum these terms.
+   *
+   * This specialization is parallelized using tbb and works for reverse
+   *   mode autodiff.
+   *
+   * ReduceFunction must define an operator() with the same signature as:
+   *   var f(int start, int end, Vec&& vmapped_subset, std::ostream* msgs,
+   * Args&&... args)
+   *
+   * `ReduceFunction` must be default constructible without any arguments
+   *
+   * Each call to `ReduceFunction` is responsible for computing the
+   *   start through end (inclusive) terms of the overall sum. All args are
+   * passed from this function through to the `ReduceFunction` instances.
+   *   However, only the start through end (inclusive) elements of the vmapped
+   * argument are passed to the `ReduceFunction` instances (as the
+   * `vmapped_subset` argument).
+   *
+   * This function distributes computation of the desired sum and the Jacobian
+   * of that sum over multiple threads by coordinating calls to `ReduceFunction`
+   * instances. Results are stored as precomputed varis in the autodiff tree.
+   *
+   * If auto partitioning is true, break work into pieces automatically,
+   *  taking grainsize as a recommended work size. The partitioning is
+   *  not deterministic nor is the order guaranteed in which partial
+   *  sums are accumulated. Due to floating point imprecisions this will likely
+   *  lead to slight differences in the accumulated results between
+   *  multiple runs. If false, break work deterministically into pieces smaller
+   *  than or equal to grainsize and accumulate all the partial sums
+   *  in the same order. This still may not achieve bitwise reproducibility.
+   *
+   * @param vmapped Vector containing one element per term of sum
+   * @param auto_partitioning Work partitioning style
+   * @param grainsize Suggested grainsize for tbb
+   * @param[in, out] msgs The print stream for warning messages
+   * @param args Shared arguments used in every sum term
+   * @return Summation of all terms
+   */
+  inline var operator()(Vec&& vmapped, bool auto_partitioning, int grainsize,
+                        std::ostream* msgs, Args&&... args) const {
+    const std::size_t num_terms = vmapped.size();
+
+    if (vmapped.empty()) {
+      return var(0.0);
+    }
+
+    const std::size_t num_vars_per_term = count_vars(vmapped[0]);
+    const std::size_t num_vars_sliced_terms = num_terms * num_vars_per_term;
+    const std::size_t num_vars_shared_terms = count_vars(args...);
+
+    vari** varis = ChainableStack::instance_->memalloc_.alloc_array<vari*>(
+        num_vars_sliced_terms + num_vars_shared_terms);
+    double* partials = ChainableStack::instance_->memalloc_.alloc_array<double>(
+        num_vars_sliced_terms + num_vars_shared_terms);
+
+    save_varis(varis, vmapped);
+    save_varis(varis + num_vars_sliced_terms, args...);
+
+    for (size_t i = 0; i < num_vars_sliced_terms; ++i) {
+      partials[i] = 0.0;
+    }
+
+    recursive_reducer worker(num_vars_per_term, num_vars_shared_terms, partials,
+                             std::forward<Vec>(vmapped), msgs,
+                             std::forward<Args>(args)...);
+
+    if (auto_partitioning) {
+      tbb::parallel_reduce(
+          tbb::blocked_range<std::size_t>(0, num_terms, grainsize), worker);
+    } else {
+      tbb::simple_partitioner partitioner;
+      tbb::parallel_deterministic_reduce(
+          tbb::blocked_range<std::size_t>(0, num_terms, grainsize), worker,
+          partitioner);
+    }
+
+    for (size_t i = 0; i < num_vars_shared_terms; ++i) {
+      partials[num_vars_sliced_terms + i] = worker.args_adjoints_(i);
+    }
+
+    return var(new precomputed_gradients_vari(
+        worker.sum_, num_vars_sliced_terms + num_vars_shared_terms, varis,
+        partials));
+  }
+};
+}  // namespace internal
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/mix/functor/reduce_sum_part1_test.cpp b/test/unit/math/mix/functor/reduce_sum_part1_test.cpp
new file mode 100644
index 00000000000..abf7bc184f3
--- /dev/null
+++ b/test/unit/math/mix/functor/reduce_sum_part1_test.cpp
@@ -0,0 +1,215 @@
+#include <stan/math/prim/meta.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/functor/reduce_sum_util.hpp>
+
+#include <limits>
+#include <vector>
+
+// Reduce sum tests are broken up into four files to avoid windows compiler
+// error
+
+TEST(MathMix_reduce_sum, grainsize_static) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::sum_lpdf;
+
+  auto f1 = [](auto&& data) {
+    return stan::math::reduce_sum_static<sum_lpdf>(data, 0, get_new_msg());
+  };
+  auto f2 = [](auto&& data) {
+    return stan::math::reduce_sum_static<sum_lpdf>(data, -1, get_new_msg());
+  };
+  auto f3 = [](auto&& data) {
+    return stan::math::reduce_sum_static<sum_lpdf>(data, 1, get_new_msg());
+  };
+  auto f4 = [](auto&& data) {
+    return stan::math::reduce_sum_static<sum_lpdf>(data, 20, get_new_msg());
+  };
+
+  std::vector<double> data(5, 10.0);
+
+  stan::test::expect_ad(f1, data);
+  stan::test::expect_ad(f2, data);
+  stan::test::expect_ad(f3, data);
+  stan::test::expect_ad(f4, data);
+}
+
+TEST(MathMix_reduce_sum, grainsize) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::sum_lpdf;
+  auto f1 = [](auto&& data) {
+    return stan::math::reduce_sum<sum_lpdf>(data, 0, get_new_msg());
+  };
+  auto f2 = [](auto&& data) {
+    return stan::math::reduce_sum<sum_lpdf>(data, -1, get_new_msg());
+  };
+  auto f3 = [](auto&& data) {
+    return stan::math::reduce_sum<sum_lpdf>(data, 1, get_new_msg());
+  };
+  auto f4 = [](auto&& data) {
+    return stan::math::reduce_sum<sum_lpdf>(data, 20, get_new_msg());
+  };
+
+  std::vector<double> data(5, 10.0);
+
+  stan::test::expect_ad(f1, data);
+  stan::test::expect_ad(f2, data);
+  stan::test::expect_ad(f3, data);
+  stan::test::expect_ad(f4, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_zero_length) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<double> data(0);
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<double> data(5, 10.0);
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_std_vector_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<std::vector<double>> data(3, std::vector<double>(2, 10.0));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_eigen_vector_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<Eigen::VectorXd> data(3, Eigen::VectorXd::Ones(2));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_eigen_row_vector_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<Eigen::RowVectorXd> data(3, Eigen::RowVectorXd::Ones(2));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_eigen_matrix_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<Eigen::MatrixXd> data(3, Eigen::MatrixXd::Ones(2, 4));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_std_vector_std_vector_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<std::vector<std::vector<double>>> data(
+      3, std::vector<std::vector<double>>(2, std::vector<double>(2, 10.0)));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_std_vector_eigen_vector_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<std::vector<Eigen::VectorXd>> data(
+      3, std::vector<Eigen::VectorXd>(2, Eigen::VectorXd::Ones(2)));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_std_vector_eigen_row_vector_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<std::vector<Eigen::RowVectorXd>> data(
+      3, std::vector<Eigen::RowVectorXd>(2, Eigen::RowVectorXd::Ones(2)));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_std_vector_eigen_matrix_double_slice) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+
+  std::vector<std::vector<Eigen::MatrixXd>> data(
+      3, std::vector<Eigen::MatrixXd>(2, Eigen::MatrixXd::Ones(2, 4)));
+
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data);
+}
+
+TEST(StanMath_reduce_sum_static, start_end_slice) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::start_end_lpdf;
+  auto start_end_static = [](auto&& arg) {
+    return stan::math::reduce_sum_static<start_end_lpdf>(arg, 1, get_new_msg(),
+                                                         arg);
+  };
+
+  auto start_end = [](auto&& arg) {
+    return stan::math::reduce_sum_static<start_end_lpdf>(arg, 1, get_new_msg(),
+                                                         arg);
+  };
+
+  std::vector<double> data(5, 1.0);
+
+  stan::test::expect_ad(start_end, data);
+  stan::test::expect_ad(start_end_static, data);
+}
+
+TEST(MathMix_reduce_sum, int_arg) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+  std::vector<double> data(2, 1.0);
+  int arg = 5.0;
+
+  stan::test::expect_ad(
+      [&](auto&& data) { return reduce_sum_static_sum_lpdf(data, arg); }, data);
+  stan::test::expect_ad(
+      [&](auto&& data) { return reduce_sum_sum_lpdf(data, arg); }, data);
+}
+
+TEST(MathMix_reduce_sum, std_vector_int_arg) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+  std::vector<double> data(2, 10.0);
+  std::vector<int> arg(2, 10);
+
+  stan::test::expect_ad(
+      [&](auto&& data) { return reduce_sum_static_sum_lpdf(data, arg); }, data);
+  stan::test::expect_ad(
+      [&](auto&& data) { return reduce_sum_sum_lpdf(data, arg); }, data);
+}
+
+TEST(MathMix_reduce_sum, double_arg) {
+  stan::math::test::expect_ad_reduce_sum_lpdf(std::vector<double>(2, 10.0),
+                                              5.0);
+}
+
+TEST(MathMix_reduce_sum, std_vector_double_arg) {
+  stan::math::test::expect_ad_reduce_sum_lpdf(std::vector<double>(2, 10.0),
+                                              std::vector<double>(2, 10.0));
+}
diff --git a/test/unit/math/mix/functor/reduce_sum_part2_test.cpp b/test/unit/math/mix/functor/reduce_sum_part2_test.cpp
new file mode 100644
index 00000000000..524671b3223
--- /dev/null
+++ b/test/unit/math/mix/functor/reduce_sum_part2_test.cpp
@@ -0,0 +1,50 @@
+#include <stan/math/prim/meta.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/functor/reduce_sum_util.hpp>
+
+#include <limits>
+#include <vector>
+
+// Reduce sum tests are broken up into four files to avoid windows compiler
+// error
+
+TEST(MathMix_reduce_sum, eigen_vector_arg) {
+  std::vector<double> data(2, 10.0);
+  Eigen::VectorXd arg = Eigen::VectorXd::Ones(2);
+  stan::math::test::expect_ad_reduce_sum_lpdf(data, arg);
+}
+
+TEST(MathMix_reduce_sum, eigen_row_vector_arg) {
+  std::vector<double> data(2, 10.0);
+  Eigen::RowVectorXd arg = Eigen::RowVectorXd::Ones(2);
+
+  stan::math::test::expect_ad_reduce_sum_lpdf(data, arg);
+}
+
+TEST(MathMix_reduce_sum, eigen_matrix_arg) {
+  std::vector<double> data(2, 10.0);
+  Eigen::MatrixXd arg = Eigen::MatrixXd::Ones(2, 2);
+
+  stan::math::test::expect_ad_reduce_sum_lpdf(data, arg);
+}
+
+TEST(MathMix_reduce_sum, std_vector_std_vector_double_arg) {
+  std::vector<double> data(2, 10.0);
+  std::vector<std::vector<double>> arg(2, std::vector<double>(2, 10.0));
+
+  stan::math::test::expect_ad_reduce_sum_lpdf(data, arg);
+}
+
+TEST(MathMix_reduce_sum, std_vector_eigen_vector_arg) {
+  std::vector<double> data(2, 10.0);
+  std::vector<Eigen::VectorXd> arg(2, Eigen::VectorXd::Ones(2));
+
+  stan::math::test::expect_ad_reduce_sum_lpdf(data, arg);
+}
+
+TEST(MathMix_reduce_sum, std_vector_eigen_row_vector_arg) {
+  std::vector<double> data(2, 10.0);
+  std::vector<Eigen::RowVectorXd> arg(2, Eigen::RowVectorXd::Ones(2));
+
+  stan::math::test::expect_ad_reduce_sum_lpdf(data, arg);
+}
diff --git a/test/unit/math/mix/functor/reduce_sum_part3_test.cpp b/test/unit/math/mix/functor/reduce_sum_part3_test.cpp
new file mode 100644
index 00000000000..52723cdeee3
--- /dev/null
+++ b/test/unit/math/mix/functor/reduce_sum_part3_test.cpp
@@ -0,0 +1,85 @@
+#include <stan/math/prim/meta.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/functor/reduce_sum_util.hpp>
+
+#include <limits>
+#include <vector>
+
+// Reduce sum tests are broken up into four files to avoid windows compiler
+// error
+
+TEST(MathMix_reduce_sum, eigen_three_args1) {
+  using stan::math::test::reduce_sum_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_int_sum_lpdf;
+  Eigen::VectorXd arg1 = Eigen::VectorXd::Ones(2);
+  Eigen::RowVectorXd arg2 = Eigen::RowVectorXd::Ones(2);
+  Eigen::MatrixXd arg3 = Eigen::MatrixXd::Ones(2, 2);
+
+  stan::test::expect_ad(reduce_sum_static_int_sum_lpdf, arg1, arg2, arg3);
+  stan::test::expect_ad(reduce_sum_int_sum_lpdf, arg1, arg2, arg3);
+}
+
+TEST(MathMix_reduce_sum, eigen_three_args2) {
+  using stan::math::test::reduce_sum_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_int_sum_lpdf;
+  double arg1 = 1.0;
+  std::vector<double> arg2(2, 1.0);
+  Eigen::MatrixXd arg3 = Eigen::MatrixXd::Ones(2, 2);
+
+  stan::test::expect_ad(reduce_sum_static_int_sum_lpdf, arg1, arg2, arg3);
+  stan::test::expect_ad(reduce_sum_int_sum_lpdf, arg1, arg2, arg3);
+}
+
+TEST(MathMix_reduce_sum, eigen_three_args3) {
+  using stan::math::test::reduce_sum_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_int_sum_lpdf;
+  double arg1 = 1.0;
+  std::vector<std::vector<double>> arg2(2, std::vector<double>(2, 1.0));
+  std::vector<Eigen::MatrixXd> arg3(2, Eigen::MatrixXd::Ones(2, 2));
+
+  stan::test::expect_ad(reduce_sum_static_int_sum_lpdf, arg1, arg2, arg3);
+  stan::test::expect_ad(reduce_sum_int_sum_lpdf, arg1, arg2, arg3);
+}
+
+TEST(MathMix_reduce_sum, eigen_three_args_with_ints1) {
+  using stan::math::test::reduce_sum_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_int_sum_lpdf;
+  Eigen::VectorXd arg1 = Eigen::VectorXd::Ones(2);
+  Eigen::RowVectorXd arg2 = Eigen::RowVectorXd::Ones(2);
+  Eigen::MatrixXd arg3 = Eigen::MatrixXd::Ones(2, 2);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_static_int_sum_lpdf(
+            1, arg1, std::vector<int>{1, 2, 3}, arg2, 3, arg3);
+      },
+      arg1, arg2, arg3);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_int_sum_lpdf(1, arg1, std::vector<int>{1, 2, 3}, arg2,
+                                       3, arg3);
+      },
+      arg1, arg2, arg3);
+}
+
+TEST(MathMix_reduce_sum, eigen_three_args_with_ints2) {
+  using stan::math::test::reduce_sum_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_int_sum_lpdf;
+  double arg1 = 1.0;
+  std::vector<double> arg2(2, 1.0);
+  Eigen::MatrixXd arg3 = Eigen::MatrixXd::Ones(2, 2);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_static_int_sum_lpdf(
+            1, arg1, std::vector<int>{1, 2, 3}, arg2, 3, arg3);
+      },
+      arg1, arg2, arg3);
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_int_sum_lpdf(1, arg1, std::vector<int>{1, 2, 3}, arg2,
+                                       3, arg3);
+      },
+      arg1, arg2, arg3);
+}
diff --git a/test/unit/math/mix/functor/reduce_sum_part4_test.cpp b/test/unit/math/mix/functor/reduce_sum_part4_test.cpp
new file mode 100644
index 00000000000..e2cfead9893
--- /dev/null
+++ b/test/unit/math/mix/functor/reduce_sum_part4_test.cpp
@@ -0,0 +1,133 @@
+#include <stan/math/prim/meta.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/mix/functor/reduce_sum_util.hpp>
+
+#include <limits>
+#include <vector>
+
+// Reduce sum tests are broken up into four files to avoid windows compiler
+// error
+
+TEST(MathMix_reduce_sum, eigen_three_args_with_ints3) {
+  using stan::math::test::reduce_sum_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_int_sum_lpdf;
+  double arg1 = 1.0;
+  std::vector<std::vector<double>> arg2(2, std::vector<double>(2, 1.0));
+  std::vector<Eigen::MatrixXd> arg3(2, Eigen::MatrixXd::Ones(2, 2));
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_static_int_sum_lpdf(
+            1, arg1, std::vector<int>{1, 2, 3}, arg2, 3, arg3);
+      },
+      arg1, arg2, arg3);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_int_sum_lpdf(1, arg1, std::vector<int>{1, 2, 3}, arg2,
+                                       3, arg3);
+      },
+      arg1, arg2, arg3);
+}
+
+TEST(MathMix_reduce_sum, eigen_three_args_with_doubles1) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+  Eigen::VectorXd arg1 = Eigen::VectorXd::Ones(2);
+  Eigen::RowVectorXd arg2 = Eigen::RowVectorXd::Ones(2);
+  Eigen::MatrixXd arg3 = Eigen::MatrixXd::Ones(2, 2);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_static_sum_lpdf(
+            std::vector<double>{1.0, 2.0, 3.0}, arg1, 3.0, arg2,
+            std::vector<double>{1.0, 2.0, 3.0}, arg3);
+      },
+      arg1, arg2, arg3);
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_sum_lpdf(std::vector<double>{1.0, 2.0, 3.0}, arg1,
+                                   3.0, arg2,
+                                   std::vector<double>{1.0, 2.0, 3.0}, arg3);
+      },
+      arg1, arg2, arg3);
+}
+
+TEST(MathMix_reduce_sum, eigen_three_args_with_doubles2) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+  double arg1 = 1.0;
+  std::vector<double> arg2(2, 1.0);
+  Eigen::MatrixXd arg3 = Eigen::MatrixXd::Ones(2, 2);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_static_sum_lpdf(
+            std::vector<double>{1.0, 2.0, 3.0}, arg1, 3.0, arg2,
+            std::vector<double>{1.0, 2.0, 3.0}, arg3);
+      },
+      arg1, arg2, arg3);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_sum_lpdf(std::vector<double>{1.0, 2.0, 3.0}, arg1,
+                                   3.0, arg2,
+                                   std::vector<double>{1.0, 2.0, 3.0}, arg3);
+      },
+      arg1, arg2, arg3);
+}
+
+TEST(MathMix_reduce_sum, eigen_three_args_with_doubles3) {
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+  double arg1 = 1.0;
+  std::vector<std::vector<double>> arg2(2, std::vector<double>(2, 1.0));
+  std::vector<Eigen::MatrixXd> arg3(2, Eigen::MatrixXd::Ones(2, 2));
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_static_sum_lpdf(
+            std::vector<double>{1.0, 2.0, 3.0}, arg1, 3.0, arg2,
+            std::vector<double>{1.0, 2.0, 3.0}, arg3);
+      },
+      arg1, arg2, arg3);
+
+  stan::test::expect_ad(
+      [&](auto&& arg1, auto&& arg2, auto&& arg3) {
+        return reduce_sum_sum_lpdf(std::vector<double>{1.0, 2.0, 3.0}, arg1,
+                                   3.0, arg2,
+                                   std::vector<double>{1.0, 2.0, 3.0}, arg3);
+      },
+      arg1, arg2, arg3);
+}
+
+#ifdef STAN_THREADS
+TEST(MathMix_reduce_sum, static_check) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::static_check_lpdf;
+
+  for (auto size : {1, 3, 6, 11}) {
+    std::vector<int> data(size, 10);
+    std::vector<double> arg(size, 10);
+
+    auto fi1 = [&](auto&&... args) {
+      return stan::math::reduce_sum_static<static_check_lpdf<1>>(
+          data, 1, get_new_msg(), args...);
+    };
+
+    auto fi2 = [&](auto&&... args) {
+      return stan::math::reduce_sum_static<static_check_lpdf<2>>(
+          data, 2, get_new_msg(), args...);
+    };
+
+    auto fi3 = [&](auto&&... args) {
+      return stan::math::reduce_sum_static<static_check_lpdf<3>>(
+          data, 3, get_new_msg(), args...);
+    };
+
+    stan::test::expect_ad(fi1, arg);
+    stan::test::expect_ad(fi2, arg);
+    stan::test::expect_ad(fi3, arg);
+  }
+}
+#endif
diff --git a/test/unit/math/mix/functor/reduce_sum_util.hpp b/test/unit/math/mix/functor/reduce_sum_util.hpp
new file mode 100644
index 00000000000..9d70fe1054c
--- /dev/null
+++ b/test/unit/math/mix/functor/reduce_sum_util.hpp
@@ -0,0 +1,32 @@
+#ifndef TEST_UNIT_MATH_MIX_FUNCTOR_REDUCE_SUM_UTIL
+#define TEST_UNIT_MATH_MIX_FUNCTOR_REDUCE_SUM_UTIL
+
+#include <stan/math.hpp>
+#include <test/unit/math/prim/functor/reduce_sum_util.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <gtest/gtest.h>
+#include <algorithm>
+#include <sstream>
+#include <tuple>
+#include <vector>
+
+namespace stan {
+namespace math {
+namespace test {
+
+template <typename T1, typename T2>
+void expect_ad_reduce_sum_lpdf(T1&& data, T2&& arg) {
+  using stan::math::test::reduce_sum_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_int_sum_lpdf;
+  using stan::math::test::reduce_sum_static_sum_lpdf;
+  using stan::math::test::reduce_sum_sum_lpdf;
+  stan::test::expect_ad(reduce_sum_static_int_sum_lpdf, arg);
+  stan::test::expect_ad(reduce_sum_static_sum_lpdf, data, arg);
+  stan::test::expect_ad(reduce_sum_int_sum_lpdf, arg);
+  stan::test::expect_ad(reduce_sum_sum_lpdf, data, arg);
+}
+
+}  // namespace test
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/unit/math/prim/functor/reduce_sum_test.cpp b/test/unit/math/prim/functor/reduce_sum_test.cpp
new file mode 100644
index 00000000000..463f29a38a9
--- /dev/null
+++ b/test/unit/math/prim/functor/reduce_sum_test.cpp
@@ -0,0 +1,268 @@
+#include <stan/math.hpp>
+#include <test/unit/math/prim/functor/reduce_sum_util.hpp>
+#include <gtest/gtest.h>
+
+#include <vector>
+#include <set>
+
+TEST(StanMathPrim_reduce_sum, value) {
+  using stan::math::test::count_lpdf;
+  using stan::math::test::get_new_msg;
+  double lambda_d = 10.0;
+  const std::size_t elems = 10000;
+  const std::size_t num_iter = 1000;
+  std::vector<int> data(elems);
+
+  for (std::size_t i = 0; i != elems; ++i)
+    data[i] = i;
+
+  std::vector<int> idata;
+  std::vector<double> vlambda_d(1, lambda_d);
+
+  double poisson_lpdf = stan::math::reduce_sum<count_lpdf<double>>(
+      data, 5, get_new_msg(), vlambda_d, idata);
+  double poisson_static_lpdf
+      = stan::math::reduce_sum_static<count_lpdf<double>>(
+          data, 5, get_new_msg(), vlambda_d, idata);
+
+  double poisson_lpdf_ref = stan::math::poisson_lpmf(data, lambda_d);
+  // NOTE:(Steve) This fails with EXPECT_DOUBLE_EQ at about 10e-7
+  EXPECT_FLOAT_EQ(poisson_lpdf, poisson_lpdf_ref);
+  EXPECT_FLOAT_EQ(poisson_static_lpdf, poisson_lpdf_ref);
+}
+
+TEST(StanMathPrim_reduce_sum, grainsize) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::sum_lpdf;
+
+  std::vector<int> data(5, 10);
+
+  EXPECT_THROW(stan::math::reduce_sum<sum_lpdf>(data, 0, get_new_msg()),
+               std::domain_error);
+
+  EXPECT_THROW(stan::math::reduce_sum<sum_lpdf>(data, -1, get_new_msg()),
+               std::domain_error);
+
+  EXPECT_NO_THROW(stan::math::reduce_sum<sum_lpdf>(data, 1, get_new_msg()));
+
+  EXPECT_NO_THROW(
+      stan::math::reduce_sum<sum_lpdf>(data, 3 * data.size(), get_new_msg()));
+  EXPECT_THROW(stan::math::reduce_sum_static<sum_lpdf>(data, 0, get_new_msg()),
+               std::domain_error);
+
+  EXPECT_THROW(stan::math::reduce_sum_static<sum_lpdf>(data, -1, get_new_msg()),
+               std::domain_error);
+
+  EXPECT_NO_THROW(
+      stan::math::reduce_sum_static<sum_lpdf>(data, 1, get_new_msg()));
+
+  EXPECT_NO_THROW(stan::math::reduce_sum_static<sum_lpdf>(data, 3 * data.size(),
+                                                          get_new_msg()));
+}
+
+TEST(StanMathPrim_reduce_sum, start_end_slice) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::start_end_lpdf;
+
+  std::vector<int> data(5, 10);
+
+  EXPECT_EQ(
+      50, stan::math::reduce_sum<start_end_lpdf>(data, 1, get_new_msg(), data));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_int_slice) {
+  stan::math::test::test_slices(50, 10);
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_double_slice) {
+  stan::math::test::test_slices(50.0, 10.0);
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_double_slice) {
+  stan::math::test::test_slices(100.0, std::vector<double>(2, 10));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_int_slice) {
+  stan::math::test::test_slices(100, std::vector<int>(2, 10));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_eigen_vector_slice) {
+  Eigen::VectorXd test_x = Eigen::VectorXd::Ones(2);
+  stan::math::test::test_slices(10.0, test_x);
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_eigen_row_vector_slice) {
+  Eigen::RowVectorXd test_x = Eigen::RowVectorXd::Ones(2);
+  stan::math::test::test_slices(10.0, test_x);
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_eigen_matrix_slice) {
+  Eigen::MatrixXd test_x = Eigen::MatrixXd::Ones(2, 2);
+  stan::math::test::test_slices(20.0, test_x);
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_std_vector_int_slice) {
+  stan::math::test::test_slices(
+      200, std::vector<std::vector<int>>(2, std::vector<int>(2, 10.0)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_std_vector_double_slice) {
+  stan::math::test::test_slices(
+      200.0, std::vector<std::vector<double>>(2, std::vector<double>(2, 10.0)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_eigen_vector_slice) {
+  stan::math::test::test_slices(
+      20.0, std::vector<Eigen::VectorXd>(2, Eigen::VectorXd::Ones(2)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_eigen_row_vector_slice) {
+  stan::math::test::test_slices(
+      20.0, std::vector<Eigen::RowVectorXd>(2, Eigen::RowVectorXd::Ones(2)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_eigen_matrix_slice) {
+  stan::math::test::test_slices(
+      40.0, std::vector<Eigen::MatrixXd>(2, Eigen::MatrixXd::Ones(2, 2)));
+}
+
+TEST(StanMathPrim_reduce_sum, no_args) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::sum_lpdf;
+
+  std::vector<double> data(0);
+  EXPECT_EQ(0.0, stan::math::reduce_sum_static<sum_lpdf>(
+                     data, 1, stan::math::test::get_new_msg()))
+      << "Failed for reduce_sum_static";
+  EXPECT_EQ(0.0, stan::math::reduce_sum<sum_lpdf>(
+                     data, 1, stan::math::test::get_new_msg()))
+      << "Failed for reduce_sum";
+}
+
+TEST(StanMathPrim_reduce_sum, int_arg) {
+  stan::math::test::test_slices(5 * (10 + 5), 10, 5);
+}
+
+TEST(StanMathPrim_reduce_sum, double_arg) {
+  stan::math::test::test_slices(5 * (10.0 + 5.0), 10.0, 5.0);
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_int_arg) {
+  stan::math::test::test_slices(5 * (10.0 + 5.0 * 10), 10.0,
+                                std::vector<int>(5, 10));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_double_arg) {
+  stan::math::test::test_slices(5 * (10 + 5 * 10), 10.0,
+                                std::vector<double>(5, 10));
+}
+
+TEST(StanMathPrim_reduce_sum, eigen_vector_arg) {
+  stan::math::test::test_slices(5 * (10 + 5), 10.0, Eigen::VectorXd::Ones(5));
+}
+
+TEST(StanMathPrim_reduce_sum, eigen_row_vector_arg) {
+  stan::math::test::test_slices(5 * (10 + 5), 10.0,
+                                Eigen::RowVectorXd::Ones(5));
+}
+
+TEST(StanMathPrim_reduce_sum, eigen_matrix_arg) {
+  stan::math::test::test_slices(5 * (10 + 5 * 5), 10.0,
+                                Eigen::MatrixXd::Ones(5, 5));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_double_arg) {
+  stan::math::test::test_slices(
+      5 * (10 + 250), 10.0,
+      std::vector<std::vector<double>>(5, std::vector<double>(5, 10.0)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_eigen_vector_arg) {
+  stan::math::test::test_slices(
+      5 * (10 + 10), 10.0,
+      std::vector<Eigen::VectorXd>(2, Eigen::VectorXd::Ones(5)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_eigen_row_vector_arg) {
+  stan::math::test::test_slices(
+      5 * (10 + 10), 10.0,
+      std::vector<Eigen::RowVectorXd>(2, Eigen::RowVectorXd::Ones(5)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_eigen_matrix_arg) {
+  stan::math::test::test_slices(
+      5 * (10 + 50), 10.0,
+      std::vector<Eigen::MatrixXd>(2, Eigen::MatrixXd::Ones(5, 5)));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_std_vector_double_arg) {
+  stan::math::test::test_slices(5 * (10 + 1250), 10.0,
+                                std::vector<std::vector<std::vector<double>>>(
+                                    5, std::vector<std::vector<double>>(
+                                           5, std::vector<double>(5, 10.0))));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_eigen_vector_arg) {
+  stan::math::test::test_slices(
+      5 * (10 + 20), 10.0,
+      std::vector<std::vector<Eigen::VectorXd>>(
+          2, std::vector<Eigen::VectorXd>(2, Eigen::VectorXd::Ones(5))));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_eigen_row_vector_arg) {
+  stan::math::test::test_slices(
+      5 * (10 + 20), 10.0,
+      std::vector<std::vector<Eigen::RowVectorXd>>(
+          2, std::vector<Eigen::RowVectorXd>(2, Eigen::RowVectorXd::Ones(5))));
+}
+
+TEST(StanMathPrim_reduce_sum, std_vector_std_vector_eigen_matrix_arg) {
+  stan::math::test::test_slices(
+      5 * (10 + 60), 10.0,
+      std::vector<std::vector<Eigen::MatrixXd>>(
+          2, std::vector<Eigen::MatrixXd>(2, Eigen::MatrixXd::Ones(5, 3))));
+}
+
+TEST(StanMathPrim_reduce_sum, sum) {
+  double answer = 5 + 5 * (1 + 1 + 5 + 5 + 5 + 5 + 25 + 10 + 10);
+  stan::math::test::test_slices(
+      answer, 1.0, 1, 1.0, std::vector<int>(5, 1), std::vector<double>(5, 1.0),
+      Eigen::VectorXd::Ones(5), Eigen::RowVectorXd::Ones(5),
+      Eigen::MatrixXd::Ones(5, 5),
+      std::vector<std::vector<double>>(2, std::vector<double>(5, 1.0)),
+      std::vector<Eigen::VectorXd>(2, Eigen::VectorXd::Ones(5)));
+}
+
+#ifdef STAN_THREADS
+std::vector<int> threading_test_global;
+struct threading_test_lpdf {
+  template <typename T1>
+  inline auto operator()(std::size_t start, std::size_t end,
+                         const std::vector<T1>&, std::ostream* msgs) const {
+    threading_test_global[start] = tbb::this_task_arena::current_thread_index();
+
+    return stan::return_type_t<T1>(0);
+  }
+};
+
+TEST(StanMathPrim_reduce_sum, threading) {
+  threading_test_global = std::vector<int>(10000, 0);
+  stan::math::reduce_sum_static<threading_test_lpdf>(threading_test_global, 1,
+                                                     nullptr);
+
+  auto uniques = std::set<int>(threading_test_global.begin(),
+                               threading_test_global.end());
+
+  EXPECT_GT(uniques.size(), 1);
+
+  threading_test_global = std::vector<int>(10000, 0);
+
+  stan::math::reduce_sum<threading_test_lpdf>(threading_test_global, 1,
+                                              nullptr);
+
+  uniques = std::set<int>(threading_test_global.begin(),
+                          threading_test_global.end());
+
+  EXPECT_GT(uniques.size(), 1);
+}
+#endif
diff --git a/test/unit/math/prim/functor/reduce_sum_util.hpp b/test/unit/math/prim/functor/reduce_sum_util.hpp
new file mode 100644
index 00000000000..231f7e32d8b
--- /dev/null
+++ b/test/unit/math/prim/functor/reduce_sum_util.hpp
@@ -0,0 +1,186 @@
+#ifndef TEST_UNIT_MATH_PRIM_FUNCTOR_REDUCE_SUM_UTIL
+#define TEST_UNIT_MATH_PRIM_FUNCTOR_REDUCE_SUM_UTIL
+
+#include <stan/math.hpp>
+#include <gtest/gtest.h>
+#include <algorithm>
+#include <sstream>
+#include <tuple>
+#include <vector>
+
+namespace stan {
+namespace math {
+namespace test {
+std::ostream* get_new_msg() {
+  std::ostream* msgs = nullptr;
+  return msgs;
+}
+
+template <typename T>
+struct count_lpdf {
+  count_lpdf() {}
+
+  // does the reduction in the sub-slice start to end
+  inline T operator()(std::size_t start, std::size_t end,
+                      const std::vector<int>& sub_slice, std::ostream* msgs,
+                      const std::vector<T>& lambda,
+                      const std::vector<int>& idata) const {
+    return stan::math::poisson_lpmf(sub_slice, lambda[0]);
+  }
+};
+
+template <typename T>
+struct nesting_count_lpdf {
+  nesting_count_lpdf() {}
+
+  // does the reduction in the sub-slice start to end
+  inline T operator()(std::size_t start, std::size_t end,
+                      const std::vector<int>& sub_slice, std::ostream* msgs,
+                      const std::vector<T>& lambda,
+                      const std::vector<int>& idata) const {
+    return stan::math::reduce_sum<count_lpdf<T>>(sub_slice, 5, msgs, lambda,
+                                                 idata);
+  }
+};
+
+struct sum_lpdf {
+  template <typename T, stan::require_stan_scalar_t<T>* = nullptr>
+  T sum_(T arg) const {
+    return arg;
+  }
+
+  template <typename EigMat, stan::require_eigen_t<EigMat>* = nullptr>
+  auto sum_(EigMat&& arg) const {
+    return stan::math::sum(arg);
+  }
+
+  template <typename Vec, stan::require_std_vector_t<Vec>* = nullptr>
+  auto sum_(Vec&& arg) const {
+    stan::scalar_type_t<Vec> sum = 0;
+    for (size_t i = 0; i < arg.size(); ++i) {
+      sum += sum_(arg[i]);
+    }
+    return sum;
+  }
+
+  template <typename T, typename... Args>
+  inline auto operator()(std::size_t start, std::size_t end, T&& sub_slice,
+                         std::ostream* msgs, Args&&... args) const {
+    using return_type = stan::return_type_t<T, Args...>;
+
+    return sum_(sub_slice)
+           + sub_slice.size()
+                 * stan::math::sum(std::vector<return_type>{
+                       return_type(sum_(std::forward<Args>(args)))...});
+  }
+};
+
+struct start_end_lpdf {
+  template <typename T1, typename T2>
+  inline auto operator()(std::size_t start, std::size_t end, T1&&,
+                         std::ostream* msgs, T2&& data) const {
+    stan::return_type_t<T1, T2> sum = 0;
+    EXPECT_GE(start, 0);
+    EXPECT_LE(end, data.size() - 1);
+    for (size_t i = start; i <= end; i++) {
+      sum += data[i];
+    }
+    return sum;
+  }
+};
+
+/**
+ * slice over the grouping variable which is a var
+ */
+template <typename T>
+struct slice_group_count_lpdf {
+  slice_group_count_lpdf() {}
+
+  // does the reduction in the sub-slice start to end
+  inline T operator()(std::size_t start, std::size_t end,
+                      const std::vector<T>& lambda_slice, std::ostream* msgs,
+                      const std::vector<int>& y,
+                      const std::vector<int>& gsidx) const {
+    const std::size_t num_groups = end - start + 1;
+    T result = 0.0;
+    for (std::size_t i = 0; i != num_groups; ++i) {
+      std::vector<int> y_group(y.begin() + gsidx[start + i],
+                               y.begin() + gsidx[start + i + 1]);
+      result += stan::math::poisson_lpmf(y_group, lambda_slice[i]);
+    }
+    return result;
+  }
+};
+
+/**
+ * basic performance test for a hierarchical model
+ */
+template <typename T>
+struct grouped_count_lpdf {
+  grouped_count_lpdf() {}
+
+  // does the reduction in the sub-slice start to end
+  template <typename VecInt1, typename VecT, typename VecInt2>
+  inline T operator()(std::size_t start, std::size_t end, VecInt1&& sub_slice,
+                      std::ostream* msgs, VecT&& lambda, VecInt2&& gidx) const {
+    const std::size_t num_terms = end - start + 1;
+    std::decay_t<VecT> lambda_slice(num_terms);
+    for (std::size_t i = 0; i != num_terms; ++i)
+      lambda_slice[i] = lambda[gidx[start + i]];
+
+    return stan::math::poisson_lpmf(sub_slice, lambda_slice);
+  }
+};
+
+template <typename T1, typename T2, typename... Args>
+void test_slices(T1 result, T2&& vec_value, Args&&... args) {
+  using stan::math::test::get_new_msg;
+  using stan::math::test::sum_lpdf;
+
+  std::vector<std::decay_t<T2>> data(5, vec_value);
+  EXPECT_EQ(result, stan::math::reduce_sum_static<sum_lpdf>(
+                        data, 1, get_new_msg(), args...))
+      << "Failed for reduce_sum_static";
+  EXPECT_EQ(result,
+            stan::math::reduce_sum<sum_lpdf>(data, 1, get_new_msg(), args...))
+      << "Failed for reduce_sum";
+}
+
+auto reduce_sum_static_int_sum_lpdf = [](auto&&... args) {
+  return stan::math::reduce_sum_static<sum_lpdf>(std::vector<int>(2, 10.0), 1,
+                                                 get_new_msg(), args...);
+};
+
+auto reduce_sum_static_sum_lpdf = [](auto&& data, auto&&... args) {
+  return stan::math::reduce_sum_static<sum_lpdf>(data, 1, get_new_msg(),
+                                                 args...);
+};
+
+auto reduce_sum_int_sum_lpdf = [](auto&&... args) {
+  return stan::math::reduce_sum<sum_lpdf>(std::vector<int>(2, 10.0), 1,
+                                          get_new_msg(), args...);
+};
+
+auto reduce_sum_sum_lpdf = [](auto&& data, auto&&... args) {
+  return stan::math::reduce_sum<sum_lpdf>(data, 1, get_new_msg(), args...);
+};
+
+template <int grainsize>
+struct static_check_lpdf {
+  template <typename T>
+  inline auto operator()(std::size_t start, std::size_t end,
+                         const std::vector<int>&, std::ostream* msgs,
+                         const std::vector<T>& data) const {
+    T sum = 0;
+    EXPECT_LE(end - start + 1, grainsize);
+    for (size_t i = start; i <= end; i++) {
+      sum += data[i];
+    }
+    return sum;
+  }
+};
+
+}  // namespace test
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/unit/math/rev/functor/reduce_sum_test.cpp b/test/unit/math/rev/functor/reduce_sum_test.cpp
new file mode 100644
index 00000000000..b2f0212489e
--- /dev/null
+++ b/test/unit/math/rev/functor/reduce_sum_test.cpp
@@ -0,0 +1,426 @@
+#include <stan/math.hpp>
+#include <test/unit/math/prim/functor/reduce_sum_util.hpp>
+#include <gtest/gtest.h>
+#include <algorithm>
+#include <sstream>
+#include <tuple>
+#include <vector>
+#include <set>
+
+TEST(StanMathRev_reduce_sum, no_args) {
+  using stan::math::var;
+  using stan::math::test::get_new_msg;
+  using stan::math::test::sum_lpdf;
+  std::vector<var> data(0);
+  EXPECT_EQ(0.0, stan::math::reduce_sum_static<sum_lpdf>(
+                     data, 1, stan::math::test::get_new_msg())
+                     .val())
+      << "Failed for reduce_sum_static";
+  EXPECT_EQ(0.0, stan::math::reduce_sum<sum_lpdf>(
+                     data, 1, stan::math::test::get_new_msg())
+                     .val())
+      << "Failed for reduce_sum";
+}
+
+TEST(StanMathRev_reduce_sum, value) {
+  using stan::math::test::count_lpdf;
+  using stan::math::test::get_new_msg;
+  double lambda_d = 10.0;
+  const std::size_t elems = 10000;
+  std::vector<int> data(elems);
+
+  for (std::size_t i = 0; i != elems; ++i)
+    data[i] = i;
+
+  std::vector<int> idata;
+  std::vector<double> vlambda_d(1, lambda_d);
+
+  double poisson_lpdf = stan::math::reduce_sum<count_lpdf<double>>(
+      data, 5, get_new_msg(), vlambda_d, idata);
+
+  double poisson_lpdf_ref = stan::math::poisson_lpmf(data, lambda_d);
+
+  EXPECT_FLOAT_EQ(poisson_lpdf, poisson_lpdf_ref)
+      << "ref value of poisson lpdf : " << poisson_lpdf_ref << std::endl
+      << "value of poisson lpdf : " << poisson_lpdf << std::endl;
+
+  double poisson_lpdf_static
+      = stan::math::reduce_sum_static<count_lpdf<double>>(
+          data, 5, get_new_msg(), vlambda_d, idata);
+
+  EXPECT_FLOAT_EQ(poisson_lpdf_static, poisson_lpdf_ref);
+}
+
+TEST(StanMathRev_reduce_sum, gradient) {
+  using stan::math::var;
+  using stan::math::test::count_lpdf;
+  using stan::math::test::get_new_msg;
+
+  double lambda_d = 10.0;
+  const std::size_t elems = 10000;
+  std::vector<int> data(elems);
+
+  for (std::size_t i = 0; i != elems; ++i)
+    data[i] = i;
+
+  var lambda_v = lambda_d;
+
+  std::vector<int> idata;
+  std::vector<var> vlambda_v(1, lambda_v);
+
+  var poisson_lpdf = stan::math::reduce_sum<count_lpdf<var>>(
+      data, 5, get_new_msg(), vlambda_v, idata);
+
+  var lambda_ref = lambda_d;
+  var poisson_lpdf_ref = stan::math::poisson_lpmf(data, lambda_ref);
+
+  EXPECT_FLOAT_EQ(value_of(poisson_lpdf), value_of(poisson_lpdf_ref));
+
+  stan::math::grad(poisson_lpdf_ref.vi_);
+  const double lambda_ref_adj = lambda_ref.adj();
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf.vi_);
+  const double lambda_adj = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj, lambda_ref_adj)
+      << "ref value of poisson lpdf : " << poisson_lpdf_ref.val() << std::endl
+      << "ref gradient wrt to lambda: " << lambda_ref_adj << std::endl
+      << "value of poisson lpdf : " << poisson_lpdf.val() << std::endl
+      << "gradient wrt to lambda: " << lambda_adj << std::endl;
+
+  var poisson_lpdf_static = stan::math::reduce_sum_static<count_lpdf<var>>(
+      data, 5, get_new_msg(), vlambda_v, idata);
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf_static.vi_);
+  const double lambda_adj_static = lambda_v.adj();
+  EXPECT_FLOAT_EQ(lambda_adj_static, lambda_ref_adj);
+  stan::math::recover_memory();
+}
+
+TEST(StanMathRev_reduce_sum, grainsize) {
+  using stan::math::var;
+  using stan::math::test::count_lpdf;
+  using stan::math::test::get_new_msg;
+
+  double lambda_d = 10.0;
+  const std::size_t elems = 10000;
+  std::vector<int> data(elems);
+
+  for (std::size_t i = 0; i != elems; ++i)
+    data[i] = i;
+
+  var lambda_v = lambda_d;
+
+  std::vector<int> idata;
+  std::vector<var> vlambda_v(1, lambda_v);
+
+  EXPECT_THROW(stan::math::reduce_sum<count_lpdf<var>>(data, 0, get_new_msg(),
+                                                       vlambda_v, idata),
+               std::domain_error);
+
+  EXPECT_THROW(stan::math::reduce_sum<count_lpdf<var>>(data, -1, get_new_msg(),
+                                                       vlambda_v, idata),
+               std::domain_error);
+
+  EXPECT_NO_THROW(stan::math::reduce_sum<count_lpdf<var>>(
+      data, 1, get_new_msg(), vlambda_v, idata));
+
+  EXPECT_NO_THROW(stan::math::reduce_sum<count_lpdf<var>>(
+      data, 2 * elems, get_new_msg(), vlambda_v, idata));
+
+  EXPECT_THROW(stan::math::reduce_sum_static<count_lpdf<var>>(
+                   data, 0, get_new_msg(), vlambda_v, idata),
+               std::domain_error);
+
+  EXPECT_THROW(stan::math::reduce_sum_static<count_lpdf<var>>(
+                   data, -1, get_new_msg(), vlambda_v, idata),
+               std::domain_error);
+
+  EXPECT_NO_THROW(stan::math::reduce_sum_static<count_lpdf<var>>(
+      data, 1, get_new_msg(), vlambda_v, idata));
+
+  EXPECT_NO_THROW(stan::math::reduce_sum_static<count_lpdf<var>>(
+      data, 2 * elems, get_new_msg(), vlambda_v, idata));
+
+  stan::math::recover_memory();
+}
+
+TEST(StanMathRev_reduce_sum, nesting_gradient) {
+  using stan::math::var;
+  using stan::math::test::get_new_msg;
+  using stan::math::test::nesting_count_lpdf;
+
+  double lambda_d = 10.0;
+  const std::size_t elems = 10000;
+  std::vector<int> data(elems);
+
+  for (std::size_t i = 0; i != elems; ++i)
+    data[i] = i;
+
+  var lambda_v = lambda_d;
+
+  std::vector<int> idata;
+  std::vector<var> vlambda_v(1, lambda_v);
+
+  var poisson_lpdf = stan::math::reduce_sum<nesting_count_lpdf<var>>(
+      data, 5, get_new_msg(), vlambda_v, idata);
+
+  var lambda_ref = lambda_d;
+  var poisson_lpdf_ref = stan::math::poisson_lpmf(data, lambda_ref);
+
+  EXPECT_FLOAT_EQ(value_of(poisson_lpdf), value_of(poisson_lpdf_ref));
+
+  stan::math::grad(poisson_lpdf_ref.vi_);
+  const double lambda_ref_adj = lambda_ref.adj();
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf.vi_);
+  const double lambda_adj = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj, lambda_ref_adj)
+      << "ref value of poisson lpdf : " << poisson_lpdf_ref.val() << std::endl
+      << "ref gradient wrt to lambda: " << lambda_ref_adj << std::endl
+      << "value of poisson lpdf : " << poisson_lpdf.val() << std::endl
+      << "gradient wrt to lambda: " << lambda_adj << std::endl;
+
+  var poisson_lpdf_static
+      = stan::math::reduce_sum_static<nesting_count_lpdf<var>>(
+          data, 5, get_new_msg(), vlambda_v, idata);
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf_static.vi_);
+  const double lambda_adj_static = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj_static, lambda_ref_adj);
+  stan::math::recover_memory();
+}
+
+TEST(StanMathRev_reduce_sum, grouped_gradient) {
+  using stan::math::var;
+  using stan::math::test::get_new_msg;
+  using stan::math::test::grouped_count_lpdf;
+
+  double lambda_d = 10.0;
+  const std::size_t groups = 10;
+  const std::size_t elems_per_group = 1000;
+  const std::size_t elems = groups * elems_per_group;
+
+  std::vector<int> data(elems);
+  std::vector<int> gidx(elems);
+
+  for (std::size_t i = 0; i != elems; ++i) {
+    data[i] = i;
+    gidx[i] = i / elems_per_group;
+  }
+
+  std::vector<var> vlambda_v;
+
+  for (std::size_t i = 0; i != groups; ++i)
+    vlambda_v.push_back(i + 0.2);
+
+  var lambda_v = vlambda_v[0];
+
+  var poisson_lpdf = stan::math::reduce_sum<grouped_count_lpdf<var>>(
+      data, 5, get_new_msg(), vlambda_v, gidx);
+
+  std::vector<var> vref_lambda_v;
+  for (std::size_t i = 0; i != elems; ++i) {
+    vref_lambda_v.push_back(vlambda_v[gidx[i]]);
+  }
+  var lambda_ref = vlambda_v[0];
+  var poisson_lpdf_ref = stan::math::poisson_lpmf(data, vref_lambda_v);
+
+  EXPECT_FLOAT_EQ(value_of(poisson_lpdf), value_of(poisson_lpdf_ref));
+
+  stan::math::grad(poisson_lpdf_ref.vi_);
+  const double lambda_ref_adj = lambda_ref.adj();
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf.vi_);
+  const double lambda_adj = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj, lambda_ref_adj)
+      << "ref value of poisson lpdf : " << poisson_lpdf_ref.val() << std::endl
+      << "ref gradient wrt to lambda: " << lambda_ref_adj << std::endl
+      << "value of poisson lpdf : " << poisson_lpdf.val() << std::endl
+      << "gradient wrt to lambda: " << lambda_adj << std::endl;
+
+  var poisson_lpdf_static
+      = stan::math::reduce_sum_static<grouped_count_lpdf<var>>(
+          data, 5, get_new_msg(), vlambda_v, gidx);
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf_static.vi_);
+  const double lambda_adj_static = lambda_v.adj();
+  EXPECT_FLOAT_EQ(lambda_adj_static, lambda_ref_adj);
+  stan::math::recover_memory();
+}
+
+TEST(StanMathRev_reduce_sum, grouped_gradient_eigen) {
+  using stan::math::var;
+  using stan::math::test::get_new_msg;
+  using stan::math::test::grouped_count_lpdf;
+
+  double lambda_d = 10.0;
+  const std::size_t groups = 10;
+  const std::size_t elems_per_group = 1000;
+  const std::size_t elems = groups * elems_per_group;
+
+  std::vector<int> data(elems);
+  std::vector<int> gidx(elems);
+
+  for (std::size_t i = 0; i != elems; ++i) {
+    data[i] = i;
+    gidx[i] = i / elems_per_group;
+  }
+
+  Eigen::Matrix<var, -1, 1> vlambda_v(groups);
+
+  for (std::size_t i = 0; i != groups; ++i)
+    vlambda_v[i] = i + 0.2;
+  var lambda_v = vlambda_v[0];
+
+  var poisson_lpdf = stan::math::reduce_sum<grouped_count_lpdf<var>>(
+      data, 5, get_new_msg(), vlambda_v, gidx);
+
+  std::vector<var> vref_lambda_v;
+  for (std::size_t i = 0; i != elems; ++i) {
+    vref_lambda_v.push_back(vlambda_v[gidx[i]]);
+  }
+  var lambda_ref = vlambda_v[0];
+
+  var poisson_lpdf_ref = stan::math::poisson_lpmf(data, vref_lambda_v);
+
+  EXPECT_FLOAT_EQ(value_of(poisson_lpdf), value_of(poisson_lpdf_ref));
+
+  stan::math::grad(poisson_lpdf_ref.vi_);
+  const double lambda_ref_adj = lambda_ref.adj();
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf.vi_);
+  const double lambda_adj = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj, lambda_ref_adj)
+      << "ref value of poisson lpdf : " << poisson_lpdf_ref.val() << std::endl
+      << "ref gradient wrt to lambda: " << lambda_ref_adj << std::endl
+      << "value of poisson lpdf : " << poisson_lpdf.val() << std::endl
+      << "gradient wrt to lambda: " << lambda_adj << std::endl;
+
+  var poisson_lpdf_static
+      = stan::math::reduce_sum_static<grouped_count_lpdf<var>>(
+          data, 5, get_new_msg(), vlambda_v, gidx);
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf_static.vi_);
+  const double lambda_adj_static = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj_static, lambda_ref_adj);
+
+  stan::math::recover_memory();
+}
+
+TEST(StanMathRev_reduce_sum, slice_group_gradient) {
+  using stan::math::var;
+  using stan::math::test::get_new_msg;
+  using stan::math::test::slice_group_count_lpdf;
+
+  double lambda_d = 10.0;
+  const std::size_t groups = 10;
+  const std::size_t elems_per_group = 1000;
+  const std::size_t elems = groups * elems_per_group;
+
+  std::vector<int> data(elems);
+  std::vector<int> gidx(elems);
+  std::vector<int> gsidx(groups + 1);
+
+  for (std::size_t i = 0, k = 0; i != groups; ++i) {
+    gsidx[i] = k;
+    for (std::size_t j = 0; j != elems_per_group; ++j, ++k) {
+      data[k] = k;
+      gidx[k] = i;
+    }
+    gsidx[i + 1] = k;
+  }
+
+  std::vector<var> vlambda_v;
+
+  for (std::size_t i = 0; i != groups; ++i)
+    vlambda_v.push_back(i + 0.2);
+
+  var lambda_v = vlambda_v[0];
+
+  var poisson_lpdf = stan::math::reduce_sum<slice_group_count_lpdf<var>>(
+      vlambda_v, 5, get_new_msg(), data, gsidx);
+
+  std::vector<var> vref_lambda_v;
+  for (std::size_t i = 0; i != elems; ++i) {
+    vref_lambda_v.push_back(vlambda_v[gidx[i]]);
+  }
+  var lambda_ref = vlambda_v[0];
+
+  var poisson_lpdf_ref = stan::math::poisson_lpmf(data, vref_lambda_v);
+
+  EXPECT_FLOAT_EQ(value_of(poisson_lpdf), value_of(poisson_lpdf_ref));
+
+  stan::math::grad(poisson_lpdf_ref.vi_);
+  const double lambda_ref_adj = lambda_ref.adj();
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf.vi_);
+  const double lambda_adj = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj, lambda_ref_adj)
+      << "ref value of poisson lpdf : " << poisson_lpdf_ref.val() << std::endl
+      << "ref gradient wrt to lambda: " << lambda_ref_adj << std::endl
+      << "value of poisson lpdf : " << poisson_lpdf.val() << std::endl
+      << "gradient wrt to lambda: " << lambda_adj << std::endl;
+
+  var poisson_lpdf_static
+      = stan::math::reduce_sum_static<slice_group_count_lpdf<var>>(
+          vlambda_v, 5, get_new_msg(), data, gsidx);
+
+  stan::math::set_zero_all_adjoints();
+  stan::math::grad(poisson_lpdf_static.vi_);
+  const double lambda_adj_static = lambda_v.adj();
+
+  EXPECT_FLOAT_EQ(lambda_adj_static, lambda_ref_adj);
+
+  stan::math::recover_memory();
+}
+
+#ifdef STAN_THREADS
+std::vector<int> threading_test_global;
+struct threading_test_lpdf {
+  template <typename T1>
+  inline auto operator()(std::size_t start, std::size_t end,
+                         const std::vector<T1>&, std::ostream* msgs) const {
+    threading_test_global[start] = tbb::this_task_arena::current_thread_index();
+
+    return stan::return_type_t<T1>(0);
+  }
+};
+
+TEST(StanMathRev_reduce_sum, threading) {
+  threading_test_global = std::vector<int>(10000, 0);
+  std::vector<stan::math::var> data(threading_test_global.size(), 0);
+  stan::math::reduce_sum_static<threading_test_lpdf>(data, 1, nullptr);
+
+  auto uniques = std::set<int>(threading_test_global.begin(),
+                               threading_test_global.end());
+
+  EXPECT_GT(uniques.size(), 1);
+
+  threading_test_global = std::vector<int>(10000, 0);
+
+  stan::math::reduce_sum<threading_test_lpdf>(data, 1, nullptr);
+
+  uniques = std::set<int>(threading_test_global.begin(),
+                          threading_test_global.end());
+
+  EXPECT_GT(uniques.size(), 1);
+
+  stan::math::recover_memory();
+}
+#endif