diff --git a/stan/math/rev/core.hpp b/stan/math/rev/core.hpp
index 5d3c91f423b..33eb934c08e 100644
--- a/stan/math/rev/core.hpp
+++ b/stan/math/rev/core.hpp
@@ -1,13 +1,16 @@
 #ifndef STAN_MATH_REV_CORE_HPP
 #define STAN_MATH_REV_CORE_HPP
 
+#include <stan/math/rev/core/accumulate_adjoints.hpp>
 #include <stan/math/rev/core/autodiffstackstorage.hpp>
 #include <stan/math/rev/core/build_vari_array.hpp>
 #include <stan/math/rev/core/chainable_alloc.hpp>
 #include <stan/math/rev/core/chainablestack.hpp>
+#include <stan/math/rev/core/count_vars.hpp>
 #include <stan/math/rev/core/init_chainablestack.hpp>
 #include <stan/math/rev/core/std_iterator_traits.hpp>
 #include <stan/math/rev/core/ddv_vari.hpp>
+#include <stan/math/rev/core/deep_copy_vars.hpp>
 #include <stan/math/rev/core/dv_vari.hpp>
 #include <stan/math/rev/core/dvd_vari.hpp>
 #include <stan/math/rev/core/dvv_vari.hpp>
@@ -63,5 +66,6 @@
 #include <stan/math/rev/core/vv_vari.hpp>
 #include <stan/math/rev/core/vvd_vari.hpp>
 #include <stan/math/rev/core/vvv_vari.hpp>
+#include <stan/math/rev/core/save_varis.hpp>
 
 #endif
