diff --git a/stan/math/opencl/kernel_generator/load.hpp b/stan/math/opencl/kernel_generator/load.hpp
index 799cdb551a0..ae0f7d1eb6d 100644
--- a/stan/math/opencl/kernel_generator/load.hpp
+++ b/stan/math/opencl/kernel_generator/load.hpp
@@ -49,7 +49,7 @@ class load_
    * Creates a deep copy of this expression.
    * @return copy of \c *this
    */
-  inline load_<T&> deep_copy() const & { return load_<T&>(a_); }
+  inline load_<T&> deep_copy() const& { return load_<T&>(a_); }
   inline load_<T> deep_copy() && { return load_<T>(std::forward<T>(a_)); }
 
   /**
diff --git a/stan/math/prim/fun/add.hpp b/stan/math/prim/fun/add.hpp
index 218fa19ad6b..d20b0d2fa8a 100644
--- a/stan/math/prim/fun/add.hpp
+++ b/stan/math/prim/fun/add.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -22,7 +23,8 @@ namespace math {
  * dimensions.
  */
 template <typename Mat1, typename Mat2,
-          typename = require_all_eigen_t<Mat1, Mat2>>
+          require_all_eigen_t<Mat1, Mat2>* = nullptr,
+          require_all_not_eigen_vt<is_var, Mat1, Mat2>* = nullptr>
 inline auto add(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("add", "m1", m1, "m2", m2);
   return (m1 + m2).eval();
@@ -37,8 +39,10 @@ inline auto add(const Mat1& m1, const Mat2& m2) {
  * @param c Scalar.
  * @return The matrix plus the scalar.
  */
-template <typename Mat, typename Scal, typename = require_eigen_t<Mat>,
-          typename = require_stan_scalar_t<Scal>>
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_not_var_t<Scal>* = nullptr,
+          require_not_eigen_vt<is_var, Mat>* = nullptr>
 inline auto add(const Mat& m, const Scal c) {
   return (m.array() + c).matrix().eval();
 }
@@ -52,8 +56,10 @@ inline auto add(const Mat& m, const Scal c) {
  * @param m Matrix.
  * @return The scalar plus the matrix.
  */
-template <typename Scal, typename Mat, typename = require_stan_scalar_t<Scal>,
-          typename = require_eigen_t<Mat>>
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_not_var_t<Scal>* = nullptr,
+          require_not_eigen_vt<is_var, Mat>* = nullptr>
 inline auto add(const Scal c, const Mat& m) {
   return (c + m.array()).matrix().eval();
 }
diff --git a/stan/math/prim/fun/subtract.hpp b/stan/math/prim/fun/subtract.hpp
index 5ba3459d57d..fe5b44906cd 100644
--- a/stan/math/prim/fun/subtract.hpp
+++ b/stan/math/prim/fun/subtract.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -20,7 +21,8 @@ namespace math {
  * @return Difference between first matrix and second matrix.
  */
 template <typename Mat1, typename Mat2,
-          typename = require_all_eigen_t<Mat1, Mat2>>
+          require_all_eigen_t<Mat1, Mat2>* = nullptr,
+          require_all_not_eigen_vt<is_var, Mat1, Mat2>* = nullptr>
 inline auto subtract(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("subtract", "m1", m1, "m2", m2);
   return (m1 - m2).eval();
@@ -36,8 +38,10 @@ inline auto subtract(const Mat1& m1, const Mat2& m2) {
  * @param m Matrix or expression.
  * @return The scalar minus the matrix.
  */
-template <typename Scal, typename Mat, typename = require_stan_scalar_t<Scal>,
-          typename = require_eigen_t<Mat>>
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_not_var_t<Scal>* = nullptr,
+          require_not_eigen_vt<is_var, Mat>* = nullptr>
 inline auto subtract(const Scal c, const Mat& m) {
   return (c - m.array()).matrix().eval();
 }
@@ -52,8 +56,10 @@ inline auto subtract(const Scal c, const Mat& m) {
  * @param c Scalar.
  * @return The matrix minus the scalar.
  */
-template <typename Mat, typename Scal, typename = require_eigen_t<Mat>,
-          typename = require_stan_scalar_t<Scal>>
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_not_var_t<Scal>* = nullptr,
+          require_not_eigen_vt<is_var, Mat>* = nullptr>
 inline auto subtract(const Mat& m, const Scal c) {
   return (m.array() - c).matrix().eval();
 }
diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp
index a6025090f21..c0c31b3b5ae 100644
--- a/stan/math/rev/fun.hpp
+++ b/stan/math/rev/fun.hpp
@@ -12,6 +12,7 @@
 #include <stan/math/rev/fun/abs.hpp>
 #include <stan/math/rev/fun/acos.hpp>
 #include <stan/math/rev/fun/acosh.hpp>
+#include <stan/math/rev/fun/add.hpp>
 #include <stan/math/rev/fun/as_bool.hpp>
 #include <stan/math/rev/fun/arg.hpp>
 #include <stan/math/rev/fun/asin.hpp>
@@ -127,6 +128,7 @@
 #include <stan/math/rev/fun/squared_distance.hpp>
 #include <stan/math/rev/fun/stan_print.hpp>
 #include <stan/math/rev/fun/step.hpp>
+#include <stan/math/rev/fun/subtract.hpp>
 #include <stan/math/rev/fun/sum.hpp>
 #include <stan/math/rev/fun/tan.hpp>
 #include <stan/math/rev/fun/tanh.hpp>
diff --git a/stan/math/rev/fun/add.hpp b/stan/math/rev/fun/add.hpp
new file mode 100644
index 00000000000..8c892d0dc56
--- /dev/null
+++ b/stan/math/rev/fun/add.hpp
@@ -0,0 +1,538 @@
+#ifndef STAN_MATH_REV_FUN_MATRIX_ADD_HPP
+#define STAN_MATH_REV_FUN_MATRIX_ADD_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/prim.hpp>
+#include <type_traits>
+#include <iostream>
+
+namespace stan {
+namespace math {
+
+/**
+ * This is a subclass of the vari class for matrix
+ * addition A + B, where A and B are N by M matrices.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and B, and pointers to
+ * the varis for A and B if A or B is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + B.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tb type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Ta, typename Tb, int R, int C>
+class add_mat_vari : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double* Bd_;
+  vari** variRefA_;
+  vari** variRefB_;
+  vari** variRefAplusB_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memeory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param  B matrix
+   **/
+  add_mat_vari(const Eigen::Matrix<Ta, R, C>& A,
+               const Eigen::Matrix<Tb, R, C>& B)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(B.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAplusB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_vi>(variRefB_, rows_, cols_) = B.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Map<matrix_d> Bd(Bd_, rows_, cols_);
+    Ad = A.val();
+    Bd = B.val();
+    Map<matrix_vi>(variRefAplusB_, rows_, cols_)
+        = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAplusB(rows_, cols_);
+    adjAplusB = Map<matrix_vi>(variRefAplusB_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAplusB;
+    Map<matrix_vi>(variRefB_, rows_, cols_).adj() += adjAplusB;
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * addition A + B, where A is an N by M matrix of
+ * doubles and and B is an N by M matrix of vars.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and B, and pointers to
+ * the varis for A and B if A or B is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + B.
+ *
+ * @tparam Tb type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Tb, int R, int C>
+class add_mat_vari<double, Tb, R, C> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double* Bd_;
+  vari** variRefB_;
+  vari** variRefAplusB_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memeory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param  B matrix
+   **/
+  add_mat_vari(const Eigen::Matrix<double, R, C>& A,
+               const Eigen::Matrix<Tb, R, C>& B)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(B.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        variRefB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAplusB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefB_, rows_, cols_) = B.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Map<matrix_d> Bd(Bd_, rows_, cols_);
+    Ad = A;
+    Bd = B.val();
+    Map<matrix_vi>(variRefAplusB_, rows_, cols_)
+        = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAplusB(rows_, cols_);
+    adjAplusB = Map<matrix_vi>(variRefAplusB_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefB_, rows_, cols_).adj() += adjAplusB;
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * addition A + B, where A is an N by M matrix of vars
+ * and B is an N by M matrix of doubles.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and B, and pointers to
+ * the varis for A and B if A or B is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + B.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Ta, int R, int C>
+class add_mat_vari<Ta, double, R, C> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double* Bd_;
+  vari** variRefA_;
+  vari** variRefAplusB_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memeory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param  B matrix
+   **/
+  add_mat_vari(const Eigen::Matrix<Ta, R, C>& A,
+               const Eigen::Matrix<double, R, C>& B)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(B.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAplusB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Map<matrix_d> Bd(Bd_, rows_, cols_);
+    Ad = A.val();
+    Bd = B;
+    Map<matrix_vi>(variRefAplusB_, rows_, cols_)
+        = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAplusB(rows_, cols_);
+    adjAplusB = Map<matrix_vi>(variRefAplusB_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAplusB;
+  }
+};
+
+/**
+ * Return the sum of two matrices
+ * @tparam Mat1 type of first matrix
+ * @tparam Mat2 type of second matrix
+ *
+ * @param[in] A matrix
+ * @param[in] B matrix
+ * @return Sum of two matrices
+ */
+template <typename Mat1, typename Mat2,
+          require_all_eigen_t<Mat1, Mat2>* = nullptr,
+          require_any_eigen_vt<is_var, Mat1, Mat2>* = nullptr>
+auto add(const Mat1& A, const Mat2& B) {
+  using Ta = value_type_t<Mat1>;
+  using Tb = value_type_t<Mat2>;
+  constexpr int R = Mat1::RowsAtCompileTime;
+  constexpr int C = Mat1::ColsAtCompileTime;
+  check_matching_dims("add", "A", A, "B", B);
+  check_not_nan("add", "A", A);
+  check_not_nan("add", "B", B);
+
+  // Manage memory in arena allocator
+  add_mat_vari<Ta, Tb, R, C>* baseVari = new add_mat_vari<Ta, Tb, R, C>(A, B);
+  Eigen::Matrix<var, R, C> AplusB_v(A.rows(), A.cols());
+  AplusB_v.vi()
+      = Eigen::Map<matrix_vi>(&baseVari->variRefAplusB_[0], A.rows(), A.cols());
+
+  return AplusB_v;
+}
+
+/**
+ * This is a subclass of the vari class for matrix
+ * addition A + c, where A is a N by M matrix and
+ * c is a scalar.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and c, and pointers to
+ * the varis for A and c if A or c is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + c.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tc type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Ta, typename Tc, int R, int C>
+class add_mat_scal_vari : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double cd_;
+  vari** variRefA_;
+  vari* variRefc_;
+  vari** variRefAplusc_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param c scalar
+   **/
+  add_mat_scal_vari(const Eigen::Matrix<Ta, R, C>& A, const Tc c)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(A.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        cd_(c.val()),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefc_(c.vi_),
+        variRefAplusc_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    double cd = c.val();
+    Ad = A.val();
+    Map<matrix_vi>(variRefAplusc_, rows_, cols_)
+        = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) {
+            return new vari(x, false);
+          });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAplusc(rows_, cols_);
+    adjAplusc = Map<matrix_vi>(variRefAplusc_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAplusc;
+    variRefc_->adj_ += adjAplusc.sum();
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * addition A + c, where A is a N by M matrix and
+ * c is a double.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and c, and pointers to
+ * the varis for A and c if A or c is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + c.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tc type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Ta, int R, int C>
+class add_mat_scal_vari<Ta, double, R, C> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double cd_;
+  vari** variRefA_;
+  vari** variRefAplusc_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param c scalar
+   **/
+  add_mat_scal_vari(const Eigen::Matrix<Ta, R, C>& A, const double c)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(A.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        cd_(c),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAplusc_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    double cd = c;
+    Ad = A.val();
+    Map<matrix_vi>(variRefAplusc_, rows_, cols_)
+        = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) {
+            return new vari(x, false);
+          });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAplusc(rows_, cols_);
+    adjAplusc = Map<matrix_vi>(variRefAplusc_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAplusc;
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * addition A + c, where A is a N by M matrix of
+ * doubles c is a scalar var.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and c, and pointers to
+ * the varis for A and c if A or c is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + c.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tc type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Tc, int R, int C>
+class add_mat_scal_vari<double, Tc, R, C> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double cd_;
+  vari* variRefc_;
+  vari** variRefAplusc_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param c scalar
+   **/
+  add_mat_scal_vari(const Eigen::Matrix<double, R, C>& A, const Tc c)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(A.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        cd_(c.val()),
+        variRefc_(c.vi_),
+        variRefAplusc_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+
+    double cd = c.val();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Ad = A.val();
+    Map<matrix_vi>(variRefAplusc_, rows_, cols_)
+        = (Ad.array() + cd).matrix().unaryExpr([](double x) {
+            return new vari(x, false);
+          });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAplusc(rows_, cols_);
+    adjAplusc = Map<matrix_vi>(variRefAplusc_, rows_, cols_).adj();
+    variRefc_->adj_ += adjAplusc.sum();
+  }
+};
+
+/**
+ * Return the sum of the specified matrix and specified scalar.
+ *
+ * @tparam Mat type of the matrix or expression
+ * @tparam Scal type of the scalar
+ * @param m Matrix or expression.
+ * @param c Scalar.
+ * @return The matrix plus the scalar.
+ */
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_any_var_t<typename value_type<Mat>::type, Scal>* = nullptr>
+auto add(const Mat& A, const Scal c) {
+  using Ta = value_type_t<Mat>;
+  constexpr int R = Mat::RowsAtCompileTime;
+  constexpr int C = Mat::ColsAtCompileTime;
+  check_not_nan("add matrix-scalar A + c", "A", A);
+  check_not_nan("add matrix-scalar A + c", "c", c);
+
+  add_mat_scal_vari<Ta, Scal, R, C>* baseVari
+      = new add_mat_scal_vari<Ta, Scal, R, C>(A, c);
+  Eigen::Matrix<var, R, C> Aplusc_v(A.rows(), A.cols());
+  Aplusc_v.vi()
+      = Eigen::Map<matrix_vi>(&baseVari->variRefAplusc_[0], A.rows(), A.cols());
+
+  return Aplusc_v;
+}
+
+/**
+ * Return the sum of the specified scalar and specified matrix.
+ *
+ * @tparam Scal type of the scalar
+ * @tparam Mat type of the matrix or expression
+ * @param c Scalar.
+ * @param m Matrix.
+ * @return The scalar plus the matrix.
+ */
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_any_var_t<typename value_type<Mat>::type, Scal>* = nullptr>
+inline auto add(const Scal c, const Mat& A) {
+  return add(A, c);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/fun/subtract.hpp b/stan/math/rev/fun/subtract.hpp
new file mode 100644
index 00000000000..98821c9d80a
--- /dev/null
+++ b/stan/math/rev/fun/subtract.hpp
@@ -0,0 +1,591 @@
+#ifndef STAN_MATH_REV_FUN_MATRIX_SUBTRACT_HPP
+#define STAN_MATH_REV_FUN_MATRIX_SUBTRACT_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/prim.hpp>
+#include <type_traits>
+
+namespace stan {
+namespace math {
+
+/**
+ * This is a subclass of the vari class for matrix
+ * subtraction A - B, whre A and B are N by M matrices.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and B, and pointers to
+ * the varis for A and B if A or B is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A - B.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tb type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Ta, typename Tb, int R, int C>
+class subtract_mat_vari : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double* Bd_;
+  vari** variRefA_;
+  vari** variRefB_;
+  vari** variRefAminusB_;
+
+  /**
+   * Constructor for subtract_mat_vari.
+   *
+   * All memeory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param  B matrix
+   **/
+  subtract_mat_vari(const Eigen::Matrix<Ta, R, C>& A,
+                    const Eigen::Matrix<Tb, R, C>& B)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(B.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAminusB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_vi>(variRefB_, rows_, cols_) = B.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Map<matrix_d> Bd(Bd_, rows_, cols_);
+    Ad = A.val();
+    Bd = B.val();
+    Map<matrix_vi>(variRefAminusB_, rows_, cols_)
+        = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAminusB(rows_, cols_);
+    adjAminusB = Map<matrix_vi>(variRefAminusB_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAminusB;
+    Map<matrix_vi>(variRefB_, rows_, cols_).adj() -= adjAminusB;
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * subtraction A - B, where A is an N by M matrix of
+ * doubles and and B is an N by M matrix of vars.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and B, and pointers to
+ * the varis for A and B if A or B is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A - B.
+ *
+ * @tparam Tb type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Tb, int R, int C>
+class subtract_mat_vari<double, Tb, R, C> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double* Bd_;
+  vari** variRefB_;
+  vari** variRefAminusB_;
+
+  /**
+   * Constructor for subtract_mat_vari.
+   *
+   * All memeory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param  B matrix
+   **/
+  subtract_mat_vari(const Eigen::Matrix<double, R, C>& A,
+                    const Eigen::Matrix<Tb, R, C>& B)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(B.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        variRefB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAminusB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefB_, rows_, cols_) = B.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Map<matrix_d> Bd(Bd_, rows_, cols_);
+    Ad = A;
+    Bd = B.val();
+    Map<matrix_vi>(variRefAminusB_, rows_, cols_)
+        = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAminusB(rows_, cols_);
+    adjAminusB = Map<matrix_vi>(variRefAminusB_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefB_, rows_, cols_).adj() -= adjAminusB;
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * subtraction A - B, whre A is an N by M matrix of vars
+ * and B is an N by M matrix of doubles.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and B, and pointers to
+ * the varis for A and B if A or B is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A - B.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ */
+
+template <typename Ta, int R, int C>
+class subtract_mat_vari<Ta, double, R, C> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double* Bd_;
+  vari** variRefA_;
+  vari** variRefAminusB_;
+
+  /**
+   * Constructor for subtract_mat_vari.
+   *
+   * All memeory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param  B matrix
+   **/
+  subtract_mat_vari(const Eigen::Matrix<Ta, R, C>& A,
+                    const Eigen::Matrix<double, R, C>& B)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(B.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAminusB_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Map<matrix_d> Bd(Bd_, rows_, cols_);
+    Ad = A.val();
+    Bd = B;
+    Map<matrix_vi>(variRefAminusB_, rows_, cols_)
+        = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); });
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAminusB(rows_, cols_);
+    adjAminusB = Map<matrix_vi>(variRefAminusB_, rows_, cols_).adj();
+    Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAminusB;
+  }
+};
+
+/**
+ * Return the sum of two matrices
+ * @tparam Mat1 type of first matrix
+ * @tparam Mat2 type of second matrix
+ *
+ * @param[in] A matrix
+ * @param[in] B matrix
+ * @return Sum of two matrices
+ */
+template <typename Mat1, typename Mat2,
+          require_all_eigen_t<Mat1, Mat2>* = nullptr,
+          require_any_eigen_vt<is_var, Mat1, Mat2>* = nullptr>
+auto subtract(const Mat1& A, const Mat2& B) {
+  using Ta = value_type_t<Mat1>;
+  using Tb = value_type_t<Mat2>;
+  constexpr int R = Mat1::RowsAtCompileTime;
+  constexpr int C = Mat1::ColsAtCompileTime;
+  check_matching_dims("subtract", "A", A, "B", B);
+  check_not_nan("subtract", "A", A);
+  check_not_nan("subtract", "B", B);
+
+  // Manage memory in arena allocator
+  subtract_mat_vari<Ta, Tb, R, C>* baseVari
+      = new subtract_mat_vari<Ta, Tb, R, C>(A, B);
+  Eigen::Matrix<var, R, C> AminusB_v(A.rows(), A.cols());
+  AminusB_v.vi() = Eigen::Map<matrix_vi>(&baseVari->variRefAminusB_[0],
+                                         A.rows(), A.cols());
+
+  return AminusB_v;
+}
+
+/**
+ * This is a subclass of the vari class for matrix
+ * subtraction A + c, where A is a N by M matrix and
+ * c is a scalar.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and c, and pointers to
+ * the varis for A and c if A or c is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + c.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tc type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ * @tparam MatMinusScal: true if computing A - c, false if
+ * computing c-A
+ */
+
+template <typename Ta, typename Tc, int R, int C, bool MatMinusScal>
+class subtract_mat_scal_vari : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double cd_;
+  vari** variRefA_;
+  vari* variRefc_;
+  vari** variRefAminusc_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param c scalar
+   **/
+  subtract_mat_scal_vari(const Eigen::Matrix<Ta, R, C>& A, const Tc c)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(A.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        cd_(c.val()),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefc_(c.vi_),
+        variRefAminusc_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    double cd = c.val();
+    Ad = A.val();
+    if constexpr (MatMinusScal) {
+      Map<matrix_vi>(variRefAminusc_, rows_, cols_)
+          = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) {
+              return new vari(x, false);
+            });
+    } else {
+      Map<matrix_vi>(variRefAminusc_, rows_, cols_)
+          = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) {
+              return new vari(x, false);
+            });
+    }
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAminusc(rows_, cols_);
+    adjAminusc = Map<matrix_vi>(variRefAminusc_, rows_, cols_).adj();
+    if constexpr (MatMinusScal) {
+      Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAminusc;
+      variRefc_->adj_ -= adjAminusc.sum();
+    } else {
+      Map<matrix_vi>(variRefA_, rows_, cols_).adj() -= adjAminusc;
+      variRefc_->adj_ += adjAminusc.sum();
+    }
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * subtraction A + c, where A is a N by M matrix and
+ * c is a double.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and c, and pointers to
+ * the varis for A and c if A or c is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + c.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tc type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ * @tparam MatMinusScal: true if computing A - c, false if
+ * computing c-A
+ */
+
+template <typename Ta, int R, int C, bool MatMinusScal>
+class subtract_mat_scal_vari<Ta, double, R, C, MatMinusScal> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double cd_;
+  vari** variRefA_;
+  vari** variRefAminusc_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param c scalar
+   **/
+  subtract_mat_scal_vari(const Eigen::Matrix<Ta, R, C>& A, const double c)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(A.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        cd_(c),
+        variRefA_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
+        variRefAminusc_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+    Map<matrix_vi>(variRefA_, rows_, cols_) = A.vi();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    double cd = c;
+    Ad = A.val();
+    if constexpr (MatMinusScal) {
+      Map<matrix_vi>(variRefAminusc_, rows_, cols_)
+          = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) {
+              return new vari(x, false);
+            });
+    } else {
+      Map<matrix_vi>(variRefAminusc_, rows_, cols_)
+          = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) {
+              return new vari(x, false);
+            });
+    }
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAminusc(rows_, cols_);
+    adjAminusc = Map<matrix_vi>(variRefAminusc_, rows_, cols_).adj();
+
+    if constexpr (MatMinusScal) {
+      Map<matrix_vi>(variRefA_, rows_, cols_).adj() += adjAminusc;
+    } else {
+      Map<matrix_vi>(variRefA_, rows_, cols_).adj() -= adjAminusc;
+    }
+  }
+};
+
+/**
+ * This is a subclass of the vari class for matrix
+ * subtraction A + c, where A is a N by M matrix of
+ * doubles c is a scalar var.
+ *
+ * The class stores the structure of each matrix,
+ * the double values of A and c, and pointers to
+ * the varis for A and c if A or c is a var. It
+ * also instatiates and stores pointers to varis
+ * for all elements of A + c.
+ *
+ * @tparam Ta type of elements in matrix A
+ * @tparam Tc type of elements in matrix B
+ * @tparam R number of rows (common to A and B)
+ * @tparam C number of columns (common to A and B)
+ * @tparam MatMinusScal: true if computing A - c, false if
+ * computing c-A
+ */
+
+template <typename Tc, int R, int C, bool MatMinusScal>
+class subtract_mat_scal_vari<double, Tc, R, C, MatMinusScal> : public vari {
+ public:
+  int rows_;
+  int cols_;
+  int size_;
+  double* Ad_;
+  double cd_;
+  vari* variRefc_;
+  vari** variRefAminusc_;
+
+  /**
+   * Constructor for add_mat_vari.
+   *
+   * All memory is allocated in ChainableStack's
+   * stack_alloc arena.
+   *
+   * It is critical for the efficiency of this object
+   * that the constructor create new varis that aren't
+   * popped onto the var_stack_, but rather are
+   * popped onto the var_nochain_stack_. This is
+   * controlled by the second argument to
+   * vari's constructor.
+   *
+   * @param A matrix
+   * @param c scalar
+   **/
+  subtract_mat_scal_vari(const Eigen::Matrix<double, R, C>& A, const Tc c)
+      : vari(0.0),
+        rows_(A.rows()),
+        cols_(A.cols()),
+        size_(A.size()),
+        Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
+        cd_(c.val()),
+        variRefc_(c.vi_),
+        variRefAminusc_(
+            ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
+    using Eigen::Map;
+
+    double cd = c.val();
+    Map<matrix_d> Ad(Ad_, rows_, cols_);
+    Ad = A.val();
+    if constexpr (MatMinusScal) {
+      Map<matrix_vi>(variRefAminusc_, rows_, cols_)
+          = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) {
+              return new vari(x, false);
+            });
+    } else {
+      Map<matrix_vi>(variRefAminusc_, rows_, cols_)
+          = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) {
+              return new vari(x, false);
+            });
+    }
+  }
+
+  virtual void chain() {
+    using Eigen::Map;
+    matrix_d adjAminusc(rows_, cols_);
+    adjAminusc = Map<matrix_vi>(variRefAminusc_, rows_, cols_).adj();
+    if constexpr (MatMinusScal) {
+      variRefc_->adj_ -= adjAminusc.sum();
+    } else {
+      variRefc_->adj_ += adjAminusc.sum();
+    }
+  }
+};
+
+/**
+ * Return the difference A - c of a matrix A and a scalar c.
+ *
+ * @tparam Mat type of the matrix or expression
+ * @tparam Scal type of the scalar
+ * @param m Matrix or expression.
+ * @param c Scalar.
+ * @return The matrix minus the scalar.
+ */
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_any_var_t<typename value_type<Mat>::type, Scal>* = nullptr>
+auto subtract(const Mat& A, const Scal c) {
+  using Ta = value_type_t<Mat>;
+  constexpr int R = Mat::RowsAtCompileTime;
+  constexpr int C = Mat::ColsAtCompileTime;
+  check_not_nan("subtract matrix-scalar A - c", "A", A);
+  check_not_nan("subtract matrix-scalar A - c", "c", c);
+
+  subtract_mat_scal_vari<Ta, Scal, R, C, true>* baseVari
+      = new subtract_mat_scal_vari<Ta, Scal, R, C, true>(A, c);
+  Eigen::Matrix<var, R, C> Aminusc_v(A.rows(), A.cols());
+  Aminusc_v.vi() = Eigen::Map<matrix_vi>(&baseVari->variRefAminusc_[0],
+                                         A.rows(), A.cols());
+
+  return Aminusc_v;
+}
+
+/**
+ * Return the difference c - A of a scalar c and a matrix A.
+ *
+ * @tparam Scal type of the scalar
+ * @tparam Mat type of the matrix or expression
+ * @param c Scalar.
+ * @param m Matrix.
+ * @return The scalar minus the matrix.
+ */
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_any_var_t<typename value_type<Mat>::type, Scal>* = nullptr>
+auto subtract(const Scal c, const Mat& A) {
+  using Ta = value_type_t<Mat>;
+  constexpr int R = Mat::RowsAtCompileTime;
+  constexpr int C = Mat::ColsAtCompileTime;
+  check_not_nan("subtract matrix-scalar A - c", "A", A);
+  check_not_nan("subtract matrix-scalar A - c", "c", c);
+
+  subtract_mat_scal_vari<Ta, Scal, R, C, false>* baseVari
+      = new subtract_mat_scal_vari<Ta, Scal, R, C, false>(A, c);
+  Eigen::Matrix<var, R, C> Aminusc_v(A.rows(), A.cols());
+  Aminusc_v.vi() = Eigen::Map<matrix_vi>(&baseVari->variRefAminusc_[0],
+                                         A.rows(), A.cols());
+
+  return Aminusc_v;
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/rev/fun/add_test.cpp b/test/unit/math/rev/fun/add_test.cpp
new file mode 100644
index 00000000000..a27ed8143f5
--- /dev/null
+++ b/test/unit/math/rev/fun/add_test.cpp
@@ -0,0 +1,78 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+
+TEST(MathRevAdd, add_v_d) {
+  using stan::math::add;
+  using stan::math::matrix_d;
+  using stan::math::matrix_v;
+  using stan::math::var;
+  matrix_d Ad = Eigen::MatrixXd::Random(4, 6);
+  matrix_d Bd = Eigen::MatrixXd::Random(4, 6);
+  var c = 3.2;
+  var d = -2.73;
+  matrix_v Av = c * Ad;
+  matrix_v Bv = d * Bd;
+  matrix_v Cvd = add(Av, Bd);
+  matrix_v Cdv = add(Ad, Bv);
+  matrix_v Cvv = add(Av, Bv);
+  matrix_v Cdvd = add(Ad, add(Bv, Eigen::MatrixXd::Random(4, 6)));
+  matrix_v Adcv = add(c, Ad);
+  matrix_v Bvcv = add(c, Bv);
+  matrix_v Bvscal = add(0.7, Bv);
+  Cvd(0, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), Ad(0, 1));
+  EXPECT_FLOAT_EQ(d.adj(), 0.0);
+
+  stan::math::set_zero_all_adjoints();
+  Cdv(2, 3).grad();
+  EXPECT_FLOAT_EQ(d.adj(), Bd(2, 3));
+  EXPECT_FLOAT_EQ(c.adj(), 0.0);
+
+  stan::math::set_zero_all_adjoints();
+  Cvv(3, 2).grad();
+  EXPECT_FLOAT_EQ(d.adj(), Bd(3, 2));
+  EXPECT_FLOAT_EQ(c.adj(), Ad(3, 2));
+
+  stan::math::set_zero_all_adjoints();
+  Cdvd(3, 4).grad();
+  EXPECT_FLOAT_EQ(d.adj(), Bd(3, 4));
+  EXPECT_FLOAT_EQ(c.adj(), 0.0);
+
+  stan::math::set_zero_all_adjoints();
+  Adcv(1, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), 1.0);
+  EXPECT_FLOAT_EQ(d.adj(), 0.0);
+
+  stan::math::set_zero_all_adjoints();
+  Bvcv(2, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), 1.0);
+  EXPECT_FLOAT_EQ(d.adj(), Bd(2, 1));
+
+  stan::math::set_zero_all_adjoints();
+  Bvscal(2, 0).grad();
+  EXPECT_FLOAT_EQ(d.adj(), Bd(2, 0));
+  EXPECT_FLOAT_EQ(c.adj(), 0.0);
+}
+
+TEST(MathRevAdd, add_check_dim) {
+  using stan::math::add;
+  using stan::math::matrix_d;
+  using stan::math::matrix_v;
+  matrix_v A, B;
+
+  A.resize(3, 2);
+  B.resize(3, 2);
+  EXPECT_NO_THROW(add(A, B));
+
+  A.resize(1, 0);
+  B.resize(1, 0);
+  EXPECT_NO_THROW(add(A, B));
+
+  A.resize(0, 0);
+  B.resize(0, 0);
+  EXPECT_NO_THROW(add(A, B));
+
+  A.resize(2, 3);
+  B.resize(3, 2);
+  EXPECT_THROW(add(A, B), std::invalid_argument);
+};
diff --git a/test/unit/math/rev/fun/subtract_test.cpp b/test/unit/math/rev/fun/subtract_test.cpp
new file mode 100644
index 00000000000..9e73d317381
--- /dev/null
+++ b/test/unit/math/rev/fun/subtract_test.cpp
@@ -0,0 +1,99 @@
+#include <stan/math/rev.hpp>
+#include <gtest/gtest.h>
+
+TEST(MathRevSubtract, subtract_v_d) {
+  using stan::math::matrix_d;
+  using stan::math::matrix_v;
+  using stan::math::subtract;
+  using stan::math::var;
+  matrix_d Ad = Eigen::MatrixXd::Random(4, 6);
+  matrix_d Bd = Eigen::MatrixXd::Random(4, 6);
+  var c(3.2);
+  var d(-2.73);
+  matrix_v Av = c * Ad;
+  matrix_v Bv = d * Bd;
+  matrix_v Cvd = subtract(Av, Bd);
+  matrix_v Cdv = subtract(Ad, Bv);
+  matrix_v Cvv = subtract(Av, Bv);
+  matrix_v Cdvd = subtract(Ad, subtract(Bv, Eigen::MatrixXd::Random(4, 6)));
+  matrix_v Cddv = subtract(Ad, subtract(Eigen::MatrixXd::Random(4, 6), Bv));
+  matrix_v Adcv = subtract(Ad, c);
+  matrix_v cvAd = subtract(c, Ad);
+  matrix_v Bvcv = subtract(Bv, c);
+  matrix_v cvBv = subtract(c, Bv);
+  matrix_v Bvscal = subtract(Bv, 0.7);
+  matrix_v scalBv = subtract(0.7, Bv);
+
+  Cvd(0, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), Ad(0, 1));
+
+  stan::math::set_zero_all_adjoints();
+  Cdv(2, 3).grad();
+  EXPECT_FLOAT_EQ(d.adj(), -Bd(2, 3));
+
+  stan::math::set_zero_all_adjoints();
+  Cvv(3, 2).grad();
+  EXPECT_FLOAT_EQ(d.adj(), -Bd(3, 2));
+  EXPECT_FLOAT_EQ(c.adj(), Ad(3, 2));
+
+  stan::math::set_zero_all_adjoints();
+  Cdvd(3, 4).grad();
+  EXPECT_FLOAT_EQ(d.adj(), -Bd(3, 4));
+
+  stan::math::set_zero_all_adjoints();
+  Cddv(2, 5).grad();
+  EXPECT_FLOAT_EQ(d.adj(), Bd(2, 5));
+
+  stan::math::set_zero_all_adjoints();
+  Adcv(1, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), -1.0);
+  EXPECT_FLOAT_EQ(d.adj(), 0.0);
+
+  stan::math::set_zero_all_adjoints();
+  cvAd(1, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), 1.0);
+  EXPECT_FLOAT_EQ(d.adj(), 0.0);
+
+  stan::math::set_zero_all_adjoints();
+  Bvcv(2, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), -1.0);
+  EXPECT_FLOAT_EQ(d.adj(), Bd(2, 1));
+
+  stan::math::set_zero_all_adjoints();
+  cvBv(2, 1).grad();
+  EXPECT_FLOAT_EQ(c.adj(), 1.0);
+  EXPECT_FLOAT_EQ(d.adj(), -Bd(2, 1));
+
+  stan::math::set_zero_all_adjoints();
+  Bvscal(2, 0).grad();
+  EXPECT_FLOAT_EQ(d.adj(), Bd(2, 0));
+  EXPECT_FLOAT_EQ(c.adj(), 0.0);
+
+  stan::math::set_zero_all_adjoints();
+  scalBv(2, 0).grad();
+  EXPECT_FLOAT_EQ(d.adj(), -Bd(2, 0));
+  EXPECT_FLOAT_EQ(c.adj(), 0.0);
+}
+
+TEST(MathRevSubtract, subtract_check_dim) {
+  using stan::math::matrix_d;
+  using stan::math::matrix_v;
+  using stan::math::subtract;
+  matrix_v A, B;
+
+  A.resize(3, 2);
+  B.resize(3, 2);
+  EXPECT_NO_THROW(subtract(A, B));
+
+  A.resize(1, 0);
+  B.resize(1, 0);
+  EXPECT_NO_THROW(subtract(A, B));
+
+  A.resize(0, 0);
+  B.resize(0, 0);
+  EXPECT_NO_THROW(subtract(A, B));
+
+  A.resize(2, 3);
+  B.resize(3, 2);
+  EXPECT_THROW(subtract(A, B), std::invalid_argument);
+};