diff --git a/stan/math/rev/core/accumulate_adjoints.hpp b/stan/math/rev/core/accumulate_adjoints.hpp
new file mode 100644
index 00000000000..98f0899119d
--- /dev/null
+++ b/stan/math/rev/core/accumulate_adjoints.hpp
@@ -0,0 +1,154 @@
+#ifndef STAN_MATH_REV_CORE_ACCUMULATE_ADJOINTS_HPP
+#define STAN_MATH_REV_CORE_ACCUMULATE_ADJOINTS_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core/var.hpp>
+
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+template <typename... Pargs>
+inline double* accumulate_adjoints(double* dest, const var& x, Pargs&&... args);
+
+template <typename VarVec, require_std_vector_vt<is_var, VarVec>* = nullptr,
+          typename... Pargs>
+inline double* accumulate_adjoints(double* dest, VarVec&& x, Pargs&&... args);
+
+template <typename VecContainer,
+          require_std_vector_st<is_var, VecContainer>* = nullptr,
+          require_std_vector_vt<is_container, VecContainer>* = nullptr,
+          typename... Pargs>
+inline double* accumulate_adjoints(double* dest, VecContainer&& x,
+                                   Pargs&&... args);
+
+template <typename EigT, require_eigen_vt<is_var, EigT>* = nullptr,
+          typename... Pargs>
+inline double* accumulate_adjoints(double* dest, EigT&& x, Pargs&&... args);
+
+template <typename Arith, require_arithmetic_t<scalar_type_t<Arith>>* = nullptr,
+          typename... Pargs>
+inline double* accumulate_adjoints(double* dest, Arith&& x, Pargs&&... args);
+
+inline double* accumulate_adjoints(double* dest);
+
+/**
+ * Accumulate adjoints from x into storage pointed to by dest,
+ *   increment the adjoint storage pointer,
+ *   recursively accumulate the adjoints of the rest of the arguments,
+ *   and return final position of storage pointer.
+ *
+ * @tparam Pargs Types of remaining arguments
+ * @param dest Pointer to where adjoints are to be accumulated
+ * @param x A var
+ * @param args Further args to accumulate over
+ * @return Final position of adjoint storage pointer
+ */
+template <typename... Pargs>
+inline double* accumulate_adjoints(double* dest, const var& x,
+                                   Pargs&&... args) {
+  *dest += x.adj();
+  return accumulate_adjoints(dest + 1, std::forward<Pargs>(args)...);
+}
+
+/**
+ * Accumulate adjoints from std::vector x into storage pointed to by dest,
+ *   increment the adjoint storage pointer,
+ *   recursively accumulate the adjoints of the rest of the arguments,
+ *   and return final position of storage pointer.
+ *
+ * @tparam Pargs Types of remaining arguments
+ * @param dest Pointer to where adjoints are to be accumulated
+ * @param x A std::vector of vars
+ * @param args Further args to accumulate over
+ * @return Final position of adjoint storage pointer
+ */
+template <typename VarVec, require_std_vector_vt<is_var, VarVec>*,
+          typename... Pargs>
+inline double* accumulate_adjoints(double* dest, VarVec&& x, Pargs&&... args) {
+  for (auto&& x_iter : x) {
+    *dest += x_iter.adj();
+    ++dest;
+  }
+  return accumulate_adjoints(dest, std::forward<Pargs>(args)...);
+}
+
+/**
+ * Accumulate adjoints from x (a std::vector of containers containing vars)
+ *   into storage pointed to by dest,
+ *   increment the adjoint storage pointer,
+ *   recursively accumulate the adjoints of the rest of the arguments,
+ *   and return final position of storage pointer.
+ *
+ * @tparam VecContainer the type of a standard container holding var
+ * containers.
+ * @tparam Pargs Types of remaining arguments
+ * @param dest Pointer to where adjoints are to be accumulated
+ * @param x A std::vector of containers holding vars
+ * @param args Further args to accumulate over
+ * @return Final position of adjoint storage pointer
+ */
+template <typename VecContainer, require_std_vector_st<is_var, VecContainer>*,
+          require_std_vector_vt<is_container, VecContainer>*, typename... Pargs>
+inline double* accumulate_adjoints(double* dest, VecContainer&& x,
+                                   Pargs&&... args) {
+  for (auto&& x_iter : x) {
+    dest = accumulate_adjoints(dest, x_iter);
+  }
+  return accumulate_adjoints(dest, std::forward<Pargs>(args)...);
+}
+
+/**
+ * Accumulate adjoints from x (an Eigen type containing vars)
+ *   into storage pointed to by dest,
+ *   increment the adjoint storage pointer,
+ *   recursively accumulate the adjoints of the rest of the arguments,
+ *   and return final position of storage pointer.
+ *
+ * @tparam EigT Type derived from `EigenBase` containing vars.
+ * @tparam Pargs Types of remaining arguments
+ * @param dest Pointer to where adjoints are to be accumulated
+ * @param x An eigen type holding vars to accumulate over
+ * @param args Further args to accumulate over
+ * @return Final position of adjoint storage pointer
+ */
+template <typename EigT, require_eigen_vt<is_var, EigT>*, typename... Pargs>
+inline double* accumulate_adjoints(double* dest, EigT&& x, Pargs&&... args) {
+  Eigen::Map<Eigen::MatrixXd>(dest, x.rows(), x.cols()) += x.adj();
+  return accumulate_adjoints(dest + x.size(), std::forward<Pargs>(args)...);
+}
+
+/**
+ * Ignore arithmetic types.
+ *
+ * Recursively accumulate the adjoints of the rest of the arguments
+ *   and return final position of adjoint storage pointer.
+ *
+ * @tparam Arith A type satisfying `std::is_arithmetic`.
+ * @tparam Pargs Types of remaining arguments
+ * @param dest Pointer to where adjoints are to be accumulated
+ * @param x An object that is either arithmetic or a container of Arithmetic
+ *  types
+ * @param args Further args to accumulate over
+ * @return Final position of adjoint storage pointer
+ */
+template <typename Arith, require_arithmetic_t<scalar_type_t<Arith>>*,
+          typename... Pargs>
+inline double* accumulate_adjoints(double* dest, Arith&& x, Pargs&&... args) {
+  return accumulate_adjoints(dest, std::forward<Pargs>(args)...);
+}
+
+/**
+ * End accumulate_adjoints recursion and return pointer
+ *
+ * @param dest Pointer
+ */
+inline double* accumulate_adjoints(double* dest) { return dest; }
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/core/count_vars.hpp b/stan/math/rev/core/count_vars.hpp
new file mode 100644
index 00000000000..b0b536a27ab
--- /dev/null
+++ b/stan/math/rev/core/count_vars.hpp
@@ -0,0 +1,153 @@
+#ifndef STAN_MATH_REV_CORE_COUNT_VARS_HPP
+#define STAN_MATH_REV_CORE_COUNT_VARS_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core/var.hpp>
+
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+template <typename VecVar, require_std_vector_vt<is_var, VecVar>* = nullptr,
+          typename... Pargs>
+inline size_t count_vars_impl(size_t count, VecVar&& x, Pargs&&... args);
+
+template <typename VecContainer,
+          require_std_vector_st<is_var, VecContainer>* = nullptr,
+          require_std_vector_vt<is_container, VecContainer>* = nullptr,
+          typename... Pargs>
+inline size_t count_vars_impl(size_t count, VecContainer&& x, Pargs&&... args);
+
+template <typename EigT, require_eigen_vt<is_var, EigT>* = nullptr,
+          typename... Pargs>
+inline size_t count_vars_impl(size_t count, EigT&& x, Pargs&&... args);
+
+template <typename... Pargs>
+inline size_t count_vars_impl(size_t count, const var& x, Pargs&&... args);
+
+template <typename Arith, require_arithmetic_t<scalar_type_t<Arith>>* = nullptr,
+          typename... Pargs>
+inline size_t count_vars_impl(size_t count, Arith& x, Pargs&&... args);
+
+inline size_t count_vars_impl(size_t count);
+
+/**
+ * Count the number of vars in x (a std::vector of vars),
+ *  add it to the running total,
+ *  count the number of vars in the remaining arguments
+ *  and return the result.
+ *
+ * @tparam VecVar type of standard container holding vars
+ * @tparam Pargs Types of remaining arguments
+ * @param[in] count The current count of the number of vars
+ * @param[in] x A std::vector holding vars.
+ * @param[in] args objects to be forwarded to recursive call of
+ * `count_vars_impl`
+ */
+template <typename VecVar, require_std_vector_vt<is_var, VecVar>*,
+          typename... Pargs>
+inline size_t count_vars_impl(size_t count, VecVar&& x, Pargs&&... args) {
+  return count_vars_impl(count + x.size(), std::forward<Pargs>(args)...);
+}
+
+/**
+ * Count the number of vars in x (a std::vector holding other containers),
+ *  add it to the running total,
+ *  count the number of vars in the remaining arguments
+ *  and return the result.
+ *
+ * @tparam VecContainer std::vector holding arguments which contain Vars
+ * @tparam Pargs Types of remaining arguments
+ * @param[in] count The current count of the number of vars
+ * @param[in] x A vector holding containers of vars
+ * @param[in] args objects to be forwarded to recursive call of
+ * `count_vars_impl`
+ */
+template <typename VecContainer, require_std_vector_st<is_var, VecContainer>*,
+          require_std_vector_vt<is_container, VecContainer>*, typename... Pargs>
+inline size_t count_vars_impl(size_t count, VecContainer&& x, Pargs&&... args) {
+  for (auto&& x_iter : x) {
+    count = count_vars_impl(count, x_iter);
+  }
+  return count_vars_impl(count, std::forward<Pargs>(args)...);
+}
+
+/**
+ * Count the number of vars in x (an eigen container),
+ *  add it to the running total,
+ *  count the number of vars in the remaining arguments
+ *  and return the result.
+ *
+ * @tparam EigT A type derived from `EigenBase`
+ * @tparam Pargs Types of remaining arguments
+ * @param[in] count The current count of the number of vars
+ * @param[in] x An Eigen container holding vars
+ * @param[in] args objects to be forwarded to recursive call of
+ * `count_vars_impl`
+ */
+template <typename EigT, require_eigen_vt<is_var, EigT>*, typename... Pargs>
+inline size_t count_vars_impl(size_t count, EigT&& x, Pargs&&... args) {
+  return count_vars_impl(count + x.size(), std::forward<Pargs>(args)...);
+}
+
+/**
+ * Add one to the running total number of vars,
+ *  count the number of vars in the remaining arguments
+ *  and return the result.
+ *
+ * @tparam Pargs Types of remaining arguments
+ * @param[in] count The current count of the number of vars
+ * @param[in] x A var
+ * @param[in] args objects to be forwarded to recursive call of
+ * `count_vars_impl`
+ */
+template <typename... Pargs>
+inline size_t count_vars_impl(size_t count, const var& x, Pargs&&... args) {
+  return count_vars_impl(count + 1, std::forward<Pargs>(args)...);
+}
+
+/**
+ * Arguments without vars contribute zero to the total number of vars.
+ *
+ * Return the running total number of vars plus the number of
+ *   vars in the remaining aruments.
+ *
+ * @tparam Arith An object that is either arithmetic or holds arithmetic
+ * types
+ * @tparam Pargs Types of remaining arguments
+ * @param[in] count The current count of the number of vars
+ * @param[in] x An arithmetic value or container
+ * @param[in] args objects to be forwarded to recursive call of
+ * `count_vars_impl`
+ */
+template <typename Arith, require_arithmetic_t<scalar_type_t<Arith>>*,
+          typename... Pargs>
+inline size_t count_vars_impl(size_t count, Arith& x, Pargs&&... args) {
+  return count_vars_impl(count, std::forward<Pargs>(args)...);
+}
+
+/**
+ * End count_vars_impl recursion and return total number of counted vars
+ */
+inline size_t count_vars_impl(size_t count) { return count; }
+}  // namespace internal
+
+/**
+ * Count the number of vars in the input argument list
+ *
+ * @tparam Pargs Types of input arguments
+ * @return Number of vars in input
+ */
+template <typename... Pargs>
+inline size_t count_vars(Pargs&&... args) {
+  return internal::count_vars_impl(0, std::forward<Pargs>(args)...);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/core/deep_copy_vars.hpp b/stan/math/rev/core/deep_copy_vars.hpp
new file mode 100644
index 00000000000..5eb1b3d8542
--- /dev/null
+++ b/stan/math/rev/core/deep_copy_vars.hpp
@@ -0,0 +1,87 @@
+#ifndef STAN_MATH_REV_CORE_DEEP_COPY_VARS_HPP
+#define STAN_MATH_REV_CORE_DEEP_COPY_VARS_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core/var.hpp>
+
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Forward arguments that do not contain vars. There
+ *   is no copying to be done.
+ *
+ * @tparam Arith an arithmetic type.
+ * @param arg For lvalue references this will be passed by reference.
+ *  Otherwise it will be moved.
+ */
+template <typename Arith, typename = require_arithmetic_t<scalar_type_t<Arith>>>
+inline decltype(auto) deep_copy_vars(Arith&& arg) {
+  return std::forward<Arith>(arg);
+}
+
+/**
+ * Copy the value of a var but reallocate a new vari
+ *
+ * @param arg A var
+ * @return A new var
+ */
+inline auto deep_copy_vars(const var& arg) {
+  return var(new vari(arg.val(), false));
+}
+
+/**
+ * Copy the vars in arg but reallocate new varis for them
+ *
+ * @tparam VarVec A variant of std::vector<var>
+ * @param arg A std::vector of vars
+ * @return A new std::vector of vars
+ */
+template <typename VarVec, require_std_vector_vt<is_var, VarVec>* = nullptr>
+inline auto deep_copy_vars(VarVec&& arg) {
+  std::vector<var> copy_vec(arg.size());
+  for (size_t i = 0; i < arg.size(); ++i) {
+    copy_vec[i] = new vari(arg[i].val(), false);
+  }
+  return copy_vec;
+}
+
+/**
+ * Copy the vars in arg but reallocate new varis for them
+ *
+ * @tparam VecContainer std::vector<T> where T is another type containing vars
+ * @param arg A std::vector of containers containing vars
+ * @return A new std::vector of containers containing vars
+ */
+template <typename VecContainer,
+          require_std_vector_st<is_var, VecContainer>* = nullptr,
+          require_std_vector_vt<is_container, VecContainer>* = nullptr>
+inline auto deep_copy_vars(VecContainer&& arg) {
+  std::vector<value_type_t<VecContainer>> copy_vec(arg.size());
+  for (size_t i = 0; i < arg.size(); ++i) {
+    copy_vec[i] = deep_copy_vars(arg[i]);
+  }
+  return copy_vec;
+}
+
+/**
+ * Copy the vars in arg but reallocate new varis for them
+ *
+ * @tparam EigT An Eigen type with var value type
+ * @param arg An Eigen container of vars
+ * @return A new Eigen container of vars
+ */
+template <typename EigT, require_eigen_vt<is_var, EigT>* = nullptr>
+inline auto deep_copy_vars(EigT&& arg) {
+  return arg.unaryExpr([](auto&& x) { return var(new vari(x.val(), false)); })
+      .eval();
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/core/save_varis.hpp b/stan/math/rev/core/save_varis.hpp
new file mode 100644
index 00000000000..d97e8a1a474
--- /dev/null
+++ b/stan/math/rev/core/save_varis.hpp
@@ -0,0 +1,150 @@
+#ifndef STAN_MATH_REV_CORE_SAVE_VARIS_HPP
+#define STAN_MATH_REV_CORE_SAVE_VARIS_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core/var.hpp>
+
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+template <typename... Pargs>
+inline vari** save_varis(vari** dest, const var& x, Pargs&&... args);
+
+template <typename VarVec, require_std_vector_vt<is_var, VarVec>* = nullptr,
+          typename... Pargs>
+inline vari** save_varis(vari** dest, VarVec&& x, Pargs&&... args);
+
+template <typename VecContainer,
+          require_std_vector_st<is_var, VecContainer>* = nullptr,
+          require_std_vector_vt<is_container, VecContainer>* = nullptr,
+          typename... Pargs>
+inline vari** save_varis(vari** dest, VecContainer&& x, Pargs&&... args);
+
+template <typename EigT, require_eigen_vt<is_var, EigT>* = nullptr,
+          typename... Pargs>
+inline vari** save_varis(vari** dest, EigT&& x, Pargs&&... args);
+
+template <typename Arith, require_arithmetic_t<scalar_type_t<Arith>>* = nullptr,
+          typename... Pargs>
+inline vari** save_varis(vari** dest, Arith&& x, Pargs&&... args);
+
+inline vari** save_varis(vari** dest);
+
+/**
+ * Save the vari pointer in x into the memory pointed to by dest,
+ *   increment the dest storage pointer,
+ *   recursively call save_varis on the rest of the arguments,
+ *   and return the final value of the dest storage pointer.
+ *
+ * @tparam Pargs Types of remaining arguments
+ * @param[in, out] dest Pointer to where vari pointers are saved
+ * @param[in] x A var
+ * @param[in] args Additional arguments to have their varis saved
+ * @return Final position of dest pointer
+ */
+template <typename... Pargs>
+inline vari** save_varis(vari** dest, const var& x, Pargs&&... args) {
+  *dest = x.vi_;
+  return save_varis(dest + 1, std::forward<Pargs>(args)...);
+}
+
+/**
+ * Save the vari pointers in x into the memory pointed to by dest,
+ *   increment the dest storage pointer,
+ *   recursively call save_varis on the rest of the arguments,
+ *   and return the final value of the dest storage pointer.
+ *
+ * @tparam VarVec A variant of std::vector<var>
+ * @tparam Pargs Types of remaining arguments
+ * @param[in, out] dest Pointer to where vari pointers are saved
+ * @param[in] x A std::vector of vars
+ * @param[in] args Additional arguments to have their varis saved
+ * @return Final position of dest pointer
+ */
+template <typename VarVec, require_std_vector_vt<is_var, VarVec>*,
+          typename... Pargs>
+inline vari** save_varis(vari** dest, VarVec&& x, Pargs&&... args) {
+  for (int i = 0; i < x.size(); ++i) {
+    dest[i] = x[i].vi_;
+  }
+  return save_varis(dest + x.size(), std::forward<Pargs>(args)...);
+}
+
+/**
+ * Save the vari pointers in x into the memory pointed to by dest,
+ *   increment the dest storage pointer,
+ *   recursively call save_varis on the rest of the arguments,
+ *   and return the final value of the dest storage pointer.
+ *
+ * @tparam VecContainer std::vector<T> where T is another type containing vars
+ * @tparam Pargs Types of remaining arguments
+ * @param[in, out] dest Pointer to where vari pointers are saved
+ * @param[in] x A std::vector of containers containing of vars
+ * @param[in] args Additional arguments to have their varis saved
+ * @return Final position of dest pointer
+ */
+template <typename VecContainer, require_std_vector_st<is_var, VecContainer>*,
+          require_std_vector_vt<is_container, VecContainer>*, typename... Pargs>
+inline vari** save_varis(vari** dest, VecContainer&& x, Pargs&&... args) {
+  for (size_t i = 0; i < x.size(); ++i) {
+    dest = save_varis(dest, x[i]);
+  }
+  return save_varis(dest, std::forward<Pargs>(args)...);
+}
+
+/**
+ * Save the vari pointers in x into the memory pointed to by dest,
+ *   increment the dest storage pointer,
+ *   recursively call save_varis on the rest of the arguments,
+ *   and return the final value of the dest storage pointer.
+ *
+ * @tparam EigT An Eigen type with var value type
+ * @tparam Pargs Types of remaining arguments
+ * @param[in, out] dest Pointer to where vari pointers are saved
+ * @param[in] x An Eigen container of vars
+ * @param[in] args Additional arguments to have their varis saved
+ * @return Final position of dest pointer
+ */
+template <typename EigT, require_eigen_vt<is_var, EigT>*, typename... Pargs>
+inline vari** save_varis(vari** dest, EigT&& x, Pargs&&... args) {
+  for (int i = 0; i < x.size(); ++i) {
+    dest[i] = x(i).vi_;
+  }
+  return save_varis(dest + x.size(), std::forward<Pargs>(args)...);
+}
+
+/**
+ * Ignore arithmetic types.
+ *
+ * Recursively call save_varis on the rest of the arguments
+ *   and return the final value of the dest storage pointer.
+ *
+ * @tparam Arith An arithmetic type
+ * @tparam Pargs Types of remaining arguments
+ * @param[in, out] dest Pointer to where vari pointers are saved
+ * @param[in] x An argument not containing vars
+ * @param[in] args Additional arguments to have their varis saved
+ * @return Final position of dest pointer
+ */
+template <typename Arith, require_arithmetic_t<scalar_type_t<Arith>>*,
+          typename... Pargs>
+inline vari** save_varis(vari** dest, Arith&& x, Pargs&&... args) {
+  return save_varis(dest, std::forward<Pargs>(args)...);
+}
+
+/**
+ * End save_varis recursion and return pointer
+ *
+ * @param dest Pointer
+ */
+inline vari** save_varis(vari** dest) { return dest; }
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/prim/functor/apply_test.cpp b/test/unit/math/prim/functor/apply_test.cpp
index 429da239920..1eabbc1f315 100644
--- a/test/unit/math/prim/functor/apply_test.cpp
+++ b/test/unit/math/prim/functor/apply_test.cpp
@@ -11,6 +11,15 @@ struct func {
   }
 };
 
+TEST(MathFunctions, apply_basic_empty) {
+  std::tuple<> x;
+
+  auto y = stan::math::apply([]() { return static_cast<double>(1.7); }, x);
+
+  EXPECT_EQ(1.7, y);
+  EXPECT_TRUE((std::is_same<double, decltype(y)>::value));
+}
+
 TEST(MathFunctions, apply_basic_double) {
   std::tuple<double> x = std::make_tuple(1.0);
 
diff --git a/test/unit/math/rev/core/accumulate_adjoints_test.cpp b/test/unit/math/rev/core/accumulate_adjoints_test.cpp
new file mode 100644
index 00000000000..644282f88eb
--- /dev/null
+++ b/test/unit/math/rev/core/accumulate_adjoints_test.cpp
@@ -0,0 +1,438 @@
+#include <gtest/gtest.h>
+#include <stan/math/rev/core.hpp>
+#include <stan/math.hpp>
+#include <vector>
+
+TEST(AgradRev_accumulate_adjoints, zero_args) {
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data());
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, int_arg) {
+  int arg = 5;
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, double_arg) {
+  double arg = 5.0;
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_int_arg) {
+  std::vector<int> arg(5, 10);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_double_arg) {
+  std::vector<double> arg(5, 10.0);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, eigen_vector_arg) {
+  Eigen::VectorXd arg = Eigen::VectorXd::Ones(5);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, eigen_row_vector_arg) {
+  Eigen::RowVectorXd arg = Eigen::RowVectorXd::Ones(5);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, eigen_matrix_arg) {
+  Eigen::MatrixXd arg = Eigen::MatrixXd::Ones(5, 5);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_std_vector_double_arg) {
+  std::vector<std::vector<double>> arg(5, std::vector<double>(5, 10.0));
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_eigen_vector_arg) {
+  std::vector<Eigen::VectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_eigen_row_vector_arg) {
+  std::vector<Eigen::RowVectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_eigen_matrix_arg) {
+  std::vector<Eigen::MatrixXd> arg(2, Eigen::MatrixXd::Ones(5, 3));
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_accumulate_adjoints, var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  var arg(5.0);
+  arg.vi_->adj_ = 1.0;
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 1;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  std::vector<var> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg[i] = 5.0;
+    arg[i].vi_->adj_ = i + 1.0;
+  }
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 5;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), i + 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, eigen_vector_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = 5.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 5;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), i + 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, eigen_row_vector_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = 5.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 5;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), i + 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, eigen_matrix_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg(5, 5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = 5.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 25;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), i + 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_std_vector_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  std::vector<var> arg_(5, var(5.0));
+  for (size_t i = 0; i < arg_.size(); ++i)
+    arg_[i].vi_->adj_ = 1.0;
+  std::vector<std::vector<var>> arg(5, arg_);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 25;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_eigen_vector_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg_(5);
+  for (size_t i = 0; i < arg_.size(); ++i) {
+    arg_(i) = 5.0;
+    arg_(i).vi_->adj_ = 1.0;
+  }
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, 1>> arg(2, arg_);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 10;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_eigen_row_vector_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg_(5);
+  for (size_t i = 0; i < arg_.size(); ++i) {
+    arg_(i) = 5.0;
+    arg_(i).vi_->adj_ = 1.0;
+  }
+  std::vector<Eigen::Matrix<var, 1, Eigen::Dynamic>> arg(2, arg_);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 10;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, std_vector_eigen_matrix_var_arg) {
+  using stan::math::var;
+  using stan::math::vari;
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg_(5, 3);
+  for (size_t i = 0; i < arg_.size(); ++i) {
+    arg_(i) = 5.0;
+    arg_(i).vi_->adj_ = 1.0;
+  }
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>> arg(2, arg_);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(storage.data(), arg);
+
+  const size_t num_vars = 30;
+
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_accumulate_adjoints, sum) {
+  using stan::math::var;
+  using stan::math::vari;
+  int arg1 = 1;
+  double arg2 = 1.0;
+  std::vector<int> arg3(5, 1);
+  std::vector<double> arg4(5, 1.0);
+  Eigen::VectorXd arg5 = Eigen::VectorXd::Ones(5);
+  Eigen::RowVectorXd arg6 = Eigen::RowVectorXd::Ones(5);
+  Eigen::MatrixXd arg7 = Eigen::MatrixXd::Ones(5, 5);
+  std::vector<std::vector<double>> arg8(2, arg4);
+  std::vector<Eigen::VectorXd> arg9(2, arg5);
+
+  var arg10(new vari(5.0));
+  arg10.vi_->adj_ = 1.0;
+  std::vector<var> arg11(5, new vari(5.0, 1.0));
+  for (size_t i = 0; i < arg11.size(); ++i) {
+    arg11[i].vi_->adj_ = 1.0;
+  }
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg12(3);
+  for (size_t i = 0; i < arg12.size(); ++i) {
+    arg12(i) = 5.0;
+    arg12(i).vi_->adj_ = 1.0;
+  }
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg13(4);
+  for (size_t i = 0; i < arg13.size(); ++i) {
+    arg13(i) = 5.0;
+    arg13(i).vi_->adj_ = 1.0;
+  }
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg14(5, 3);
+  for (size_t i = 0; i < arg14.size(); ++i) {
+    arg14(i) = 5.0;
+    arg14(i).vi_->adj_ = 1.0;
+  }
+  std::vector<std::vector<var>> arg15(2, arg11);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, 1>> arg16(2, arg12);
+  std::vector<Eigen::Matrix<var, 1, Eigen::Dynamic>> arg17(2, arg13);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>> arg18(2,
+                                                                        arg14);
+
+  Eigen::VectorXd storage = Eigen::VectorXd::Zero(1000);
+  double* ptr = stan::math::accumulate_adjoints(
+      storage.data(), arg1, arg18, arg17, arg2, arg16, arg3, arg15, arg4, arg14,
+      arg5, arg13, arg12, arg6, arg11, arg7, arg10, arg8, arg9);
+
+  const size_t num_vars = 82;
+  for (int i = 0; i < num_vars; ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 1.0);
+  }
+  for (int i = num_vars; i < storage.size(); ++i) {
+    EXPECT_FLOAT_EQ(storage(i), 0.0);
+  }
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+  stan::math::recover_memory();
+}
diff --git a/test/unit/math/rev/core/count_vars_test.cpp b/test/unit/math/rev/core/count_vars_test.cpp
new file mode 100644
index 00000000000..ec9a6b96e38
--- /dev/null
+++ b/test/unit/math/rev/core/count_vars_test.cpp
@@ -0,0 +1,159 @@
+#include <gtest/gtest.h>
+#include <stan/math/rev/core.hpp>
+#include <stan/math.hpp>
+#include <vector>
+
+using stan::math::var;
+
+TEST(AgradRev_count_vars, int_arg) {
+  int arg = 5;
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, double_arg) {
+  double arg = 5.0;
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_int_arg) {
+  std::vector<int> arg(5, 10);
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_double_arg) {
+  std::vector<double> arg(5, 10.0);
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, eigen_vector_arg) {
+  Eigen::VectorXd arg = Eigen::VectorXd::Ones(5);
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, eigen_row_vector_arg) {
+  Eigen::RowVectorXd arg = Eigen::RowVectorXd::Ones(5);
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, eigen_matrix_arg) {
+  Eigen::MatrixXd arg = Eigen::MatrixXd::Ones(5, 5);
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_std_vector_double_arg) {
+  std::vector<std::vector<double>> arg(5, std::vector<double>(5, 10.0));
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_eigen_vector_arg) {
+  std::vector<Eigen::VectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_eigen_row_vector_arg) {
+  std::vector<Eigen::RowVectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_eigen_matrix_arg) {
+  std::vector<Eigen::MatrixXd> arg(2, Eigen::MatrixXd::Ones(5, 3));
+
+  EXPECT_EQ(0, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, var_arg) {
+  var arg = 5.0;
+
+  EXPECT_EQ(1, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_var_arg) {
+  std::vector<var> arg(5);
+
+  EXPECT_EQ(5, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, eigen_vector_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg(5);
+
+  EXPECT_EQ(5, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, eigen_row_vector_var_arg) {
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg(5);
+
+  EXPECT_EQ(5, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, eigen_matrix_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg(5, 5);
+
+  EXPECT_EQ(25, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_std_vector_var_arg) {
+  std::vector<std::vector<var>> arg(5, std::vector<var>(5));
+
+  EXPECT_EQ(25, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_eigen_vector_var_arg) {
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, 1>> arg(
+      2, Eigen::Matrix<var, Eigen::Dynamic, 1>(5));
+
+  EXPECT_EQ(10, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_eigen_row_vector_var_arg) {
+  std::vector<Eigen::Matrix<var, 1, Eigen::Dynamic>> arg(
+      2, Eigen::Matrix<var, 1, Eigen::Dynamic>(5));
+
+  EXPECT_EQ(10, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, std_vector_eigen_matrix_var_arg) {
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>> arg(
+      2, Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>(5, 3));
+
+  EXPECT_EQ(30, stan::math::count_vars(arg));
+}
+
+TEST(AgradRev_count_vars, zero_args) { EXPECT_EQ(0, stan::math::count_vars()); }
+
+TEST(AgradRev_count_vars, sum) {
+  int arg1 = 1;
+  double arg2 = 1.0;
+  std::vector<int> arg3(5, 1);
+  std::vector<double> arg4(5, 1.0);
+  Eigen::VectorXd arg5 = Eigen::VectorXd::Ones(5);
+  Eigen::RowVectorXd arg6 = Eigen::RowVectorXd::Ones(5);
+  Eigen::MatrixXd arg7 = Eigen::MatrixXd::Ones(5, 5);
+  std::vector<std::vector<double>> arg8(2, arg4);
+  std::vector<Eigen::VectorXd> arg9(2, arg5);
+
+  var arg10 = 1.0;
+  std::vector<var> arg11(5, 1.0);
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg12(3);
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg13(4);
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg14(5, 3);
+  std::vector<std::vector<var>> arg15(2, arg11);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, 1>> arg16(2, arg12);
+  std::vector<Eigen::Matrix<var, 1, Eigen::Dynamic>> arg17(2, arg13);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>> arg18(2,
+                                                                        arg14);
+
+  EXPECT_EQ(
+      1 + 5 + 3 + 4 + 15 + 2 * 5 + 2 * 3 + 2 * 4 + 30,
+      count_vars(arg1, arg18, arg17, arg2, arg16, arg3, arg15, arg4, arg14,
+                 arg5, arg13, arg12, arg6, arg11, arg7, arg10, arg8, arg9));
+}
diff --git a/test/unit/math/rev/core/deep_copy_vars_test.cpp b/test/unit/math/rev/core/deep_copy_vars_test.cpp
new file mode 100644
index 00000000000..ac9e067a589
--- /dev/null
+++ b/test/unit/math/rev/core/deep_copy_vars_test.cpp
@@ -0,0 +1,307 @@
+#include <gtest/gtest.h>
+#include <stan/math/rev/core.hpp>
+#include <stan/math.hpp>
+#include <vector>
+
+using stan::math::var;
+using stan::math::vari;
+
+TEST(AgradRev_deep_copy_vars, int_arg) {
+  int arg = 5;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, double_arg) {
+  double arg = 5.0;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_int_arg) {
+  std::vector<int> arg(5, 10);
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_double_arg) {
+  std::vector<double> arg(5, 10.0);
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, eigen_vector_arg) {
+  Eigen::VectorXd arg = Eigen::VectorXd::Ones(5);
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, eigen_row_vector_arg) {
+  Eigen::RowVectorXd arg = Eigen::RowVectorXd::Ones(5);
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, eigen_matrix_arg) {
+  Eigen::MatrixXd arg = Eigen::MatrixXd::Ones(5, 5);
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_std_vector_double_arg) {
+  std::vector<std::vector<double>> arg(5, std::vector<double>(5, 10.0));
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_eigen_vector_arg) {
+  std::vector<Eigen::VectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_eigen_row_vector_arg) {
+  std::vector<Eigen::RowVectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_eigen_matrix_arg) {
+  std::vector<Eigen::MatrixXd> arg(2, Eigen::MatrixXd::Ones(5, 3));
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  EXPECT_EQ(&out, &arg);
+}
+
+TEST(AgradRev_deep_copy_vars, var_arg) {
+  var arg(5.0);
+  arg.vi_->adj_ = 2.0;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  out.grad();
+
+  EXPECT_EQ(out.vi_->adj_, 1.0);
+  EXPECT_EQ(arg.vi_->adj_, 2.0);
+  EXPECT_EQ(out, arg);
+  EXPECT_NE(out.vi_, arg.vi_);
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_var_arg) {
+  std::vector<var> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i)
+    arg[i] = i + 1.0;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i) {
+    stan::math::set_zero_all_adjoints();
+    arg[i].vi_->adj_ = 2.0;
+
+    out[i].grad();
+
+    EXPECT_EQ(out[i].vi_->adj_, 1.0);
+    EXPECT_EQ(arg[i].vi_->adj_, 2.0);
+    EXPECT_EQ(out[i], arg[i]);
+    EXPECT_NE(out[i].vi_, arg[i].vi_);
+  }
+}
+
+TEST(AgradRev_deep_copy_vars, eigen_vector_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = i + 1.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i) {
+    stan::math::set_zero_all_adjoints();
+    arg(i).vi_->adj_ = 2.0;
+
+    out(i).grad();
+
+    EXPECT_EQ(out(i).vi_->adj_, 1.0);
+    EXPECT_EQ(arg(i).vi_->adj_, 2.0);
+    EXPECT_EQ(out(i), arg(i));
+    EXPECT_NE(out(i).vi_, arg(i).vi_);
+  }
+}
+
+TEST(AgradRev_deep_copy_vars, eigen_row_vector_var_arg) {
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = i + 1.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i) {
+    stan::math::set_zero_all_adjoints();
+    arg(i).vi_->adj_ = 2.0;
+
+    out(i).grad();
+
+    EXPECT_EQ(out(i).vi_->adj_, 1.0);
+    EXPECT_EQ(arg(i).vi_->adj_, 2.0);
+    EXPECT_EQ(out(i), arg(i));
+    EXPECT_NE(out(i).vi_, arg(i).vi_);
+  }
+}
+
+TEST(AgradRev_deep_copy_vars, eigen_matrix_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg(5, 5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = i + 1.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i) {
+    stan::math::set_zero_all_adjoints();
+    arg(i).vi_->adj_ = 2.0;
+
+    out(i).grad();
+
+    EXPECT_EQ(out(i).vi_->adj_, 1.0);
+    EXPECT_EQ(arg(i).vi_->adj_, 2.0);
+    EXPECT_EQ(out(i), arg(i));
+    EXPECT_NE(out(i).vi_, arg(i).vi_);
+  }
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_std_vector_var_arg) {
+  std::vector<var> arg_(5);
+  std::vector<std::vector<var>> arg(5, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i][j] = i * arg[i].size() + j + 5.0;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i)
+    for (int j = 0; j < arg[i].size(); ++j) {
+      stan::math::set_zero_all_adjoints();
+      arg[i][j].vi_->adj_ = 2.0;
+
+      out[i][j].grad();
+
+      EXPECT_EQ(out[i][j].vi_->adj_, 1.0);
+      EXPECT_EQ(arg[i][j].vi_->adj_, 2.0);
+      EXPECT_EQ(out[i][j], arg[i][j]);
+      EXPECT_NE(out[i][j].vi_, arg[i][j].vi_);
+    }
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_eigen_vector_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg_(5);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, 1>> arg(2, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i](j) = i * arg[i].size() + j + 5.0;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i)
+    for (int j = 0; j < arg[i].size(); ++j) {
+      stan::math::set_zero_all_adjoints();
+      arg[i](j).vi_->adj_ = 2.0;
+
+      out[i](j).grad();
+
+      EXPECT_EQ(out[i](j).vi_->adj_, 1.0);
+      EXPECT_EQ(arg[i](j).vi_->adj_, 2.0);
+      EXPECT_EQ(out[i](j), arg[i](j));
+      EXPECT_NE(out[i](j).vi_, arg[i](j).vi_);
+    }
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_eigen_row_vector_var_arg) {
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg_(5);
+  std::vector<Eigen::Matrix<var, 1, Eigen::Dynamic>> arg(2, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i](j) = i * arg[i].size() + j + 5.0;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i)
+    for (int j = 0; j < arg[i].size(); ++j) {
+      stan::math::set_zero_all_adjoints();
+      arg[i](j).vi_->adj_ = 2.0;
+
+      out[i](j).grad();
+
+      EXPECT_EQ(out[i](j).vi_->adj_, 1.0);
+      EXPECT_EQ(arg[i](j).vi_->adj_, 2.0);
+      EXPECT_EQ(out[i](j), arg[i](j));
+      EXPECT_NE(out[i](j).vi_, arg[i](j).vi_);
+    }
+}
+
+TEST(AgradRev_deep_copy_vars, std_vector_eigen_matrix_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg_(5, 3);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>> arg(2, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i](j) = i * arg[i].size() + j + 5.0;
+
+  decltype(stan::math::deep_copy_vars(arg)) out
+      = stan::math::deep_copy_vars(arg);
+
+  for (int i = 0; i < arg.size(); ++i)
+    for (int j = 0; j < arg[i].size(); ++j) {
+      stan::math::set_zero_all_adjoints();
+      arg[i](j).vi_->adj_ = 2.0;
+
+      out[i](j).grad();
+
+      EXPECT_EQ(out[i](j).vi_->adj_, 1.0);
+      EXPECT_EQ(arg[i](j).vi_->adj_, 2.0);
+      EXPECT_EQ(out[i](j), arg[i](j));
+      EXPECT_NE(out[i](j).vi_, arg[i](j).vi_);
+    }
+}
diff --git a/test/unit/math/rev/core/save_varis_test.cpp b/test/unit/math/rev/core/save_varis_test.cpp
new file mode 100644
index 00000000000..e7740d15cbc
--- /dev/null
+++ b/test/unit/math/rev/core/save_varis_test.cpp
@@ -0,0 +1,425 @@
+#include <stan/math.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+using stan::math::var;
+using stan::math::vari;
+
+TEST(AgradRev_save_varis, zero_args) {
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data());
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_save_varis, int_arg) {
+  int arg = 5;
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+  stan::math::recover_memory();
+}
+
+TEST(AgradRev_save_varis, double_arg) {
+  double arg = 5.0;
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, std_vector_int_arg) {
+  std::vector<int> arg(5, 10);
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, std_vector_double_arg) {
+  std::vector<double> arg(5, 10.0);
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, eigen_vector_arg) {
+  Eigen::VectorXd arg = Eigen::VectorXd::Ones(5);
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, eigen_row_vector_arg) {
+  Eigen::RowVectorXd arg = Eigen::RowVectorXd::Ones(5);
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, eigen_matrix_arg) {
+  Eigen::MatrixXd arg = Eigen::MatrixXd::Ones(5, 5);
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, std_vector_std_vector_double_arg) {
+  std::vector<std::vector<double>> arg(5, std::vector<double>(5, 10.0));
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, std_vector_eigen_vector_arg) {
+  std::vector<Eigen::VectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, std_vector_eigen_row_vector_arg) {
+  std::vector<Eigen::RowVectorXd> arg(2, Eigen::VectorXd::Ones(5));
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, std_vector_eigen_matrix_arg) {
+  std::vector<Eigen::MatrixXd> arg(2, Eigen::MatrixXd::Ones(5, 3));
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  for (int i = 0; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data());
+}
+
+TEST(AgradRev_save_varis, var_arg) {
+  var arg(5.0);
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  for (int i = 0; i < num_vars; ++i)
+    EXPECT_EQ(storage[i], arg.vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, std_vector_var_arg) {
+  std::vector<var> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i)
+    arg[i] = 5.0;
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  for (int i = 0; i < num_vars; ++i)
+    EXPECT_EQ(storage[i], arg[i].vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, eigen_vector_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = 5.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  for (int i = 0; i < num_vars; ++i)
+    EXPECT_EQ(storage[i], arg(i).vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, eigen_row_vector_var_arg) {
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg(5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = 5.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  for (int i = 0; i < num_vars; ++i)
+    EXPECT_EQ(storage[i], arg(i).vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, eigen_matrix_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg(5, 5);
+  for (size_t i = 0; i < arg.size(); ++i) {
+    arg(i) = 5.0;
+    arg(i).vi_->adj_ = i + 1.0;
+  }
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  for (int i = 0; i < num_vars; ++i)
+    EXPECT_EQ(storage[i], arg(i).vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, std_vector_std_vector_var_arg) {
+  std::vector<var> arg_(5);
+  std::vector<std::vector<var>> arg(5, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i][j] = 5.0;
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  EXPECT_EQ(arg.size() * arg[0].size(), num_vars);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      EXPECT_EQ(storage[i * arg[0].size() + j], arg[i][j].vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, std_vector_eigen_vector_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg_(5);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, 1>> arg(2, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i](j) = 5.0;
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  EXPECT_EQ(arg.size() * arg[0].size(), num_vars);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      EXPECT_EQ(storage[i * arg[0].size() + j], arg[i](j).vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, std_vector_eigen_row_vector_var_arg) {
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg_(5);
+  std::vector<Eigen::Matrix<var, 1, Eigen::Dynamic>> arg(2, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i](j) = 5.0;
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  EXPECT_EQ(arg.size() * arg[0].size(), num_vars);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      EXPECT_EQ(storage[i * arg[0].size() + j], arg[i](j).vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, std_vector_eigen_matrix_var_arg) {
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg_(5, 3);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>> arg(2, arg_);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      arg[i](j) = 5.0;
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(storage.data(), arg);
+
+  size_t num_vars = stan::math::count_vars(arg);
+
+  EXPECT_EQ(arg.size() * arg[0].size(), num_vars);
+  for (size_t i = 0; i < arg.size(); ++i)
+    for (size_t j = 0; j < arg[i].size(); ++j)
+      EXPECT_EQ(storage[i * arg[0].size() + j], arg[i](j).vi_);
+  for (int i = num_vars; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}
+
+TEST(AgradRev_save_varis, sum) {
+  int arg1 = 1;
+  double arg2 = 1.0;
+  std::vector<int> arg3(5, 1);
+  std::vector<double> arg4(5, 1.0);
+  Eigen::VectorXd arg5 = Eigen::VectorXd::Ones(5);
+  Eigen::RowVectorXd arg6 = Eigen::RowVectorXd::Ones(5);
+  Eigen::MatrixXd arg7 = Eigen::MatrixXd::Ones(5, 5);
+  std::vector<std::vector<double>> arg8(2, arg4);
+  std::vector<Eigen::VectorXd> arg9(2, arg5);
+
+  var arg10(new vari(5.0));
+  arg10.vi_->adj_ = 1.0;
+  std::vector<var> arg11(5, new vari(5.0, 1.0));
+  for (size_t i = 0; i < arg11.size(); ++i)
+    arg11[i].vi_->adj_ = 1.0;
+  Eigen::Matrix<var, Eigen::Dynamic, 1> arg12(3);
+  for (size_t i = 0; i < arg12.size(); ++i) {
+    arg12(i) = 5.0;
+    arg12(i).vi_->adj_ = 1.0;
+  }
+  Eigen::Matrix<var, 1, Eigen::Dynamic> arg13(4);
+  for (size_t i = 0; i < arg13.size(); ++i) {
+    arg13(i) = 5.0;
+    arg13(i).vi_->adj_ = 1.0;
+  }
+  Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> arg14(5, 3);
+  for (size_t i = 0; i < arg14.size(); ++i) {
+    arg14(i) = 5.0;
+    arg14(i).vi_->adj_ = 1.0;
+  }
+  std::vector<std::vector<var>> arg15(2, arg11);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, 1>> arg16(2, arg12);
+  std::vector<Eigen::Matrix<var, 1, Eigen::Dynamic>> arg17(2, arg13);
+  std::vector<Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>> arg18(2,
+                                                                        arg14);
+
+  std::vector<vari*> storage(1000, nullptr);
+  vari** ptr = stan::math::save_varis(
+      storage.data(), arg1, arg18, arg17, arg2, arg16, arg3, arg15, arg4, arg14,
+      arg5, arg13, arg12, arg6, arg11, arg7, arg10, arg8, arg9);
+
+  size_t num_vars = stan::math::count_vars(
+      arg1, arg18, arg17, arg2, arg16, arg3, arg15, arg4, arg14, arg5, arg13,
+      arg12, arg6, arg11, arg7, arg10, arg8, arg9);
+
+  int total = 0;
+  for (int i = 0; i < arg18.size(); ++i)
+    for (int j = 0; j < arg18[i].size(); ++j) {
+      EXPECT_EQ(storage[total], arg18[i](j).vi_);
+      total++;
+    }
+  for (int i = 0; i < arg17.size(); ++i)
+    for (int j = 0; j < arg17[i].size(); ++j) {
+      EXPECT_EQ(storage[total], arg17[i](j).vi_);
+      total++;
+    }
+  for (int i = 0; i < arg16.size(); ++i)
+    for (int j = 0; j < arg16[i].size(); ++j) {
+      EXPECT_EQ(storage[total], arg16[i](j).vi_);
+      total++;
+    }
+  for (int i = 0; i < arg15.size(); ++i)
+    for (int j = 0; j < arg15[i].size(); ++j) {
+      EXPECT_EQ(storage[total], arg15[i][j].vi_);
+      total++;
+    }
+  for (int i = 0; i < arg14.size(); ++i) {
+    EXPECT_EQ(storage[total], arg14(i).vi_);
+    total++;
+  }
+  for (int i = 0; i < arg13.size(); ++i) {
+    EXPECT_EQ(storage[total], arg13(i).vi_);
+    total++;
+  }
+  for (int i = 0; i < arg12.size(); ++i) {
+    EXPECT_EQ(storage[total], arg12(i).vi_);
+    total++;
+  }
+  for (int i = 0; i < arg11.size(); ++i) {
+    EXPECT_EQ(storage[total], arg11[i].vi_);
+    total++;
+  }
+  EXPECT_EQ(storage[total], arg10.vi_);
+  total++;
+  EXPECT_EQ(total, num_vars);
+
+  for (int i = total; i < storage.size(); ++i)
+    EXPECT_EQ(storage[i], nullptr);
+
+  EXPECT_EQ(ptr, storage.data() + num_vars);
+}