From 687ecce287b2cd2413a0d4ef710e1ac50363d7fc Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 15:32:09 +0200
Subject: [PATCH 1/9] backport batch 1

---
 stan/math/opencl/matrix_cl.hpp                |  81 ++++--
 .../opencl/prim/bernoulli_logit_glm_lpmf.hpp  |  18 +-
 .../math/opencl/prim/bernoulli_logit_lpmf.hpp |   2 +-
 stan/math/opencl/prim/bernoulli_lpmf.hpp      |   2 +-
 .../prim/categorical_logit_glm_lpmf.hpp       |  10 +-
 .../prim/neg_binomial_2_log_glm_lpmf.hpp      |  24 +-
 stan/math/opencl/prim/normal_id_glm_lpdf.hpp  |  18 +-
 stan/math/opencl/prim/num_elements.hpp        |   3 +-
 .../opencl/prim/ordered_logistic_glm_lpmf.hpp |  14 +-
 .../math/opencl/prim/poisson_log_glm_lpmf.hpp |  21 +-
 stan/math/opencl/prim/poisson_log_lpmf.hpp    |   8 +-
 stan/math/opencl/prim/poisson_lpmf.hpp        |   6 +-
 stan/math/opencl/prim/std_normal_lpdf.hpp     |   2 +-
 stan/math/opencl/zeros_strict_tri.hpp         |   6 +-
 .../prim/fun/offset_multiplier_constrain.hpp  | 264 ++++++++++++++++--
 stan/math/prim/fun/size_zero.hpp              |   3 +-
 stan/math/prim/functor/coupled_ode_system.hpp |  10 +-
 stan/math/prim/functor/ode_rk45.hpp           |   6 +-
 stan/math/prim/functor/reduce_sum.hpp         |   3 +-
 stan/math/prim/prob/bernoulli_lpmf.hpp        |   9 +-
 stan/math/prim/prob/beta_proportion_lpdf.hpp  |  22 +-
 stan/math/prim/prob/binomial_logit_lpmf.hpp   |  26 +-
 stan/math/prim/prob/cauchy_lpdf.hpp           |  22 +-
 stan/math/prim/prob/chi_square_cdf.hpp        |   4 +-
 stan/math/prim/prob/chi_square_lccdf.hpp      |   4 +-
 stan/math/prim/prob/chi_square_lcdf.hpp       |   4 +-
 stan/math/prim/prob/chi_square_lpdf.hpp       |  17 +-
 .../prim/prob/double_exponential_lpdf.hpp     |  18 +-
 stan/math/prim/prob/exp_mod_normal_lpdf.hpp   |  22 +-
 stan/math/prim/prob/exponential_lccdf.hpp     |  28 +-
 stan/math/prim/prob/exponential_lpdf.hpp      |  15 +-
 stan/math/prim/prob/frechet_lpdf.hpp          |  17 +-
 stan/math/prim/prob/gamma_cdf.hpp             |   4 +-
 stan/math/prim/prob/gamma_lccdf.hpp           |   4 +-
 stan/math/prim/prob/gamma_lcdf.hpp            |   4 +-
 stan/math/prim/prob/gamma_lpdf.hpp            |  17 +-
 stan/math/prim/prob/gumbel_lpdf.hpp           |  17 +-
 37 files changed, 456 insertions(+), 299 deletions(-)

diff --git a/stan/math/opencl/matrix_cl.hpp b/stan/math/opencl/matrix_cl.hpp
index ded00e3a821..6c5470e055a 100644
--- a/stan/math/opencl/matrix_cl.hpp
+++ b/stan/math/opencl/matrix_cl.hpp
@@ -11,7 +11,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/vec_concat.hpp>
-#include <CL/cl2.hpp>
+#include <CL/opencl.hpp>
 #include <algorithm>
 #include <iostream>
 #include <string>
@@ -31,7 +31,11 @@ namespace math {
  *  @{
  */
 
-template <typename, typename = void>
+// forward declare
+template <typename T>
+class arena_matrix_cl;
+
+template <typename>
 class matrix_cl;
 
 /**
@@ -39,7 +43,7 @@ class matrix_cl;
  * @tparam T an arithmetic type for the type stored in the OpenCL buffer.
  */
 template <typename T>
-class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
+class matrix_cl : public matrix_cl_base {
  private:
   cl::Buffer buffer_cl_;  // Holds the allocated memory on the device
   int rows_{0};           // Number of rows.
@@ -55,8 +59,6 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   // Forward declare the methods that work in place on the matrix
   template <matrix_cl_view matrix_view = matrix_cl_view::Entire>
   inline void zeros_strict_tri();
-  template <TriangularMapCL triangular_map = TriangularMapCL::LowerToUpper>
-  inline void triangular_transpose();
 
   int rows() const { return rows_; }
 
@@ -197,7 +199,9 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
     if (A.size() == 0) {
       return;
     }
-    initialize_buffer(A);
+    buffer_cl_ = cl::Buffer(opencl_context.context(), CL_MEM_READ_WRITE,
+                            sizeof(T) * this->size());
+    initialize_buffer_cl(A);
   }
 
   /**
@@ -212,6 +216,13 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
         write_events_(std::move(A.write_events_)),
         read_events_(std::move(A.read_events_)) {}
 
+  /**
+   * Constructor from `arena_matrix_cl`.
+   * @param A matrix_cl to move
+   */
+  // defined in rev/arena_matrix_cl.hpp
+  matrix_cl(const arena_matrix_cl<T>& A);  // NOLINT(runtime/explicit)
+
   /**
    * Constructor for the matrix_cl that creates a copy of a std::vector of Eigen
    * matrices on the OpenCL device. Each matrix is flattened into one column
@@ -272,7 +283,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   matrix_cl(const int rows, const int cols,
             matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(rows), cols_(cols), view_(partial_view) {
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     cl::Context& ctx = opencl_context.context();
@@ -310,7 +321,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(A.rows()), cols_(A.cols()), view_(partial_view) {
     using Mat_type = std::decay_t<ref_type_for_opencl_t<Mat>>;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer_no_heap_if<
@@ -415,8 +426,10 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
    * @tparam Expr type of the expression
    * @param expression expression
    */
+  // defined in kernel_generator/matrix_cl_conversion.hpp
   template <typename Expr,
-            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr>
+            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr,
+            require_not_matrix_cl_t<Expr>* = nullptr>
   matrix_cl(const Expr& expression);  // NOLINT(runtime/explicit)
 
   /**
@@ -438,13 +451,19 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
    */
   matrix_cl<T>& operator=(const matrix_cl<T>& a) {
     this->view_ = a.view();
-    this->rows_ = a.rows();
-    this->cols_ = a.cols();
     if (a.size() == 0) {
+      this->rows_ = a.rows();
+      this->cols_ = a.cols();
       return *this;
     }
     this->wait_for_read_write_events();
-    initialize_buffer(a);
+    if (this->size() != a.size()) {
+      buffer_cl_ = cl::Buffer(opencl_context.context(), CL_MEM_READ_WRITE,
+                              sizeof(T) * a.size());
+    }
+    this->rows_ = a.rows();
+    this->cols_ = a.cols();
+    initialize_buffer_cl(a);
     return *this;
   }
 
@@ -454,10 +473,20 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
    * @tparam Expr type of the expression
    * @param expression expression
    */
+  // defined in kernel_generator/matrix_cl_conversion.hpp
   template <typename Expr,
-            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr>
+            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr,
+            require_not_matrix_cl_t<Expr>* = nullptr>
   matrix_cl<T>& operator=(const Expr& expression);
 
+  /**
+   * Assignment of `arena_matrix_cl<T>`.
+   * @tparam Expr type of the expression
+   * @param expression expression
+   */
+  // defined in rev/arena_matrix_cl.hpp
+  matrix_cl<T>& operator=(const arena_matrix_cl<T>& other);
+
   /**
    * Evaluates `this`. This is a no-op.
    * @return `*this`
@@ -489,7 +518,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(const T* A) {
     cl::Event transfer_event;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -509,7 +538,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(T* A) {
     cl::Event transfer_event;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -548,7 +577,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
    */
   template <bool No_heap, typename U, std::enable_if_t<No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer(obj.data());
@@ -558,7 +587,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool No_heap, typename U, std::enable_if_t<!No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
     using U_val = std::decay_t<ref_type_for_opencl_t<U>>;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     auto* obj_heap = new U_val(std::move(obj));
@@ -582,15 +611,12 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
    * size of given matrix.
    * @param A matrix_cl
    */
-  void initialize_buffer(const matrix_cl<T>& A) {
-    cl::Context& ctx = opencl_context.context();
-    cl::CommandQueue queue = opencl_context.queue();
+  void initialize_buffer_cl(const matrix_cl<T>& A) {
     try {
-      buffer_cl_ = cl::Buffer(ctx, CL_MEM_READ_WRITE, sizeof(T) * this->size());
       cl::Event cstr_event;
-      queue.enqueueCopyBuffer(A.buffer(), this->buffer(), 0, 0,
-                              A.size() * sizeof(T), &A.write_events(),
-                              &cstr_event);
+      opencl_context.queue().enqueueCopyBuffer(A.buffer(), this->buffer(), 0, 0,
+                                               A.size() * sizeof(T),
+                                               &A.write_events(), &cstr_event);
       this->add_write_event(cstr_event);
       A.add_read_event(cstr_event);
     } catch (const cl::Error& e) {
@@ -621,13 +647,6 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
     delete static_cast<U*>(container);
   }
 };
-
-template <typename T>
-using matrix_cl_prim = matrix_cl<T, require_arithmetic_t<T>>;
-
-template <typename T>
-using matrix_cl_fp = matrix_cl<T, require_floating_point_t<T>>;
-
 /** @}*/
 
 }  // namespace math
diff --git a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
index 76a342408df..9de5422b293 100644
--- a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
@@ -66,13 +66,14 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
 
   if (N == 0) {
@@ -120,7 +121,7 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
       logp_expr, calc_if<need_theta_derivative>(theta_derivative_expr),
       calc_if<need_theta_derivative_sum>(colwise_sum(theta_derivative_expr)));
 
-  T_partials_return logp = sum(from_matrix_cl<Eigen::Dynamic, 1>(logp_cl));
+  T_partials_return logp = sum(from_matrix_cl(logp_cl));
   if (!std::isfinite(logp)) {
     results(check_cl(function, "Vector of dependent variables", y_val,
                      "in the interval [0, 1]"),
@@ -152,14 +153,13 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
     // transposition of a vector can be done without copying
     const matrix_cl<double> theta_derivative_transpose_cl(
         theta_derivative_cl.buffer(), 1, theta_derivative_cl.rows());
-    matrix_cl<double>& edge3_partials
-        = forward_as<matrix_cl<double>&>(ops_partials.edge3_.partials_);
     matrix_cl<double> edge3_partials_transpose_cl
         = theta_derivative_transpose_cl * x_val;
-    edge3_partials = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
-                                       edge3_partials_transpose_cl.cols(), 1);
+    ops_partials.edge3_.partials_
+        = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
+                            edge3_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      edge3_partials.add_write_event(
+      ops_partials.edge3_.partials_.add_write_event(
           edge3_partials_transpose_cl.write_events().back());
     }
   }
diff --git a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
index 8db546a1853..accc947ab17 100644
--- a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
@@ -36,7 +36,7 @@ return_type_t<T_prob_cl> bernoulli_logit_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/bernoulli_lpmf.hpp b/stan/math/opencl/prim/bernoulli_lpmf.hpp
index b6335842b60..ddbfb5e9074 100644
--- a/stan/math/opencl/prim/bernoulli_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_lpmf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_prob_cl> bernoulli_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
index df3ecc809f8..baf4974f60d 100644
--- a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
@@ -3,6 +3,7 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/opencl/prim/size.hpp>
+#include <stan/math/opencl/rev/arena_matrix_cl.hpp>
 #include <stan/math/opencl/rev/operands_and_partials.hpp>
 #include <stan/math/opencl/matrix_cl.hpp>
 #include <stan/math/opencl/copy.hpp>
@@ -60,10 +61,10 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
   static const char* function = "categorical_logit_glm_lpmf";
   if (is_y_vector) {
     check_size_match(function, "Rows of ", "x", N_instances, "size of ", "y",
-                     size(y));
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "beta", N_classes, "size of ",
-                   "alpha", size(alpha));
+                   "alpha", math::size(alpha));
   check_size_match(function, "Columns of ", "x", N_attributes, "Rows of",
                    "beta", beta.rows());
 
@@ -150,8 +151,9 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
     try {
       opencl_kernels::categorical_logit_glm_beta_derivative(
           cl::NDRange(local_size * N_attributes), cl::NDRange(local_size),
-          forward_as<matrix_cl<double>>(ops_partials.edge3_.partials_), temp,
-          y_val_cl, x_val, N_instances, N_attributes, N_classes, is_y_vector);
+          forward_as<arena_matrix_cl<double>>(ops_partials.edge3_.partials_),
+          temp, y_val_cl, x_val, N_instances, N_attributes, N_classes,
+          is_y_vector);
     } catch (const cl::Error& e) {
       check_opencl_error(function, e);
     }
diff --git a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
index 1f4aad75e88..1bfec8258b6 100644
--- a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
@@ -82,17 +82,18 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_phi_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "phi",
-                     size(phi));
+                     math::size(phi));
   }
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (N == 0) {
     return 0;
@@ -149,7 +150,7 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     check_opencl_error(function, e);
   }
 
-  T_partials_return logp = sum(from_matrix_cl<Dynamic, 1>(logp_cl));
+  T_partials_return logp = sum(from_matrix_cl(logp_cl));
   if (!std::isfinite(logp)) {
     results(
         check_cl(function, "Vector of dependent variables", y_val,
@@ -188,14 +189,13 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     // transposition of a vector can be done without copying
     const matrix_cl<double> theta_derivative_transpose_cl(
         theta_derivative_cl.buffer(), 1, theta_derivative_cl.rows());
-    matrix_cl<double>& edge3_partials
-        = forward_as<matrix_cl<double>&>(ops_partials.edge3_.partials_);
     matrix_cl<double> edge3_partials_transpose_cl
         = theta_derivative_transpose_cl * x_val;
-    edge3_partials = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
-                                       edge3_partials_transpose_cl.cols(), 1);
+    ops_partials.edge3_.partials_
+        = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
+                            edge3_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      edge3_partials.add_write_event(
+      ops_partials.edge3_.partials_.add_write_event(
           edge3_partials_transpose_cl.write_events().back());
     }
   }
@@ -205,7 +205,7 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     } else {
       forward_as<internal::broadcast_array<double>>(
           ops_partials.edge2_.partials_)[0]
-          = sum(from_matrix_cl<Dynamic, 1>(theta_derivative_sum_cl));
+          = sum(from_matrix_cl(theta_derivative_sum_cl));
     }
   }
   if (!is_constant_all<T_phi_cl>::value) {
@@ -214,7 +214,7 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     } else {
       forward_as<internal::broadcast_array<double>>(
           ops_partials.edge4_.partials_)[0]
-          = sum(from_matrix_cl<Dynamic, 1>(phi_derivative_cl));
+          = sum(from_matrix_cl(phi_derivative_cl));
     }
   }
   return ops_partials.build(logp);
diff --git a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
index 786a6f6cd2a..b999c25e535 100644
--- a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
+++ b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
@@ -74,17 +74,18 @@ normal_id_glm_lpdf(const T_y_cl& y, const T_x_cl& x, const T_alpha_cl& alpha,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_sigma_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "sigma",
-                     size(sigma));
+                     math::size(sigma));
   }
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (!include_summand<propto, T_y_cl, T_x_cl, T_alpha_cl, T_beta_cl,
                        T_sigma_cl>::value) {
@@ -171,14 +172,13 @@ normal_id_glm_lpdf(const T_y_cl& y, const T_x_cl& x, const T_alpha_cl& alpha,
     // transposition of a vector can be done without copying
     const matrix_cl<double> mu_derivative_transpose_cl(
         mu_derivative_cl.buffer(), 1, mu_derivative_cl.rows());
-    matrix_cl<double>& edge4_partials
-        = forward_as<matrix_cl<double>&>(ops_partials.edge4_.partials_);
     matrix_cl<double> edge4_partials_transpose_cl
         = mu_derivative_transpose_cl * x_val;
-    edge4_partials = matrix_cl<double>(edge4_partials_transpose_cl.buffer(),
-                                       edge4_partials_transpose_cl.cols(), 1);
+    ops_partials.edge4_.partials_
+        = matrix_cl<double>(edge4_partials_transpose_cl.buffer(),
+                            edge4_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      edge4_partials.add_write_event(
+      ops_partials.edge4_.partials_.add_write_event(
           edge4_partials_transpose_cl.write_events().back());
     }
   }
diff --git a/stan/math/opencl/prim/num_elements.hpp b/stan/math/opencl/prim/num_elements.hpp
index a5d965290d3..aedae05f616 100644
--- a/stan/math/opencl/prim/num_elements.hpp
+++ b/stan/math/opencl/prim/num_elements.hpp
@@ -3,6 +3,7 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/opencl/prim/size.hpp>
 
 namespace stan {
 namespace math {
@@ -16,7 +17,7 @@ namespace math {
 template <typename T,
           require_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr>
 size_t num_elements(const T& m) {
-  return size(m);
+  return math::size(m);
 }
 
 }  // namespace math
diff --git a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
index 71726353e90..5f5912b2328 100644
--- a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
@@ -63,19 +63,19 @@ return_type_t<T_x, T_beta, T_cuts> ordered_logistic_glm_lpmf(
 
   const size_t N_instances = x.rows();
   const size_t N_attributes = x.cols();
-  const size_t N_classes = size(cuts) + 1;
+  const size_t N_classes = math::size(cuts) + 1;
 
   if (is_y_vector) {
     check_size_match(function, "Rows of ", "x", N_instances, "rows of ", "y",
-                     size(y));
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x", N_attributes, "Size of",
-                   "beta", size(beta));
+                   "beta", math::size(beta));
 
   const auto& cuts_val = eval(value_of(cuts));
   if (N_classes >= 2) {
-    auto cuts_head = block_zero_based(cuts_val, 0, 0, size(cuts) - 1, 1);
-    auto cuts_tail = block_zero_based(cuts_val, 1, 0, size(cuts) - 1, 1);
+    auto cuts_head = block_zero_based(cuts_val, 0, 0, math::size(cuts) - 1, 1);
+    auto cuts_tail = block_zero_based(cuts_val, 1, 0, math::size(cuts) - 1, 1);
     check_cl(function, "Cuts", cuts_head, "ordered and finite")
         = cuts_head < cuts_tail && isfinite(cuts_head) && isfinite(cuts_tail);
   } else {
@@ -140,8 +140,8 @@ return_type_t<T_x, T_beta, T_cuts> ordered_logistic_glm_lpmf(
         edge2_partials_transpose.buffer(), edge2_partials_transpose.cols(),
         edge2_partials_transpose.rows());
     if (beta.rows() != 0) {
-      forward_as<matrix_cl<double>>(ops_partials.edge2_.partials_)
-          .add_write_event(edge2_partials_transpose.write_events().back());
+      ops_partials.edge2_.partials_.add_write_event(
+          edge2_partials_transpose.write_events().back());
     }
   }
   if (!is_constant_all<T_cuts>::value) {
diff --git a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
index 47ea9dc69ae..c6b730e7644 100644
--- a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
@@ -66,13 +66,14 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (N == 0) {
     return 0;
@@ -108,9 +109,8 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
   results(theta_derivative_cl, theta_derivative_sum_cl, logp_cl) = expressions(
       theta_derivative_expr, colwise_sum(theta_derivative_expr), logp_expr);
 
-  double theta_derivative_sum
-      = sum(from_matrix_cl<Dynamic, 1>(theta_derivative_sum_cl));
-  logp += sum(from_matrix_cl<Dynamic, 1>(logp_cl));
+  double theta_derivative_sum = sum(from_matrix_cl(theta_derivative_sum_cl));
+  logp += sum(from_matrix_cl(logp_cl));
   if (!std::isfinite(theta_derivative_sum)) {
     results(check_cl(function, "Vector of dependent variables", y_val,
                      "nonnegative"),
@@ -142,14 +142,13 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
     // transposition of a vector can be done without copying
     const matrix_cl<double> theta_derivative_transpose_cl(
         theta_derivative_cl.buffer(), 1, theta_derivative_cl.rows());
-    matrix_cl<double>& edge3_partials
-        = forward_as<matrix_cl<double>&>(ops_partials.edge3_.partials_);
     matrix_cl<double> edge3_partials_transpose_cl
         = theta_derivative_transpose_cl * x_val;
-    edge3_partials = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
-                                       edge3_partials_transpose_cl.cols(), 1);
+    ops_partials.edge3_.partials_
+        = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
+                            edge3_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      edge3_partials.add_write_event(
+      ops_partials.edge3_.partials_.add_write_event(
           edge3_partials_transpose_cl.write_events().back());
     }
   }
diff --git a/stan/math/opencl/prim/poisson_log_lpmf.hpp b/stan/math/opencl/prim/poisson_log_lpmf.hpp
index 515a311a690..31ed3dfaecd 100644
--- a/stan/math/opencl/prim/poisson_log_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_lpmf.hpp
@@ -39,7 +39,7 @@ return_type_t<T_log_rate_cl> poisson_log_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n, "Log rate parameter",
                          alpha);
-  const size_t N = is_n_vector ? size(n) : size(alpha);
+  const size_t N = is_n_vector ? math::size(n) : math::size(alpha);
   if (N == 0) {
     return 0.0;
   }
@@ -60,8 +60,8 @@ return_type_t<T_log_rate_cl> poisson_log_lpmf(const T_n_cl& n,
       = check_cl(function, "Log rate parameter", alpha_val, "not nan");
   auto alpha_not_nan = !isnan(alpha_val);
 
-  auto return_log_zero = colwise_max(
-      constant(0, N, 1) + (isinf(alpha_val) && (alpha_val > 0 || n != 0)));
+  auto return_log_zero
+      = colwise_max(cast<char>(isinf(alpha_val) && (alpha_val > 0 || n != 0)));
   auto exp_alpha = exp(alpha_val);
 
   auto logp1 = elt_multiply(n, alpha_val);
@@ -72,7 +72,7 @@ return_type_t<T_log_rate_cl> poisson_log_lpmf(const T_n_cl& n,
 
   auto deriv = n - exp_alpha;
 
-  matrix_cl<int> return_log_zero_cl;
+  matrix_cl<char> return_log_zero_cl;
   matrix_cl<double> logp_cl;
   matrix_cl<double> deriv_cl;
 
diff --git a/stan/math/opencl/prim/poisson_lpmf.hpp b/stan/math/opencl/prim/poisson_lpmf.hpp
index 896ef4d2209..2878e9d2849 100644
--- a/stan/math/opencl/prim/poisson_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_lpmf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_rate_cl> poisson_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n, "Rate parameter",
                          lambda);
-  const size_t N = is_n_vector ? size(n) : size(lambda);
+  const size_t N = is_n_vector ? math::size(n) : math::size(lambda);
   if (N == 0) {
     return 0.0;
   }
@@ -60,7 +60,7 @@ return_type_t<T_rate_cl> poisson_lpmf(const T_n_cl& n,
   auto lambda_nonnegative = 0.0 <= lambda_val;
 
   auto return_log_zero = colwise_max(
-      constant(0, N, 1) + (isinf(lambda_val) || (lambda_val == 0 && n != 0)));
+      cast<char>(isinf(lambda_val) || (lambda_val == 0 && n != 0)));
 
   auto logp1 = multiply_log(n, lambda_val);
   auto logp2 = static_select<include_summand<propto, T_rate_cl>::value>(
@@ -70,7 +70,7 @@ return_type_t<T_rate_cl> poisson_lpmf(const T_n_cl& n,
 
   auto deriv = elt_divide(n, lambda_val) - 1.0;
 
-  matrix_cl<int> return_log_zero_cl;
+  matrix_cl<char> return_log_zero_cl;
   matrix_cl<double> logp_cl;
   matrix_cl<double> deriv_cl;
 
diff --git a/stan/math/opencl/prim/std_normal_lpdf.hpp b/stan/math/opencl/prim/std_normal_lpdf.hpp
index b9782db683e..2bf4de436c1 100644
--- a/stan/math/opencl/prim/std_normal_lpdf.hpp
+++ b/stan/math/opencl/prim/std_normal_lpdf.hpp
@@ -35,7 +35,7 @@ inline return_type_t<T_y_cl> std_normal_lpdf(const T_y_cl& y) {
   using std::isfinite;
   using std::isnan;
 
-  const size_t N = size(y);
+  const size_t N = math::size(y);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/zeros_strict_tri.hpp b/stan/math/opencl/zeros_strict_tri.hpp
index 66ff85c64f4..3ff0b92842a 100644
--- a/stan/math/opencl/zeros_strict_tri.hpp
+++ b/stan/math/opencl/zeros_strict_tri.hpp
@@ -9,7 +9,7 @@
 #include <stan/math/opencl/kernels/fill_strict_tri.hpp>
 #include <stan/math/prim/meta.hpp>
 
-#include <CL/cl2.hpp>
+#include <CL/opencl.hpp>
 
 namespace stan {
 namespace math {
@@ -29,7 +29,7 @@ namespace math {
  */
 template <typename T>
 template <matrix_cl_view matrix_view>
-inline void matrix_cl<T, require_arithmetic_t<T>>::zeros_strict_tri() try {
+inline void matrix_cl<T>::zeros_strict_tri() try {
   if (matrix_view == matrix_cl_view::Entire) {
     invalid_argument(
         "zeros_strict_tri", "matrix_view",
@@ -40,7 +40,7 @@ inline void matrix_cl<T, require_arithmetic_t<T>>::zeros_strict_tri() try {
         "zeros_strict_tri", "matrix_view",
         "matrix_cl_view::Diagonal is not a valid template parameter value", "");
   }
-  if (size() == 0) {
+  if (this->size() == 0) {
     return;
   }
   this->view_ = both(this->view_, invert(matrix_view));
diff --git a/stan/math/prim/fun/offset_multiplier_constrain.hpp b/stan/math/prim/fun/offset_multiplier_constrain.hpp
index 5f8a984a949..85a58ce4771 100644
--- a/stan/math/prim/fun/offset_multiplier_constrain.hpp
+++ b/stan/math/prim/fun/offset_multiplier_constrain.hpp
@@ -7,6 +7,9 @@
 #include <stan/math/prim/fun/identity_constrain.hpp>
 #include <stan/math/prim/fun/multiply_log.hpp>
 #include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/sum.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <cmath>
 
 namespace stan {
@@ -35,19 +38,27 @@ namespace math {
  * @throw std::domain_error if sigma <= 0
  * @throw std::domain_error if mu is not finite
  */
-template <typename T, typename M, typename S>
-inline return_type_t<T, M, S> offset_multiplier_constrain(const T& x,
-                                                          const M& mu,
-                                                          const S& sigma) {
-  check_finite("offset_multiplier_constrain", "offset", mu);
-  if (sigma == 1) {
-    if (mu == 0) {
-      return identity_constrain(x);
-    }
-    return mu + x;
+template <typename T, typename M, typename S,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              T, M, S>* = nullptr>
+inline auto offset_multiplier_constrain(const T& x, const M& mu,
+                                        const S& sigma) {
+  const auto& mu_ref = to_ref(mu);
+  const auto& sigma_ref = to_ref(sigma);
+  if (is_matrix<T>::value && is_matrix<M>::value) {
+    check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
+  }
+  if (is_matrix<T>::value && is_matrix<S>::value) {
+    check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
+  } else if (is_matrix<M>::value && is_matrix<S>::value) {
+    check_matching_dims("offset_multiplier_constrain", "mu", mu, "sigma",
+                        sigma);
   }
-  check_positive_finite("offset_multiplier_constrain", "multiplier", sigma);
-  return fma(sigma, x, mu);
+
+  check_finite("offset_multiplier_constrain", "offset", value_of_rec(mu_ref));
+  check_positive_finite("offset_multiplier_constrain", "multiplier",
+                        value_of_rec(sigma_ref));
+  return fma(sigma_ref, x, mu_ref);
 }
 
 /**
@@ -76,22 +87,223 @@ inline return_type_t<T, M, S> offset_multiplier_constrain(const T& x,
  * @throw std::domain_error if sigma <= 0
  * @throw std::domain_error if mu is not finite
  */
+template <typename T, typename M, typename S,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              T, M, S>* = nullptr>
+inline auto offset_multiplier_constrain(const T& x, const M& mu, const S& sigma,
+                                        return_type_t<T, M, S>& lp) {
+  const auto& mu_ref = to_ref(mu);
+  const auto& sigma_ref = to_ref(sigma);
+  if (is_matrix<T>::value && is_matrix<M>::value) {
+    check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
+  }
+  if (is_matrix<T>::value && is_matrix<S>::value) {
+    check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
+  } else if (is_matrix<M>::value && is_matrix<S>::value) {
+    check_matching_dims("offset_multiplier_constrain", "mu", mu, "sigma",
+                        sigma);
+  }
+
+  check_finite("offset_multiplier_constrain", "offset", value_of_rec(mu_ref));
+  check_positive_finite("offset_multiplier_constrain", "multiplier",
+                        value_of_rec(sigma_ref));
+  if (math::size(sigma_ref) == 1) {
+    lp += sum(multiply_log(math::size(x), sigma_ref));
+  } else {
+    lp += sum(log(sigma_ref));
+  }
+  return fma(sigma_ref, x, mu_ref);
+}
+
+/**
+ * Overload for array of x and non-array mu and sigma
+ */
+template <typename T, typename M, typename S,
+          require_all_not_std_vector_t<M, S>* = nullptr>
+inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
+                                        const S& sigma) {
+  std::vector<
+      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu, sigma))>>
+      ret;
+  ret.reserve(x.size());
+  const auto& mu_ref = to_ref(mu);
+  const auto& sigma_ref = to_ref(sigma);
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma_ref));
+  }
+  return ret;
+}
+
+/**
+ * Overload for array of x and non-array mu and sigma with lp
+ */
+template <typename T, typename M, typename S,
+          require_all_not_std_vector_t<M, S>* = nullptr>
+inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
+                                        const S& sigma,
+                                        return_type_t<T, M, S>& lp) {
+  std::vector<
+      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu, sigma, lp))>>
+      ret;
+  ret.reserve(x.size());
+  const auto& mu_ref = to_ref(mu);
+  const auto& sigma_ref = to_ref(sigma);
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma_ref, lp));
+  }
+  return ret;
+}
+
+/**
+ * Overload for array of x and sigma and non-array mu
+ */
+template <typename T, typename M, typename S,
+          require_not_std_vector_t<M>* = nullptr>
+inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
+                                        const std::vector<S>& sigma) {
+  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
+  std::vector<
+      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu, sigma[0]))>>
+      ret;
+  ret.reserve(x.size());
+  const auto& mu_ref = to_ref(mu);
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma[i]));
+  }
+  return ret;
+}
+
+/**
+ * Overload for array of x and sigma and non-array mu with lp
+ */
+template <typename T, typename M, typename S,
+          require_not_std_vector_t<M>* = nullptr>
+inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
+                                        const std::vector<S>& sigma,
+                                        return_type_t<T, M, S>& lp) {
+  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
+  std::vector<plain_type_t<decltype(
+      offset_multiplier_constrain(x[0], mu, sigma[0], lp))>>
+      ret;
+  ret.reserve(x.size());
+  const auto& mu_ref = to_ref(mu);
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma[i], lp));
+  }
+  return ret;
+}
+
+/**
+ * Overload for array of x and mu and non-array sigma
+ */
+template <typename T, typename M, typename S,
+          require_not_std_vector_t<S>* = nullptr>
+inline auto offset_multiplier_constrain(const std::vector<T>& x,
+                                        const std::vector<M>& mu,
+                                        const S& sigma) {
+  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
+  std::vector<
+      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu[0], sigma))>>
+      ret;
+  ret.reserve(x.size());
+  const auto& sigma_ref = to_ref(sigma);
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma_ref));
+  }
+  return ret;
+}
+
+/**
+ * Overload for array of x and mu and non-array sigma with lp
+ */
+template <typename T, typename M, typename S,
+          require_not_std_vector_t<S>* = nullptr>
+inline auto offset_multiplier_constrain(const std::vector<T>& x,
+                                        const std::vector<M>& mu,
+                                        const S& sigma,
+                                        return_type_t<T, M, S>& lp) {
+  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
+  std::vector<plain_type_t<decltype(
+      offset_multiplier_constrain(x[0], mu[0], sigma, lp))>>
+      ret;
+  ret.reserve(x.size());
+  const auto& sigma_ref = to_ref(sigma);
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma_ref, lp));
+  }
+  return ret;
+}
+
+/**
+ * Overload for array of x, mu, and sigma
+ */
+template <typename T, typename M, typename S>
+inline auto offset_multiplier_constrain(const std::vector<T>& x,
+                                        const std::vector<M>& mu,
+                                        const std::vector<S>& sigma) {
+  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
+  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
+  std::vector<plain_type_t<decltype(
+      offset_multiplier_constrain(x[0], mu[0], sigma[0]))>>
+      ret;
+  ret.reserve(x.size());
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma[i]));
+  }
+  return ret;
+}
+
+/**
+ * Overload for array of x, mu, and sigma with lp
+ */
 template <typename T, typename M, typename S>
-inline return_type_t<T, M, S> offset_multiplier_constrain(const T& x,
-                                                          const M& mu,
-                                                          const S& sigma,
-                                                          T& lp) {
-  using std::log;
-  check_finite("offset_multiplier_constrain", "offset", mu);
-  if (sigma == 1) {
-    if (mu == 0) {
-      return identity_constrain(x);
-    }
-    return mu + x;
+inline auto offset_multiplier_constrain(const std::vector<T>& x,
+                                        const std::vector<M>& mu,
+                                        const std::vector<S>& sigma,
+                                        return_type_t<T, M, S>& lp) {
+  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
+  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
+  std::vector<plain_type_t<decltype(
+      offset_multiplier_constrain(x[0], mu[0], sigma[0], lp))>>
+      ret;
+  ret.reserve(x.size());
+  for (size_t i = 0; i < x.size(); ++i) {
+    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma[i], lp));
+  }
+  return ret;
+}
+
+/**
+ * Return the linearly transformed value for the specified unconstrained input
+ * and specified offset and multiplier. If the `Jacobian` parameter is `true`,
+ * the log density accumulator is incremented with the log absolute Jacobian
+ * determinant of the transform.  All of the transforms are specified with their
+ * Jacobians in the *Stan Reference Manual* chapter Constraint Transforms.
+ *
+ * @tparam Jacobian if `true`, increment log density accumulator with log
+ * absolute Jacobian determinant of constraining transform
+ * @tparam T A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
+ * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
+ * @tparam M A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
+ * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
+ * @tparam S A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
+ * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
+ * @param[in] x Unconstrained scalar input
+ * @param[in] mu offset of constrained output
+ * @param[in] sigma multiplier of constrained output
+ * @param[in, out] lp log density accumulator
+ * @return linear transformed value corresponding to inputs
+ * @throw std::domain_error if sigma <= 0
+ * @throw std::domain_error if mu is not finite
+ */
+template <bool Jacobian, typename T, typename M, typename S>
+inline auto offset_multiplier_constrain(const T& x, const M& mu, const S& sigma,
+                                        return_type_t<T, M, S>& lp) {
+  if (Jacobian) {
+    return offset_multiplier_constrain(x, mu, sigma, lp);
+  } else {
+    return offset_multiplier_constrain(x, mu, sigma);
   }
-  check_positive_finite("offset_multiplier_constrain", "multiplier", sigma);
-  lp += multiply_log(size(x), sigma);
-  return fma(sigma, x, mu);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/size_zero.hpp b/stan/math/prim/fun/size_zero.hpp
index 8af469eee9f..c6aa81da3fa 100644
--- a/stan/math/prim/fun/size_zero.hpp
+++ b/stan/math/prim/fun/size_zero.hpp
@@ -1,6 +1,7 @@
 #ifndef STAN_MATH_PRIM_FUN_SIZE_ZERO_HPP
 #define STAN_MATH_PRIM_FUN_SIZE_ZERO_HPP
 
+#include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <utility>
 
@@ -16,7 +17,7 @@ namespace math {
  */
 template <typename T>
 inline bool size_zero(const T& x) {
-  return !size(x);
+  return !math::size(x);
 }
 
 /**
diff --git a/stan/math/prim/functor/coupled_ode_system.hpp b/stan/math/prim/functor/coupled_ode_system.hpp
index a4dbe69603d..5298b53a570 100644
--- a/stan/math/prim/functor/coupled_ode_system.hpp
+++ b/stan/math/prim/functor/coupled_ode_system.hpp
@@ -67,9 +67,9 @@ struct coupled_ode_system_impl<true, F, T_y0, Args...> {
 
     dz_dt.resize(y.size());
 
-    Eigen::VectorXd f_y_t
-        = apply([&](const Args&... args) { return f_(t, y, msgs_, args...); },
-                args_tuple_);
+    Eigen::VectorXd f_y_t = math::apply(
+        [&](const Args&... args) { return f_(t, y, msgs_, args...); },
+        args_tuple_);
 
     check_size_match("coupled_ode_system", "dy_dt", f_y_t.size(), "states",
                      y.size());
@@ -110,8 +110,8 @@ struct coupled_ode_system
                      const Eigen::Matrix<T_y0, Eigen::Dynamic, 1>& y0,
                      std::ostream* msgs, const Args&... args)
       : coupled_ode_system_impl<
-            std::is_arithmetic<return_type_t<T_y0, Args...>>::value, F, T_y0,
-            Args...>(f, y0, msgs, args...) {}
+          std::is_arithmetic<return_type_t<T_y0, Args...>>::value, F, T_y0,
+          Args...>(f, y0, msgs, args...) {}
 };
 
 }  // namespace math
diff --git a/stan/math/prim/functor/ode_rk45.hpp b/stan/math/prim/functor/ode_rk45.hpp
index 8a90bd1e2c7..4850b1169c6 100644
--- a/stan/math/prim/functor/ode_rk45.hpp
+++ b/stan/math/prim/functor/ode_rk45.hpp
@@ -79,7 +79,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   std::tuple<ref_type_t<Args>...> args_ref_tuple(args...);
 
-  apply(
+  math::apply(
       [&](const auto&... args_ref) {
         // Code from https://stackoverflow.com/a/17340003
         std::vector<int> unused_temp{
@@ -102,7 +102,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   using return_t = return_type_t<T_y0, T_t0, T_ts, Args...>;
   // creates basic or coupled system by template specializations
-  auto&& coupled_system = apply(
+  auto&& coupled_system = math::apply(
       [&](const auto&... args_ref) {
         return coupled_ode_system<F, T_y0_t0, ref_type_t<Args>...>(f, y0, msgs,
                                                                    args_ref...);
@@ -128,7 +128,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
       observer_initial_recorded = true;
       return;
     }
-    apply(
+    math::apply(
         [&](const auto&... args_ref) {
           y.emplace_back(ode_store_sensitivities(
               f, coupled_state, y0, t0, ts[time_index], msgs, args_ref...));
diff --git a/stan/math/prim/functor/reduce_sum.hpp b/stan/math/prim/functor/reduce_sum.hpp
index 868de060f32..ff0b4f5d592 100644
--- a/stan/math/prim/functor/reduce_sum.hpp
+++ b/stan/math/prim/functor/reduce_sum.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 
 #include <tbb/task_arena.h>
 #include <tbb/parallel_reduce.h>
@@ -80,7 +81,7 @@ struct reduce_sum_impl<ReduceFunction, require_arithmetic_t<ReturnType>,
         sub_slice.emplace_back(vmapped_[i]);
       }
 
-      sum_ += apply(
+      sum_ += math::apply(
           [&](auto&&... args) {
             return ReduceFunction()(sub_slice, r.begin(), r.end() - 1, &msgs_,
                                     args...);
diff --git a/stan/math/prim/prob/bernoulli_lpmf.hpp b/stan/math/prim/prob/bernoulli_lpmf.hpp
index 4e5876f7359..9b54f02d839 100644
--- a/stan/math/prim/prob/bernoulli_lpmf.hpp
+++ b/stan/math/prim/prob/bernoulli_lpmf.hpp
@@ -42,7 +42,8 @@ return_type_t<T_prob> bernoulli_lpmf(const T_n& n, const T_prob& theta) {
   const T_n_ref n_ref = to_ref(n);
   const T_theta_ref theta_ref = to_ref(theta);
   check_bounded(function, "n", n_ref, 0, 1);
-  check_bounded(function, "Probability parameter", theta_ref, 0.0, 1.0);
+  check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0,
+                1.0);
 
   if (size_zero(n, theta)) {
     return 0.0;
@@ -58,7 +59,7 @@ return_type_t<T_prob> bernoulli_lpmf(const T_n& n, const T_prob& theta) {
   scalar_seq_view<T_theta_ref> theta_vec(theta_ref);
   size_t N = max_size(n, theta);
 
-  if (size(theta) == 1) {
+  if (math::size(theta) == 1) {
     size_t sum = 0;
     for (size_t n = 0; n < N; n++) {
       sum += n_vec.val(n);
@@ -83,8 +84,8 @@ return_type_t<T_prob> bernoulli_lpmf(const T_n& n, const T_prob& theta) {
       logp += (N - sum) * log1m_theta;
 
       if (!is_constant_all<T_prob>::value) {
-        ops_partials.edge1_.partials_[0] += sum / theta_dbl;
-        ops_partials.edge1_.partials_[0] += (N - sum) / (theta_dbl - 1);
+        ops_partials.edge1_.partials_[0] += sum * inv(theta_dbl);
+        ops_partials.edge1_.partials_[0] += (N - sum) * inv(theta_dbl - 1);
       }
     }
   } else {
diff --git a/stan/math/prim/prob/beta_proportion_lpdf.hpp b/stan/math/prim/prob/beta_proportion_lpdf.hpp
index f3bf5ff409b..523fdc1c138 100644
--- a/stan/math/prim/prob/beta_proportion_lpdf.hpp
+++ b/stan/math/prim/prob/beta_proportion_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
@@ -51,8 +52,6 @@ return_type_t<T_y, T_loc, T_prec> beta_proportion_lpdf(const T_y& y,
                                                        const T_loc& mu,
                                                        const T_prec& kappa) {
   using T_partials_return = partials_return_t<T_y, T_loc, T_prec>;
-  using T_partials_return_kappa = return_type_t<T_prec>;
-  using T_partials_array = Eigen::Array<T_partials_return, Eigen::Dynamic, 1>;
   using std::log;
   using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
   using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
@@ -68,23 +67,14 @@ return_type_t<T_y, T_loc, T_prec> beta_proportion_lpdf(const T_y& y,
   T_mu_ref mu_ref = mu;
   T_kappa_ref kappa_ref = kappa;
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
-  const auto& kappa_col = as_column_vector_or_scalar(kappa_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& mu_arr = as_array_or_scalar(mu_col);
-  const auto& kappa_arr
-      = promote_scalar<T_partials_return_kappa>(as_array_or_scalar(kappa_col));
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
-  ref_type_t<decltype(value_of(kappa_arr))> kappa_val = value_of(kappa_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
+  decltype(auto) kappa_val = to_ref(as_value_column_array_or_scalar(kappa_ref));
 
   check_positive(function, "Location parameter", mu_val);
   check_less(function, "Location parameter", mu_val, 1.0);
   check_positive_finite(function, "Precision parameter", kappa_val);
-  check_bounded(function, "Random variable", y_val, 0, 1);
+  check_bounded(function, "Random variable", value_of(y_val), 0, 1);
 
   if (!include_summand<propto, T_y, T_loc, T_prec>::value) {
     return 0;
@@ -99,7 +89,7 @@ return_type_t<T_y, T_loc, T_prec> beta_proportion_lpdf(const T_y& y,
   size_t N = max_size(y, mu, kappa);
   T_partials_return logp(0);
   if (include_summand<propto, T_prec>::value) {
-    logp += sum(lgamma(kappa_val)) * N / size(kappa);
+    logp += sum(lgamma(kappa_val)) * N / math::size(kappa);
   }
   if (include_summand<propto, T_loc, T_prec>::value) {
     logp -= sum(lgamma(mukappa) + lgamma(kappa_val - mukappa)) * N
diff --git a/stan/math/prim/prob/binomial_logit_lpmf.hpp b/stan/math/prim/prob/binomial_logit_lpmf.hpp
index 3c4993b9495..08e42e2873a 100644
--- a/stan/math/prim/prob/binomial_logit_lpmf.hpp
+++ b/stan/math/prim/prob/binomial_logit_lpmf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/binomial_coefficient_log.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
@@ -33,7 +34,9 @@ namespace math {
  * @throw std::domain_error if N is negative or probability parameter is invalid
  * @throw std::invalid_argument if vector sizes do not match
  */
-template <bool propto, typename T_n, typename T_N, typename T_prob>
+template <bool propto, typename T_n, typename T_N, typename T_prob,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              T_n, T_N, T_prob>* = nullptr>
 return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
                                           const T_prob& alpha) {
   using T_partials_return = partials_return_t<T_n, T_N, T_prob>;
@@ -49,19 +52,11 @@ return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
   T_N_ref N_ref = N;
   T_alpha_ref alpha_ref = alpha;
 
-  const auto& n_col = as_column_vector_or_scalar(n_ref);
-  const auto& N_col = as_column_vector_or_scalar(N_ref);
-  const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
+  decltype(auto) n_val = to_ref(as_value_column_array_or_scalar(n_ref));
+  decltype(auto) N_val = to_ref(as_value_column_array_or_scalar(N_ref));
+  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
 
-  const auto& n_arr = as_array_or_scalar(n_col);
-  const auto& N_arr = as_array_or_scalar(N_col);
-  const auto& alpha_arr = as_array_or_scalar(alpha_col);
-
-  ref_type_t<decltype(value_of(n_arr))> n_val = value_of(n_arr);
-  ref_type_t<decltype(value_of(N_arr))> N_val = value_of(N_arr);
-  ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
-
-  check_bounded(function, "Successes variable", n_val, 0, N_val);
+  check_bounded(function, "Successes variable", value_of(n_val), 0, N_val);
   check_nonnegative(function, "Population size parameter", N_val);
   check_finite(function, "Probability parameter", alpha_val);
 
@@ -92,10 +87,11 @@ return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
       ops_partials.edge1_.partials_
           = n_val * inv_logit_neg_alpha - (N_val - n_val) * inv_logit_alpha;
     } else {
-      T_partials_return sum_n = sum(n_val) * maximum_size / size(n);
+      T_partials_return sum_n = sum(n_val) * maximum_size / math::size(n);
       ops_partials.edge1_.partials_[0] = forward_as<T_partials_return>(
           sum_n * inv_logit_neg_alpha
-          - (sum(N_val) * maximum_size / size(N) - sum_n) * inv_logit_alpha);
+          - (sum(N_val) * maximum_size / math::size(N) - sum_n)
+                * inv_logit_alpha);
     }
   }
 
diff --git a/stan/math/prim/prob/cauchy_lpdf.hpp b/stan/math/prim/prob/cauchy_lpdf.hpp
index d1da16262b8..6e16581e786 100644
--- a/stan/math/prim/prob/cauchy_lpdf.hpp
+++ b/stan/math/prim/prob/cauchy_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/log1p.hpp>
@@ -42,7 +43,6 @@ template <bool propto, typename T_y, typename T_loc, typename T_scale,
 return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
                                                const T_scale& sigma) {
   using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
-  using T_partials_array = Eigen::Array<T_partials_return, Eigen::Dynamic, 1>;
   using std::log;
   using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
   using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
@@ -65,17 +65,9 @@ return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
   operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref> ops_partials(
       y_ref, mu_ref, sigma_ref);
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
-  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& mu_arr = as_array_or_scalar(mu_col);
-  const auto& sigma_arr = as_array_or_scalar(sigma_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
-  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
+  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
   check_not_nan(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
   check_positive_finite(function, "Scale parameter", sigma_val);
@@ -92,7 +84,7 @@ return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
     logp -= N * LOG_PI;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
@@ -101,8 +93,8 @@ return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
     const auto& y_minus_mu_squared
         = to_ref_if<!is_constant_all<T_scale>::value>(square(y_minus_mu));
     if (!is_constant_all<T_y, T_loc>::value) {
-      const auto& mu_deriv = to_ref_if<(!is_constant_all<T_y>::value
-                                        && !is_constant_all<T_loc>::value)>(
+      auto mu_deriv = to_ref_if<(!is_constant_all<T_y>::value
+                                 && !is_constant_all<T_loc>::value)>(
           2 * y_minus_mu / (sigma_squared + y_minus_mu_squared));
       if (!is_constant_all<T_y>::value) {
         if (is_vector<T_y>::value) {
diff --git a/stan/math/prim/prob/chi_square_cdf.hpp b/stan/math/prim/prob/chi_square_cdf.hpp
index ea439f86bf6..8e8f6cf5148 100644
--- a/stan/math/prim/prob/chi_square_cdf.hpp
+++ b/stan/math/prim/prob/chi_square_cdf.hpp
@@ -69,9 +69,9 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lccdf.hpp b/stan/math/prim/prob/chi_square_lccdf.hpp
index 8347406a373..cb6bead9799 100644
--- a/stan/math/prim/prob/chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/chi_square_lccdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> chi_square_lccdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lcdf.hpp b/stan/math/prim/prob/chi_square_lcdf.hpp
index dc7c196a454..5fefad3fe3a 100644
--- a/stan/math/prim/prob/chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/chi_square_lcdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> chi_square_lcdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lpdf.hpp b/stan/math/prim/prob/chi_square_lpdf.hpp
index 7bac96d7e19..885d9a3c988 100644
--- a/stan/math/prim/prob/chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/chi_square_lpdf.hpp
@@ -4,7 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
-#include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
@@ -54,17 +54,10 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {
                          "Degrees of freedom parameter", nu);
   T_y_ref y_ref = y;
   T_nu_ref nu_ref = nu;
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& nu_col = as_column_vector_or_scalar(nu_ref);
 
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& nu_arr = as_array_or_scalar(nu_col);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) nu_val = to_ref(as_value_column_array_or_scalar(nu_ref));
 
-  ref_type_if_t<include_summand<propto, T_y>::value, decltype(value_of(y_arr))>
-      y_val = value_of(y_arr);
-  ref_type_if_t<include_summand<propto, T_dof>::value,
-                decltype(value_of(nu_arr))>
-      nu_val = value_of(nu_arr);
   check_nonnegative(function, "Random variable", y_val);
   check_positive_finite(function, "Degrees of freedom parameter", nu_val);
 
@@ -81,12 +74,12 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {
 
   T_partials_return logp(0);
   if (include_summand<propto, T_dof>::value) {
-    logp -= sum(nu_val * HALF_LOG_TWO + lgamma(half_nu)) * N / size(nu);
+    logp -= sum(nu_val * HALF_LOG_TWO + lgamma(half_nu)) * N / math::size(nu);
   }
   logp += sum((half_nu - 1.0) * log_y);
 
   if (include_summand<propto, T_y>::value) {
-    logp -= 0.5 * sum(y_val) * N / size(y);
+    logp -= 0.5 * sum(y_val) * N / math::size(y);
   }
 
   operands_and_partials<T_y_ref, T_nu_ref> ops_partials(y_ref, nu_ref);
diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp
index 2587a1ed86f..3db919362d2 100644
--- a/stan/math/prim/prob/double_exponential_lpdf.hpp
+++ b/stan/math/prim/prob/double_exponential_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/fabs.hpp>
 #include <stan/math/prim/fun/inv.hpp>
@@ -41,7 +42,6 @@ template <bool propto, typename T_y, typename T_loc, typename T_scale,
 return_type_t<T_y, T_loc, T_scale> double_exponential_lpdf(
     const T_y& y, const T_loc& mu, const T_scale& sigma) {
   using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
-  using T_partials_array = Eigen::Array<T_partials_return, Eigen::Dynamic, 1>;
   using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
   using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
   using T_sigma_ref = ref_type_if_t<!is_constant<T_scale>::value, T_scale>;
@@ -63,17 +63,9 @@ return_type_t<T_y, T_loc, T_scale> double_exponential_lpdf(
   operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref> ops_partials(
       y_ref, mu_ref, sigma_ref);
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
-  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& mu_arr = as_array_or_scalar(mu_col);
-  const auto& sigma_arr = as_array_or_scalar(sigma_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
-  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
+  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
 
   check_finite(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
@@ -91,7 +83,7 @@ return_type_t<T_y, T_loc, T_scale> double_exponential_lpdf(
     logp -= N * LOG_TWO;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   logp -= sum(scaled_diff);
 
diff --git a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
index d92de214ae5..ba3eddbcbde 100644
--- a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/erfc.hpp>
 #include <stan/math/prim/fun/exp.hpp>
@@ -43,20 +44,11 @@ return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_lpdf(
   T_sigma_ref sigma_ref = sigma;
   T_lambda_ref lambda_ref = lambda;
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
-  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
-  const auto& lambda_col = as_column_vector_or_scalar(lambda_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& mu_arr = as_array_or_scalar(mu_col);
-  const auto& sigma_arr = as_array_or_scalar(sigma_col);
-  const auto& lambda_arr = as_array_or_scalar(lambda_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
-  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
-  ref_type_t<decltype(value_of(lambda_arr))> lambda_val = value_of(lambda_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
+  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
+  decltype(auto) lambda_val
+      = to_ref(as_value_column_array_or_scalar(lambda_ref));
 
   check_not_nan(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
@@ -87,7 +79,7 @@ return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_lpdf(
     logp -= LOG_TWO * N;
   }
   if (include_summand<propto, T_inv_scale>::value) {
-    logp += sum(log(lambda_val)) * N / size(lambda);
+    logp += sum(log(lambda_val)) * N / math::size(lambda);
   }
   const auto& log_erfc_calc = log(erfc_calc);
   logp
diff --git a/stan/math/prim/prob/exponential_lccdf.hpp b/stan/math/prim/prob/exponential_lccdf.hpp
index f7d74bea979..c48c90ab9d3 100644
--- a/stan/math/prim/prob/exponential_lccdf.hpp
+++ b/stan/math/prim/prob/exponential_lccdf.hpp
@@ -3,8 +3,8 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
-#include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
@@ -14,7 +14,9 @@
 namespace stan {
 namespace math {
 
-template <typename T_y, typename T_inv_scale>
+template <typename T_y, typename T_inv_scale,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              T_y, T_inv_scale>* = nullptr>
 return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
                                                   const T_inv_scale& beta) {
   using T_partials_return = partials_return_t<T_y, T_inv_scale>;
@@ -26,14 +28,8 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
   T_y_ref y_ref = y;
   T_beta_ref beta_ref = beta;
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& beta_arr = as_array_or_scalar(beta_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
+  auto y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  auto beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
 
   check_nonnegative(function, "Random variable", y_val);
   check_positive_finite(function, "Inverse scale parameter", beta_val);
@@ -51,13 +47,11 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     using beta_val_array = Eigen::Array<beta_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_y>::value && !is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = T_partials_array::Constant(
-          size(y), -forward_as<beta_val_scalar>(beta_val));
+          math::size(y), -forward_as<beta_val_scalar>(beta_val));
     } else if (is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = -forward_as<beta_val_array>(beta_val);
     } else {
-      forward_as<internal::broadcast_array<T_partials_return>>(
-          ops_partials.edge1_.partials_)
-          = -forward_as<beta_val_scalar>(beta_val);
+      ops_partials.edge1_.partials_[0] = -sum(beta_val);
     }
   }
   if (!is_constant_all<T_inv_scale>::value) {
@@ -65,13 +59,11 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     using y_val_array = Eigen::Array<y_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_inv_scale>::value && !is_vector<T_y>::value) {
       ops_partials.edge2_.partials_ = T_partials_array::Constant(
-          size(beta), -forward_as<y_val_scalar>(y_val));
+          math::size(beta), -forward_as<y_val_scalar>(y_val));
     } else if (is_vector<T_y>::value) {
       ops_partials.edge2_.partials_ = -forward_as<y_val_array>(y_val);
     } else {
-      forward_as<internal::broadcast_array<T_partials_return>>(
-          ops_partials.edge2_.partials_)
-          = -forward_as<y_val_scalar>(y_val);
+      ops_partials.edge2_.partials_[0] = -sum(y_val);
     }
   }
   return ops_partials.build(ccdf_log);
diff --git a/stan/math/prim/prob/exponential_lpdf.hpp b/stan/math/prim/prob/exponential_lpdf.hpp
index 6f75b1a31ae..18f3bbfd0da 100644
--- a/stan/math/prim/prob/exponential_lpdf.hpp
+++ b/stan/math/prim/prob/exponential_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
@@ -61,14 +62,8 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
   T_y_ref y_ref = y;
   T_beta_ref beta_ref = beta;
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& beta_arr = as_array_or_scalar(beta_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
 
   check_nonnegative(function, "Random variable", y_val);
   check_positive_finite(function, "Inverse scale parameter", beta_val);
@@ -81,7 +76,7 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
 
   T_partials_return logp(0.0);
   if (include_summand<propto, T_inv_scale>::value) {
-    logp = sum(log(beta_val)) * max_size(y, beta) / size(beta);
+    logp = sum(log(beta_val)) * max_size(y, beta) / math::size(beta);
   }
   if (include_summand<propto, T_y, T_inv_scale>::value) {
     logp -= sum(beta_val * y_val);
@@ -92,7 +87,7 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
     using beta_val_array = Eigen::Array<beta_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_y>::value && !is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = T_partials_array::Constant(
-          size(y), -forward_as<beta_val_scalar>(beta_val));
+          math::size(y), -forward_as<beta_val_scalar>(beta_val));
     } else if (is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = -forward_as<beta_val_array>(beta_val);
     } else {
diff --git a/stan/math/prim/prob/frechet_lpdf.hpp b/stan/math/prim/prob/frechet_lpdf.hpp
index a6658972221..b305f0a8f20 100644
--- a/stan/math/prim/prob/frechet_lpdf.hpp
+++ b/stan/math/prim/prob/frechet_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
@@ -42,17 +43,9 @@ return_type_t<T_y, T_shape, T_scale> frechet_lpdf(const T_y& y,
   T_sigma_ref sigma_ref = sigma;
   using std::pow;
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
-  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& alpha_arr = as_array_or_scalar(alpha_col);
-  const auto& sigma_arr = as_array_or_scalar(sigma_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
-  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
+  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
 
   check_positive(function, "Random variable", y_val);
   check_positive_finite(function, "Shape parameter", alpha_val);
@@ -79,7 +72,7 @@ return_type_t<T_y, T_shape, T_scale> frechet_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, sigma);
   T_partials_return logp = -sum(sigma_div_y_pow_alpha);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp -= sum((alpha_val + 1.0) * log_y) * N / max_size(y, alpha);
diff --git a/stan/math/prim/prob/gamma_cdf.hpp b/stan/math/prim/prob/gamma_cdf.hpp
index 2918ebb9f22..fb9f88f788c 100644
--- a/stan/math/prim/prob/gamma_cdf.hpp
+++ b/stan/math/prim/prob/gamma_cdf.hpp
@@ -81,9 +81,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_cdf(const T_y& y,
   using std::pow;
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp
index 12fbd24a9f2..c031c1a7e63 100644
--- a/stan/math/prim/prob/gamma_lccdf.hpp
+++ b/stan/math/prim/prob/gamma_lccdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lcdf.hpp b/stan/math/prim/prob/gamma_lcdf.hpp
index b40c75bc7f3..19c15050873 100644
--- a/stan/math/prim/prob/gamma_lcdf.hpp
+++ b/stan/math/prim/prob/gamma_lcdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lpdf.hpp b/stan/math/prim/prob/gamma_lpdf.hpp
index e23ed36ab8f..50e7a6b65c8 100644
--- a/stan/math/prim/prob/gamma_lpdf.hpp
+++ b/stan/math/prim/prob/gamma_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
@@ -63,17 +64,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lpdf(const T_y& y,
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
-  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& alpha_arr = as_array_or_scalar(alpha_col);
-  const auto& beta_arr = as_array_or_scalar(beta_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
-  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
+  decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
 
   check_not_nan(function, "Random variable", y_val);
   check_positive_finite(function, "Shape parameter", alpha_val);
@@ -99,7 +92,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, beta);
   T_partials_return logp(0.0);
   if (include_summand<propto, T_shape>::value) {
-    logp = -sum(lgamma(alpha_val)) * N / size(alpha);
+    logp = -sum(lgamma(alpha_val)) * N / math::size(alpha);
   }
   const auto& log_y = to_ref_if<is_constant_all<T_shape>::value>(log(y_val));
   if (include_summand<propto, T_shape, T_inv_scale>::value) {
diff --git a/stan/math/prim/prob/gumbel_lpdf.hpp b/stan/math/prim/prob/gumbel_lpdf.hpp
index a0f97e7d111..22814b8982b 100644
--- a/stan/math/prim/prob/gumbel_lpdf.hpp
+++ b/stan/math/prim/prob/gumbel_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
@@ -49,17 +50,9 @@ return_type_t<T_y, T_loc, T_scale> gumbel_lpdf(const T_y& y, const T_loc& mu,
   T_mu_ref mu_ref = mu;
   T_beta_ref beta_ref = beta;
 
-  const auto& y_col = as_column_vector_or_scalar(y_ref);
-  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
-  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
-
-  const auto& y_arr = as_array_or_scalar(y_col);
-  const auto& mu_arr = as_array_or_scalar(mu_col);
-  const auto& beta_arr = as_array_or_scalar(beta_col);
-
-  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
-  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
-  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
+  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
+  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
+  decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
 
   check_not_nan(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
@@ -85,7 +78,7 @@ return_type_t<T_y, T_loc, T_scale> gumbel_lpdf(const T_y& y, const T_loc& mu,
   size_t N = max_size(y, mu, beta);
   T_partials_return logp = -sum(y_minus_mu_over_beta + exp_y_m_mu_over_beta);
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(beta_val)) * N / size(beta);
+    logp -= sum(log(beta_val)) * N / math::size(beta);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {

From 48e2301ca07fa49134079b1851dc9f9208aff0bc Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 15:48:58 +0200
Subject: [PATCH 2/9] Backport batch 2

---
 stan/math/prim/prob/inv_chi_square_cdf.hpp         |  4 ++--
 stan/math/prim/prob/inv_chi_square_lccdf.hpp       |  4 ++--
 stan/math/prim/prob/inv_chi_square_lpdf.hpp        |  4 ++--
 stan/math/prim/prob/inv_gamma_cdf.hpp              |  4 ++--
 stan/math/prim/prob/inv_gamma_lccdf.hpp            |  4 ++--
 stan/math/prim/prob/inv_gamma_lcdf.hpp             |  4 ++--
 stan/math/prim/prob/inv_gamma_lpdf.hpp             |  2 +-
 stan/math/prim/prob/logistic_lpdf.hpp              |  2 +-
 stan/math/prim/prob/lognormal_lpdf.hpp             |  4 ++--
 stan/math/prim/prob/normal_lpdf.hpp                |  2 +-
 stan/math/prim/prob/normal_sufficient_lpdf.hpp     |  4 ++--
 stan/math/prim/prob/ordered_logistic_lpmf.hpp      |  2 +-
 stan/math/prim/prob/ordered_probit_lpmf.hpp        |  2 +-
 stan/math/prim/prob/pareto_lpdf.hpp                |  2 +-
 stan/math/prim/prob/pareto_type_2_lpdf.hpp         |  4 ++--
 stan/math/prim/prob/poisson_log_lpmf.hpp           |  4 ++--
 stan/math/prim/prob/poisson_lpmf.hpp               |  4 ++--
 stan/math/prim/prob/rayleigh_lpdf.hpp              |  4 ++--
 stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp  |  4 ++--
 .../math/prim/prob/scaled_inv_chi_square_lccdf.hpp |  4 ++--
 stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp |  4 ++--
 stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp | 14 +++++++-------
 stan/math/prim/prob/skew_normal_lpdf.hpp           |  2 +-
 stan/math/prim/prob/student_t_cdf.hpp              |  6 +++---
 stan/math/prim/prob/student_t_lccdf.hpp            |  6 +++---
 stan/math/prim/prob/student_t_lcdf.hpp             |  6 +++---
 stan/math/prim/prob/student_t_lpdf.hpp             |  4 ++--
 stan/math/prim/prob/uniform_lcdf.hpp               |  4 ++--
 stan/math/prim/prob/uniform_lpdf.hpp               |  4 ++--
 stan/math/prim/prob/von_mises_lpdf.hpp             |  2 +-
 stan/math/prim/prob/weibull_lpdf.hpp               |  2 +-
 stan/math/rev/functor/coupled_ode_system.hpp       |  7 ++++---
 stan/math/rev/functor/cvodes_integrator.hpp        |  9 +++++----
 stan/math/rev/functor/ode_adams.hpp                |  3 ++-
 stan/math/rev/functor/ode_bdf.hpp                  |  3 ++-
 stan/math/rev/functor/reduce_sum.hpp               |  6 +++---
 36 files changed, 77 insertions(+), 73 deletions(-)

diff --git a/stan/math/prim/prob/inv_chi_square_cdf.hpp b/stan/math/prim/prob/inv_chi_square_cdf.hpp
index 1f1ff29c386..429dbfadcc7 100644
--- a/stan/math/prim/prob/inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_cdf.hpp
@@ -69,9 +69,9 @@ return_type_t<T_y, T_dof> inv_chi_square_cdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/inv_chi_square_lccdf.hpp b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
index 394e7956f97..e7c6a1e9c88 100644
--- a/stan/math/prim/prob/inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> inv_chi_square_lccdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/inv_chi_square_lpdf.hpp b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
index 68ed6bfb274..1b41ceda34a 100644
--- a/stan/math/prim/prob/inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
@@ -87,11 +87,11 @@ return_type_t<T_y, T_dof> inv_chi_square_lpdf(const T_y& y, const T_dof& nu) {
   size_t N = max_size(y, nu);
   T_partials_return logp = -sum((half_nu + 1.0) * log_y);
   if (include_summand<propto, T_dof>::value) {
-    logp -= (sum(nu_val) * HALF_LOG_TWO + sum(lgamma(half_nu))) * N / size(nu);
+    logp -= (sum(nu_val) * HALF_LOG_TWO + sum(lgamma(half_nu))) * N / math::size(nu);
   }
   if (include_summand<propto, T_y>::value) {
     const auto& inv_y = to_ref_if<!is_constant_all<T_y>::value>(inv(y_val));
-    logp -= 0.5 * sum(inv_y) * N / size(y);
+    logp -= 0.5 * sum(inv_y) * N / math::size(y);
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_ = (0.5 * inv_y - half_nu - 1.0) * inv_y;
     }
diff --git a/stan/math/prim/prob/inv_gamma_cdf.hpp b/stan/math/prim/prob/inv_gamma_cdf.hpp
index 89271caa42d..3ac21f757d4 100644
--- a/stan/math/prim/prob/inv_gamma_cdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_cdf.hpp
@@ -80,9 +80,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_cdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lccdf.hpp b/stan/math/prim/prob/inv_gamma_lccdf.hpp
index 22212e5adff..c731eccd0e3 100644
--- a/stan/math/prim/prob/inv_gamma_lccdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lccdf.hpp
@@ -66,9 +66,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lccdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lcdf.hpp b/stan/math/prim/prob/inv_gamma_lcdf.hpp
index 27e6800b444..7ad55a35302 100644
--- a/stan/math/prim/prob/inv_gamma_lcdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lcdf.hpp
@@ -66,9 +66,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lcdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lpdf.hpp b/stan/math/prim/prob/inv_gamma_lpdf.hpp
index 96af312d941..69f20d8aaf1 100644
--- a/stan/math/prim/prob/inv_gamma_lpdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lpdf.hpp
@@ -90,7 +90,7 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lpdf(const T_y& y,
 
   size_t N = max_size(y, alpha, beta);
   if (include_summand<propto, T_shape>::value) {
-    logp -= sum(lgamma(alpha_val)) * N / size(alpha);
+    logp -= sum(lgamma(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_shape, T_scale>::value) {
     const auto& log_beta
diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp
index 4016e423324..46854bd3be6 100644
--- a/stan/math/prim/prob/logistic_lpdf.hpp
+++ b/stan/math/prim/prob/logistic_lpdf.hpp
@@ -73,7 +73,7 @@ return_type_t<T_y, T_loc, T_scale> logistic_lpdf(const T_y& y, const T_loc& mu,
   T_partials_return logp = -sum(y_minus_mu_div_sigma)
                            - 2.0 * sum(log1p(exp(-y_minus_mu_div_sigma)));
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_scale>::value) {
diff --git a/stan/math/prim/prob/lognormal_lpdf.hpp b/stan/math/prim/prob/lognormal_lpdf.hpp
index edb61e3688d..facf5e3e1e5 100644
--- a/stan/math/prim/prob/lognormal_lpdf.hpp
+++ b/stan/math/prim/prob/lognormal_lpdf.hpp
@@ -81,10 +81,10 @@ return_type_t<T_y, T_loc, T_scale> lognormal_lpdf(const T_y& y, const T_loc& mu,
   T_partials_return logp
       = N * NEG_LOG_SQRT_TWO_PI - 0.5 * sum(square(logy_m_mu) * inv_sigma_sq);
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y>::value) {
-    logp -= sum(log_y) * N / size(y);
+    logp -= sum(log_y) * N / math::size(y);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp
index 31ca456c4ec..f0b92624450 100644
--- a/stan/math/prim/prob/normal_lpdf.hpp
+++ b/stan/math/prim/prob/normal_lpdf.hpp
@@ -92,7 +92,7 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
     logp += NEG_LOG_SQRT_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_scale, T_loc>::value) {
diff --git a/stan/math/prim/prob/normal_sufficient_lpdf.hpp b/stan/math/prim/prob/normal_sufficient_lpdf.hpp
index f4087fdd158..f741d72c11d 100644
--- a/stan/math/prim/prob/normal_sufficient_lpdf.hpp
+++ b/stan/math/prim/prob/normal_sufficient_lpdf.hpp
@@ -117,7 +117,7 @@ return_type_t<T_y, T_s, T_loc, T_scale> normal_sufficient_lpdf(
   size_t N = max_size(y_bar, s_squared, n_obs, mu, sigma);
   T_partials_return logp = -sum(cons_expr / (2 * sigma_squared));
   if (include_summand<propto>::value) {
-    logp += NEG_LOG_SQRT_TWO_PI * sum(n_obs_val) * N / size(n_obs);
+    logp += NEG_LOG_SQRT_TWO_PI * sum(n_obs_val) * N / math::size(n_obs);
   }
   if (include_summand<propto, T_scale>::value) {
     logp -= sum(n_obs_val * log(sigma_val)) * N / max_size(n_obs, sigma);
@@ -151,7 +151,7 @@ return_type_t<T_y, T_s, T_loc, T_scale> normal_sufficient_lpdf(
       } else {
         forward_as<internal::broadcast_array<T_partials_return>>(
             ops_partials.edge2_.partials_)
-            = -0.5 / sigma_squared * N / size(sigma);
+            = -0.5 / sigma_squared * N / math::size(sigma);
       }
     }
   }
diff --git a/stan/math/prim/prob/ordered_logistic_lpmf.hpp b/stan/math/prim/prob/ordered_logistic_lpmf.hpp
index b9a90dabb1f..35516fcffcf 100644
--- a/stan/math/prim/prob/ordered_logistic_lpmf.hpp
+++ b/stan/math/prim/prob/ordered_logistic_lpmf.hpp
@@ -86,7 +86,7 @@ return_type_t<T_loc, T_cut> ordered_logistic_lpmf(const T_y& y,
   T_cut_ref c_ref = c;
   vector_seq_view<T_cut_ref> c_vec(c_ref);
   int K = c_vec[0].size() + 1;
-  int N = size(lambda);
+  int N = math::size(lambda);
   int C_l = size_mvt(c);
 
   check_consistent_sizes(function, "Integers", y, "Locations", lambda);
diff --git a/stan/math/prim/prob/ordered_probit_lpmf.hpp b/stan/math/prim/prob/ordered_probit_lpmf.hpp
index b7a89401737..18bbfdd0174 100644
--- a/stan/math/prim/prob/ordered_probit_lpmf.hpp
+++ b/stan/math/prim/prob/ordered_probit_lpmf.hpp
@@ -54,7 +54,7 @@ return_type_t<T_loc, T_cut> ordered_probit_lpmf(const T_y& y,
   static const char* function = "ordered_probit";
 
   check_nonzero_size(function, "Cut-points", c);
-  int N = size(lambda);
+  int N = math::size(lambda);
   int C_l = size_mvt(c);
   vector_seq_view<T_cut> c_vec(c);
   int K = c_vec[0].size() + 1;
diff --git a/stan/math/prim/prob/pareto_lpdf.hpp b/stan/math/prim/prob/pareto_lpdf.hpp
index 9ef1b53ded0..c39b86f6a14 100644
--- a/stan/math/prim/prob/pareto_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_lpdf.hpp
@@ -70,7 +70,7 @@ return_type_t<T_y, T_scale, T_shape> pareto_lpdf(const T_y& y,
   size_t N = max_size(y, y_min, alpha);
   T_partials_return logp(0);
   if (include_summand<propto, T_shape>::value) {
-    logp = sum(log(alpha_val)) * N / size(alpha);
+    logp = sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp -= sum(alpha_val * log_y + log_y) * N / max_size(alpha, y);
diff --git a/stan/math/prim/prob/pareto_type_2_lpdf.hpp b/stan/math/prim/prob/pareto_type_2_lpdf.hpp
index efb213cc3c3..b6beef26c1d 100644
--- a/stan/math/prim/prob/pareto_type_2_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_type_2_lpdf.hpp
@@ -74,10 +74,10 @@ return_type_t<T_y, T_loc, T_scale, T_shape> pareto_type_2_lpdf(
   size_t N = max_size(y, mu, lambda, alpha);
   T_partials_return logp(0.0);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(lambda_val)) * N / size(lambda);
+    logp -= sum(log(lambda_val)) * N / math::size(lambda);
   }
   if (include_summand<propto, T_y, T_loc, T_scale, T_shape>::value) {
     logp -= sum((alpha_val + 1.0) * log1p_scaled_diff);
diff --git a/stan/math/prim/prob/poisson_log_lpmf.hpp b/stan/math/prim/prob/poisson_log_lpmf.hpp
index 0568c1dcc89..1f927685252 100644
--- a/stan/math/prim/prob/poisson_log_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_log_lpmf.hpp
@@ -78,10 +78,10 @@ return_type_t<T_log_rate> poisson_log_lpmf(const T_n& n,
 
   T_partials_return logp = sum(n_val * alpha_val);
   if (include_summand<propto, T_log_rate>::value) {
-    logp -= sum(exp_alpha) * N / size(alpha);
+    logp -= sum(exp_alpha) * N / math::size(alpha);
   }
   if (include_summand<propto>::value) {
-    logp -= sum(lgamma(n_val + 1.0)) * N / size(n);
+    logp -= sum(lgamma(n_val + 1.0)) * N / math::size(n);
   }
 
   if (!is_constant_all<T_log_rate>::value) {
diff --git a/stan/math/prim/prob/poisson_lpmf.hpp b/stan/math/prim/prob/poisson_lpmf.hpp
index f9bd2320a81..bf11e44ef0c 100644
--- a/stan/math/prim/prob/poisson_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_lpmf.hpp
@@ -72,10 +72,10 @@ return_type_t<T_rate> poisson_lpmf(const T_n& n, const T_rate& lambda) {
 
   T_partials_return logp = stan::math::sum(multiply_log(n_val, lambda_val));
   if (include_summand<propto, T_rate>::value) {
-    logp -= sum(lambda_val) * N / size(lambda);
+    logp -= sum(lambda_val) * N / math::size(lambda);
   }
   if (include_summand<propto>::value) {
-    logp -= sum(lgamma(n_val + 1.0)) * N / size(n);
+    logp -= sum(lgamma(n_val + 1.0)) * N / math::size(n);
   }
 
   if (!is_constant_all<T_rate>::value) {
diff --git a/stan/math/prim/prob/rayleigh_lpdf.hpp b/stan/math/prim/prob/rayleigh_lpdf.hpp
index 7582d3dd7ca..7d43a4a12a9 100644
--- a/stan/math/prim/prob/rayleigh_lpdf.hpp
+++ b/stan/math/prim/prob/rayleigh_lpdf.hpp
@@ -61,10 +61,10 @@ return_type_t<T_y, T_scale> rayleigh_lpdf(const T_y& y, const T_scale& sigma) {
   size_t N = max_size(y, sigma);
   T_partials_return logp = -0.5 * sum(square(y_over_sigma));
   if (include_summand<propto, T_scale>::value) {
-    logp -= 2.0 * sum(log(sigma_val)) * N / size(sigma);
+    logp -= 2.0 * sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y>::value) {
-    logp += sum(log(y_val)) * N / size(y);
+    logp += sum(log(y_val)) * N / math::size(y);
   }
 
   if (!is_constant_all<T_y, T_scale>::value) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
index f441852d649..7ec39b46e56 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
@@ -78,9 +78,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_cdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
index 83b31549806..52c122b6f4a 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lccdf(
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
index c25413b2831..a5f480a438c 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lcdf(
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
index 5187e4b03f5..8c03961ed3f 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
@@ -86,7 +86,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
                 T_partials_return, T_dof>
-      half_nu(size(nu));
+      half_nu(math::size(nu));
   for (size_t i = 0; i < stan::math::size(nu); i++) {
     if (include_summand<propto, T_dof, T_y, T_scale>::value) {
       half_nu[i] = 0.5 * nu_vec.val(i);
@@ -95,7 +95,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y>::value, T_partials_return,
                 T_y>
-      log_y(size(y));
+      log_y(math::size(y));
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (include_summand<propto, T_dof, T_y>::value) {
       log_y[i] = log(y_vec.val(i));
@@ -104,7 +104,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
                 T_partials_return, T_y>
-      inv_y(size(y));
+      inv_y(math::size(y));
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (include_summand<propto, T_dof, T_y, T_scale>::value) {
       inv_y[i] = 1.0 / y_vec.val(i);
@@ -113,7 +113,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_scale>::value,
                 T_partials_return, T_scale>
-      log_s(size(s));
+      log_s(math::size(s));
   for (size_t i = 0; i < stan::math::size(s); i++) {
     if (include_summand<propto, T_dof, T_scale>::value) {
       log_s[i] = log(s_vec.val(i));
@@ -121,11 +121,11 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
   }
 
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      log_half_nu(size(nu));
+      log_half_nu(math::size(nu));
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      lgamma_half_nu(size(nu));
+      lgamma_half_nu(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_half_nu_over_two(size(nu));
+      digamma_half_nu_over_two(math::size(nu));
   for (size_t i = 0; i < stan::math::size(nu); i++) {
     if (include_summand<propto, T_dof>::value) {
       lgamma_half_nu[i] = lgamma(half_nu[i]);
diff --git a/stan/math/prim/prob/skew_normal_lpdf.hpp b/stan/math/prim/prob/skew_normal_lpdf.hpp
index e7de3a48dd0..92d89406c48 100644
--- a/stan/math/prim/prob/skew_normal_lpdf.hpp
+++ b/stan/math/prim/prob/skew_normal_lpdf.hpp
@@ -85,7 +85,7 @@ return_type_t<T_y, T_loc, T_scale, T_shape> skew_normal_lpdf(
     logp -= HALF_LOG_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y, T_loc, T_scale>::value) {
     logp -= sum(square(y_minus_mu_over_sigma)) * 0.5 * N
diff --git a/stan/math/prim/prob/student_t_cdf.hpp b/stan/math/prim/prob/student_t_cdf.hpp
index 38232b2ea0a..87c50cccf44 100644
--- a/stan/math/prim/prob/student_t_cdf.hpp
+++ b/stan/math/prim/prob/student_t_cdf.hpp
@@ -65,11 +65,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_cdf(const T_y& y,
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lccdf.hpp b/stan/math/prim/prob/student_t_lccdf.hpp
index 5c9e39865d2..cfc6e5608ef 100644
--- a/stan/math/prim/prob/student_t_lccdf.hpp
+++ b/stan/math/prim/prob/student_t_lccdf.hpp
@@ -65,11 +65,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lccdf(
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lcdf.hpp b/stan/math/prim/prob/student_t_lcdf.hpp
index 401e938cc5b..84158993b47 100644
--- a/stan/math/prim/prob/student_t_lcdf.hpp
+++ b/stan/math/prim/prob/student_t_lcdf.hpp
@@ -67,11 +67,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lcdf(const T_y& y,
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lpdf.hpp b/stan/math/prim/prob/student_t_lpdf.hpp
index 3c39c228a42..7f203582914 100644
--- a/stan/math/prim/prob/student_t_lpdf.hpp
+++ b/stan/math/prim/prob/student_t_lpdf.hpp
@@ -120,10 +120,10 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lpdf(const T_y& y,
   if (include_summand<propto, T_dof>::value) {
     logp += (sum(lgamma(half_nu + 0.5)) - sum(lgamma(half_nu))
              - 0.5 * sum(log(nu_val)))
-            * N / size(nu);
+            * N / math::size(nu);
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/uniform_lcdf.hpp b/stan/math/prim/prob/uniform_lcdf.hpp
index 0962208e046..67e9876e3ee 100644
--- a/stan/math/prim/prob/uniform_lcdf.hpp
+++ b/stan/math/prim/prob/uniform_lcdf.hpp
@@ -73,7 +73,7 @@ return_type_t<T_y, T_low, T_high> uniform_lcdf(const T_y& y, const T_low& alpha,
   if (!is_constant_all<T_y>::value) {
     if (!is_vector<T_y>::value && is_vector<T_high>::value
         && !is_vector<T_low>::value) {
-      ops_partials.edge1_.partials_ = size(beta) * inv(y_minus_alpha);
+      ops_partials.edge1_.partials_ = math::size(beta) * inv(y_minus_alpha);
     } else {
       ops_partials.edge1_.partials_ = inv(y_minus_alpha);
     }
@@ -85,7 +85,7 @@ return_type_t<T_y, T_low, T_high> uniform_lcdf(const T_y& y, const T_low& alpha,
   if (!is_constant_all<T_high>::value) {
     if (is_vector<T_y>::value && !is_vector<T_low>::value
         && !is_vector<T_high>::value) {
-      ops_partials.edge3_.partials_ = inv(-b_minus_a) * size(y);
+      ops_partials.edge3_.partials_ = inv(-b_minus_a) * math::size(y);
     } else {
       ops_partials.edge3_.partials_ = inv(-b_minus_a);
     }
diff --git a/stan/math/prim/prob/uniform_lpdf.hpp b/stan/math/prim/prob/uniform_lpdf.hpp
index 73c875a17a3..837b912b4e7 100644
--- a/stan/math/prim/prob/uniform_lpdf.hpp
+++ b/stan/math/prim/prob/uniform_lpdf.hpp
@@ -102,7 +102,7 @@ return_type_t<T_y, T_low, T_high> uniform_lpdf(const T_y& y, const T_low& alpha,
     if (!is_constant_all<T_high>::value) {
       if (is_vector<T_y>::value && !is_vector<T_low>::value
           && !is_vector<T_high>::value) {
-        ops_partials.edge3_.partials_ = -inv_beta_minus_alpha * size(y);
+        ops_partials.edge3_.partials_ = -inv_beta_minus_alpha * math::size(y);
       } else {
         ops_partials.edge3_.partials_ = -inv_beta_minus_alpha;
       }
@@ -110,7 +110,7 @@ return_type_t<T_y, T_low, T_high> uniform_lpdf(const T_y& y, const T_low& alpha,
     if (!is_constant_all<T_low>::value) {
       if (is_vector<T_y>::value && !is_vector<T_low>::value
           && !is_vector<T_high>::value) {
-        ops_partials.edge2_.partials_ = inv_beta_minus_alpha * size(y);
+        ops_partials.edge2_.partials_ = inv_beta_minus_alpha * math::size(y);
       } else {
         ops_partials.edge2_.partials_ = std::move(inv_beta_minus_alpha);
       }
diff --git a/stan/math/prim/prob/von_mises_lpdf.hpp b/stan/math/prim/prob/von_mises_lpdf.hpp
index 8fe6650ce2c..910eef74347 100644
--- a/stan/math/prim/prob/von_mises_lpdf.hpp
+++ b/stan/math/prim/prob/von_mises_lpdf.hpp
@@ -70,7 +70,7 @@ return_type_t<T_y, T_loc, T_scale> von_mises_lpdf(T_y const& y, T_loc const& mu,
     logp -= LOG_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log_modified_bessel_first_kind(0, kappa_val)) * N / size(kappa);
+    logp -= sum(log_modified_bessel_first_kind(0, kappa_val)) * N / math::size(kappa);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/weibull_lpdf.hpp b/stan/math/prim/prob/weibull_lpdf.hpp
index e355beeac2f..72e4bbc0294 100644
--- a/stan/math/prim/prob/weibull_lpdf.hpp
+++ b/stan/math/prim/prob/weibull_lpdf.hpp
@@ -96,7 +96,7 @@ return_type_t<T_y, T_shape, T_scale> weibull_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, sigma);
   T_partials_return logp = -sum(y_div_sigma_pow_alpha);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp += sum((alpha_val - 1.0) * log_y) * N / max_size(alpha, y);
diff --git a/stan/math/rev/functor/coupled_ode_system.hpp b/stan/math/rev/functor/coupled_ode_system.hpp
index 83783d4d03f..3ea68004a65 100644
--- a/stan/math/rev/functor/coupled_ode_system.hpp
+++ b/stan/math/rev/functor/coupled_ode_system.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/functor/coupled_ode_system.hpp>
 #include <stan/math/rev/functor/cvodes_utils.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
@@ -125,7 +126,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
       y_vars.coeffRef(n) = z[n];
 
     Eigen::Matrix<var, Eigen::Dynamic, 1> f_y_t_vars
-        = apply([&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
+        = math::apply(([&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
                 local_args_tuple_);
 
     check_size_match("coupled_ode_system", "dy_dt", f_y_t_vars.size(), "states",
@@ -140,7 +141,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
       // memset was faster than Eigen setZero
       memset(args_adjoints_.data(), 0, sizeof(double) * num_args_vars);
 
-      apply(
+      math::apply((
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },
@@ -148,7 +149,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
 
       // The vars here do not live on the nested stack so must be zero'd
       // separately
-      apply([&](auto&&... args) { zero_adjoints(args...); }, local_args_tuple_);
+      math::apply(([&](auto&&... args) { zero_adjoints(args...); }, local_args_tuple_);
 
       // No need to zero adjoints after last sweep
       if (i + 1 < N_) {
diff --git a/stan/math/rev/functor/cvodes_integrator.hpp b/stan/math/rev/functor/cvodes_integrator.hpp
index 4e40354540b..a174205d238 100644
--- a/stan/math/rev/functor/cvodes_integrator.hpp
+++ b/stan/math/rev/functor/cvodes_integrator.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/rev/functor/cvodes_utils.hpp>
 #include <stan/math/rev/functor/coupled_ode_system.hpp>
 #include <stan/math/rev/functor/ode_store_sensitivities.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <cvodes/cvodes.h>
@@ -103,7 +104,7 @@ class cvodes_integrator {
     const Eigen::VectorXd y_vec = Eigen::Map<const Eigen::VectorXd>(y, N_);
 
     Eigen::VectorXd dy_dt_vec
-        = apply([&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
+        = math::apply(([&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
                 value_of_args_tuple_);
 
     check_size_match("cvodes_integrator", "dy_dt", dy_dt_vec.size(), "states",
@@ -121,7 +122,7 @@ class cvodes_integrator {
     Eigen::MatrixXd Jfy;
 
     auto f_wrapped = [&](const Eigen::Matrix<var, Eigen::Dynamic, 1>& y) {
-      return apply([&](auto&&... args) { return f_(t, y, msgs_, args...); },
+      return math::apply(([&](auto&&... args) { return f_(t, y, msgs_, args...); },
                    value_of_args_tuple_);
     };
 
@@ -211,7 +212,7 @@ class cvodes_integrator {
 
     // Code from: https://stackoverflow.com/a/17340003 . Should probably do
     // something better
-    apply(
+    math::apply((
         [&](auto&&... args) {
           std::vector<int> unused_temp{
               0, (check_finite(function_name, "ode parameters and data", args),
@@ -326,7 +327,7 @@ class cvodes_integrator {
           }
         }
 
-        y.emplace_back(apply(
+        y.emplace_back(math::apply((
             [&](auto&&... args) {
               return ode_store_sensitivities(f_, coupled_state_, y0_, t0_,
                                              ts_[n], msgs_, args...);
diff --git a/stan/math/rev/functor/ode_adams.hpp b/stan/math/rev/functor/ode_adams.hpp
index ee0bdafbbd5..204adab6bce 100644
--- a/stan/math/rev/functor/ode_adams.hpp
+++ b/stan/math/rev/functor/ode_adams.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/functor/cvodes_integrator.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <ostream>
 #include <vector>
 
@@ -54,7 +55,7 @@ ode_adams_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                    long int max_num_steps,  // NOLINT(runtime/int)
                    std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return apply(
+  return math::apply((
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_ADAMS, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/ode_bdf.hpp b/stan/math/rev/functor/ode_bdf.hpp
index a07af2e3339..75c8b6bc197 100644
--- a/stan/math/rev/functor/ode_bdf.hpp
+++ b/stan/math/rev/functor/ode_bdf.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/functor/cvodes_integrator.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <ostream>
 #include <vector>
@@ -55,7 +56,7 @@ ode_bdf_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                  long int max_num_steps,  // NOLINT(runtime/int)
                  std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return apply(
+  return math::apply((
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_BDF, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/reduce_sum.hpp b/stan/math/rev/functor/reduce_sum.hpp
index dfe51dee2a7..bfecd4017e2 100644
--- a/stan/math/rev/functor/reduce_sum.hpp
+++ b/stan/math/rev/functor/reduce_sum.hpp
@@ -113,7 +113,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
         // scope. In this case no need for zeroing adjoints, since the
         // fresh copy has all adjoints set to zero.
         local_args_tuple_scope_.stack_.execute([&]() {
-          apply(
+          math::apply((
               [&](auto&&... args) {
                 local_args_tuple_scope_.args_tuple_holder_ = std::make_unique<
                     typename scoped_args_tuple::args_tuple_t>(
@@ -140,7 +140,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
       }
 
       // Perform calculation
-      var sub_sum_v = apply(
+      var sub_sum_v = math::apply((
           [&](auto&&... args) {
             return ReduceFunction()(local_sub_slice, r.begin(), r.end() - 1,
                                     &msgs_, args...);
@@ -158,7 +158,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
                           std::move(local_sub_slice));
 
       // Accumulate adjoints of shared_arguments
-      apply(
+      math::apply((
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },

From 9b8d4096921d55d836743f25aca8251755036254 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 15:55:13 +0200
Subject: [PATCH 3/9] Backport batch 3

---
 stan/math/rev/functor/adj_jac_apply.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stan/math/rev/functor/adj_jac_apply.hpp b/stan/math/rev/functor/adj_jac_apply.hpp
index 66826b7a7b9..5d81bc6ad61 100644
--- a/stan/math/rev/functor/adj_jac_apply.hpp
+++ b/stan/math/rev/functor/adj_jac_apply.hpp
@@ -509,7 +509,7 @@ struct adj_jac_vari : public vari {
     internal::build_y_adj(y_vi_, M_, y_adj);
     auto y_adj_jacs = f_.multiply_adjoint_jacobian(is_var_, y_adj);
 
-    apply([this](auto&&... args) { this->accumulate_adjoints(args...); },
+    math::apply([this](auto&&... args) { this->accumulate_adjoints(args...); },
           y_adj_jacs);
   }
 };

From 6cae73e8f97993764742afe933eaa4205d8ed858 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 16:04:14 +0200
Subject: [PATCH 4/9] Revert additional backports

---
 stan/math/opencl/matrix_cl.hpp                | 81 +++++++------------
 .../opencl/prim/bernoulli_logit_glm_lpmf.hpp  | 14 ++--
 .../prim/categorical_logit_glm_lpmf.hpp       |  6 +-
 .../prim/neg_binomial_2_log_glm_lpmf.hpp      | 18 ++---
 stan/math/opencl/prim/normal_id_glm_lpdf.hpp  | 12 +--
 .../opencl/prim/ordered_logistic_glm_lpmf.hpp |  4 +-
 .../math/opencl/prim/poisson_log_glm_lpmf.hpp | 17 ++--
 stan/math/opencl/prim/poisson_log_lpmf.hpp    |  6 +-
 stan/math/opencl/prim/poisson_lpmf.hpp        |  4 +-
 9 files changed, 71 insertions(+), 91 deletions(-)

diff --git a/stan/math/opencl/matrix_cl.hpp b/stan/math/opencl/matrix_cl.hpp
index 6c5470e055a..561694cda6f 100644
--- a/stan/math/opencl/matrix_cl.hpp
+++ b/stan/math/opencl/matrix_cl.hpp
@@ -11,7 +11,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/vec_concat.hpp>
-#include <CL/opencl.hpp>
+#include <CL/cl2.hpp>
 #include <algorithm>
 #include <iostream>
 #include <string>
@@ -31,11 +31,7 @@ namespace math {
  *  @{
  */
 
-// forward declare
-template <typename T>
-class arena_matrix_cl;
-
-template <typename>
+template <typename, typename = void>
 class matrix_cl;
 
 /**
@@ -43,7 +39,7 @@ class matrix_cl;
  * @tparam T an arithmetic type for the type stored in the OpenCL buffer.
  */
 template <typename T>
-class matrix_cl : public matrix_cl_base {
+class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
  private:
   cl::Buffer buffer_cl_;  // Holds the allocated memory on the device
   int rows_{0};           // Number of rows.
@@ -59,6 +55,8 @@ class matrix_cl : public matrix_cl_base {
   // Forward declare the methods that work in place on the matrix
   template <matrix_cl_view matrix_view = matrix_cl_view::Entire>
   inline void zeros_strict_tri();
+  template <TriangularMapCL triangular_map = TriangularMapCL::LowerToUpper>
+  inline void triangular_transpose();
 
   int rows() const { return rows_; }
 
@@ -199,9 +197,7 @@ class matrix_cl : public matrix_cl_base {
     if (A.size() == 0) {
       return;
     }
-    buffer_cl_ = cl::Buffer(opencl_context.context(), CL_MEM_READ_WRITE,
-                            sizeof(T) * this->size());
-    initialize_buffer_cl(A);
+    initialize_buffer(A);
   }
 
   /**
@@ -216,13 +212,6 @@ class matrix_cl : public matrix_cl_base {
         write_events_(std::move(A.write_events_)),
         read_events_(std::move(A.read_events_)) {}
 
-  /**
-   * Constructor from `arena_matrix_cl`.
-   * @param A matrix_cl to move
-   */
-  // defined in rev/arena_matrix_cl.hpp
-  matrix_cl(const arena_matrix_cl<T>& A);  // NOLINT(runtime/explicit)
-
   /**
    * Constructor for the matrix_cl that creates a copy of a std::vector of Eigen
    * matrices on the OpenCL device. Each matrix is flattened into one column
@@ -283,7 +272,7 @@ class matrix_cl : public matrix_cl_base {
   matrix_cl(const int rows, const int cols,
             matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(rows), cols_(cols), view_(partial_view) {
-    if (this->size() == 0) {
+    if this->size() == 0) {
       return;
     }
     cl::Context& ctx = opencl_context.context();
@@ -321,7 +310,7 @@ class matrix_cl : public matrix_cl_base {
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(A.rows()), cols_(A.cols()), view_(partial_view) {
     using Mat_type = std::decay_t<ref_type_for_opencl_t<Mat>>;
-    if (this->size() == 0) {
+    if this->size() == 0) {
       return;
     }
     initialize_buffer_no_heap_if<
@@ -426,10 +415,8 @@ class matrix_cl : public matrix_cl_base {
    * @tparam Expr type of the expression
    * @param expression expression
    */
-  // defined in kernel_generator/matrix_cl_conversion.hpp
   template <typename Expr,
-            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr,
-            require_not_matrix_cl_t<Expr>* = nullptr>
+            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr>
   matrix_cl(const Expr& expression);  // NOLINT(runtime/explicit)
 
   /**
@@ -451,19 +438,13 @@ class matrix_cl : public matrix_cl_base {
    */
   matrix_cl<T>& operator=(const matrix_cl<T>& a) {
     this->view_ = a.view();
+    this->rows_ = a.rows();
+    this->cols_ = a.cols();
     if (a.size() == 0) {
-      this->rows_ = a.rows();
-      this->cols_ = a.cols();
       return *this;
     }
     this->wait_for_read_write_events();
-    if (this->size() != a.size()) {
-      buffer_cl_ = cl::Buffer(opencl_context.context(), CL_MEM_READ_WRITE,
-                              sizeof(T) * a.size());
-    }
-    this->rows_ = a.rows();
-    this->cols_ = a.cols();
-    initialize_buffer_cl(a);
+    initialize_buffer(a);
     return *this;
   }
 
@@ -473,20 +454,10 @@ class matrix_cl : public matrix_cl_base {
    * @tparam Expr type of the expression
    * @param expression expression
    */
-  // defined in kernel_generator/matrix_cl_conversion.hpp
   template <typename Expr,
-            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr,
-            require_not_matrix_cl_t<Expr>* = nullptr>
+            require_all_kernel_expressions_and_none_scalar_t<Expr>* = nullptr>
   matrix_cl<T>& operator=(const Expr& expression);
 
-  /**
-   * Assignment of `arena_matrix_cl<T>`.
-   * @tparam Expr type of the expression
-   * @param expression expression
-   */
-  // defined in rev/arena_matrix_cl.hpp
-  matrix_cl<T>& operator=(const arena_matrix_cl<T>& other);
-
   /**
    * Evaluates `this`. This is a no-op.
    * @return `*this`
@@ -518,7 +489,7 @@ class matrix_cl : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(const T* A) {
     cl::Event transfer_event;
-    if (this->size() == 0) {
+    if this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -538,7 +509,7 @@ class matrix_cl : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(T* A) {
     cl::Event transfer_event;
-    if (this->size() == 0) {
+    if this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -577,7 +548,7 @@ class matrix_cl : public matrix_cl_base {
    */
   template <bool No_heap, typename U, std::enable_if_t<No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
-    if (this->size() == 0) {
+    if this->size() == 0) {
       return;
     }
     initialize_buffer(obj.data());
@@ -587,7 +558,7 @@ class matrix_cl : public matrix_cl_base {
   template <bool No_heap, typename U, std::enable_if_t<!No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
     using U_val = std::decay_t<ref_type_for_opencl_t<U>>;
-    if (this->size() == 0) {
+    if this->size() == 0) {
       return;
     }
     auto* obj_heap = new U_val(std::move(obj));
@@ -611,12 +582,15 @@ class matrix_cl : public matrix_cl_base {
    * size of given matrix.
    * @param A matrix_cl
    */
-  void initialize_buffer_cl(const matrix_cl<T>& A) {
+  void initialize_buffer(const matrix_cl<T>& A) {
+    cl::Context& ctx = opencl_context.context();
+    cl::CommandQueue queue = opencl_context.queue();
     try {
+      buffer_cl_ = cl::Buffer(ctx, CL_MEM_READ_WRITE, sizeof(T) * this->size());
       cl::Event cstr_event;
-      opencl_context.queue().enqueueCopyBuffer(A.buffer(), this->buffer(), 0, 0,
-                                               A.size() * sizeof(T),
-                                               &A.write_events(), &cstr_event);
+      queue.enqueueCopyBuffer(A.buffer(), this->buffer(), 0, 0,
+                              A.size() * sizeof(T), &A.write_events(),
+                              &cstr_event);
       this->add_write_event(cstr_event);
       A.add_read_event(cstr_event);
     } catch (const cl::Error& e) {
@@ -647,6 +621,13 @@ class matrix_cl : public matrix_cl_base {
     delete static_cast<U*>(container);
   }
 };
+
+template <typename T>
+using matrix_cl_prim = matrix_cl<T, require_arithmetic_t<T>>;
+
+template <typename T>
+using matrix_cl_fp = matrix_cl<T, require_floating_point_t<T>>;
+
 /** @}*/
 
 }  // namespace math
diff --git a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
index 9de5422b293..47228286ac2 100644
--- a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
@@ -66,8 +66,7 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
-                     math::size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
                    math::size(beta));
@@ -121,7 +120,7 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
       logp_expr, calc_if<need_theta_derivative>(theta_derivative_expr),
       calc_if<need_theta_derivative_sum>(colwise_sum(theta_derivative_expr)));
 
-  T_partials_return logp = sum(from_matrix_cl(logp_cl));
+  T_partials_return logp = sum(from_matrix_cl<Eigen::Dynamic, 1>(logp_cl));
   if (!std::isfinite(logp)) {
     results(check_cl(function, "Vector of dependent variables", y_val,
                      "in the interval [0, 1]"),
@@ -153,13 +152,14 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
     // transposition of a vector can be done without copying
     const matrix_cl<double> theta_derivative_transpose_cl(
         theta_derivative_cl.buffer(), 1, theta_derivative_cl.rows());
+    matrix_cl<double>& edge3_partials
+        = forward_as<matrix_cl<double>&>(ops_partials.edge3_.partials_);
     matrix_cl<double> edge3_partials_transpose_cl
         = theta_derivative_transpose_cl * x_val;
-    ops_partials.edge3_.partials_
-        = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
-                            edge3_partials_transpose_cl.cols(), 1);
+    edge3_partials = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
+                                       edge3_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      ops_partials.edge3_.partials_.add_write_event(
+      edge3_partials.add_write_event(
           edge3_partials_transpose_cl.write_events().back());
     }
   }
diff --git a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
index baf4974f60d..b6be840c7ec 100644
--- a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
@@ -3,7 +3,6 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/opencl/prim/size.hpp>
-#include <stan/math/opencl/rev/arena_matrix_cl.hpp>
 #include <stan/math/opencl/rev/operands_and_partials.hpp>
 #include <stan/math/opencl/matrix_cl.hpp>
 #include <stan/math/opencl/copy.hpp>
@@ -151,9 +150,8 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
     try {
       opencl_kernels::categorical_logit_glm_beta_derivative(
           cl::NDRange(local_size * N_attributes), cl::NDRange(local_size),
-          forward_as<arena_matrix_cl<double>>(ops_partials.edge3_.partials_),
-          temp, y_val_cl, x_val, N_instances, N_attributes, N_classes,
-          is_y_vector);
+          forward_as<matrix_cl<double>>(ops_partials.edge3_.partials_), temp,
+          y_val_cl, x_val, N_instances, N_attributes, N_classes, is_y_vector);
     } catch (const cl::Error& e) {
       check_opencl_error(function, e);
     }
diff --git a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
index 1bfec8258b6..35a3e734ab0 100644
--- a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
@@ -82,8 +82,7 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
-                     math::size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x", M, "size of ", "beta",
                    math::size(beta));
@@ -150,7 +149,7 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     check_opencl_error(function, e);
   }
 
-  T_partials_return logp = sum(from_matrix_cl(logp_cl));
+  T_partials_return logp = sum(from_matrix_cl<Dynamic, 1>(logp_cl));
   if (!std::isfinite(logp)) {
     results(
         check_cl(function, "Vector of dependent variables", y_val,
@@ -189,13 +188,14 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     // transposition of a vector can be done without copying
     const matrix_cl<double> theta_derivative_transpose_cl(
         theta_derivative_cl.buffer(), 1, theta_derivative_cl.rows());
+    matrix_cl<double>& edge3_partials
+        = forward_as<matrix_cl<double>&>(ops_partials.edge3_.partials_);
     matrix_cl<double> edge3_partials_transpose_cl
         = theta_derivative_transpose_cl * x_val;
-    ops_partials.edge3_.partials_
-        = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
-                            edge3_partials_transpose_cl.cols(), 1);
+    edge3_partials = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
+                                       edge3_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      ops_partials.edge3_.partials_.add_write_event(
+      edge3_partials.add_write_event(
           edge3_partials_transpose_cl.write_events().back());
     }
   }
@@ -205,7 +205,7 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     } else {
       forward_as<internal::broadcast_array<double>>(
           ops_partials.edge2_.partials_)[0]
-          = sum(from_matrix_cl(theta_derivative_sum_cl));
+          = sum(from_matrix_cl<Dynamic, 1>(theta_derivative_sum_cl));
     }
   }
   if (!is_constant_all<T_phi_cl>::value) {
@@ -214,7 +214,7 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
     } else {
       forward_as<internal::broadcast_array<double>>(
           ops_partials.edge4_.partials_)[0]
-          = sum(from_matrix_cl(phi_derivative_cl));
+          = sum(from_matrix_cl<Dynamic, 1>(phi_derivative_cl));
     }
   }
   return ops_partials.build(logp);
diff --git a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
index b999c25e535..2d98f29a163 100644
--- a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
+++ b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
@@ -74,8 +74,7 @@ normal_id_glm_lpdf(const T_y_cl& y, const T_x_cl& x, const T_alpha_cl& alpha,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
-                     math::size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
                    math::size(beta));
@@ -172,13 +171,14 @@ normal_id_glm_lpdf(const T_y_cl& y, const T_x_cl& x, const T_alpha_cl& alpha,
     // transposition of a vector can be done without copying
     const matrix_cl<double> mu_derivative_transpose_cl(
         mu_derivative_cl.buffer(), 1, mu_derivative_cl.rows());
+    matrix_cl<double>& edge4_partials
+        = forward_as<matrix_cl<double>&>(ops_partials.edge4_.partials_);
     matrix_cl<double> edge4_partials_transpose_cl
         = mu_derivative_transpose_cl * x_val;
-    ops_partials.edge4_.partials_
-        = matrix_cl<double>(edge4_partials_transpose_cl.buffer(),
-                            edge4_partials_transpose_cl.cols(), 1);
+    edge4_partials = matrix_cl<double>(edge4_partials_transpose_cl.buffer(),
+                                       edge4_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      ops_partials.edge4_.partials_.add_write_event(
+      edge4_partials.add_write_event(
           edge4_partials_transpose_cl.write_events().back());
     }
   }
diff --git a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
index 5f5912b2328..5fb69cd6047 100644
--- a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
@@ -140,8 +140,8 @@ return_type_t<T_x, T_beta, T_cuts> ordered_logistic_glm_lpmf(
         edge2_partials_transpose.buffer(), edge2_partials_transpose.cols(),
         edge2_partials_transpose.rows());
     if (beta.rows() != 0) {
-      ops_partials.edge2_.partials_.add_write_event(
-          edge2_partials_transpose.write_events().back());
+      forward_as<matrix_cl<double>>(ops_partials.edge2_.partials_)
+          .add_write_event(edge2_partials_transpose.write_events().back());
     }
   }
   if (!is_constant_all<T_cuts>::value) {
diff --git a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
index c6b730e7644..16e4774eb32 100644
--- a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
@@ -66,8 +66,7 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y",
-                     math::size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
                    math::size(beta));
@@ -109,8 +108,9 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
   results(theta_derivative_cl, theta_derivative_sum_cl, logp_cl) = expressions(
       theta_derivative_expr, colwise_sum(theta_derivative_expr), logp_expr);
 
-  double theta_derivative_sum = sum(from_matrix_cl(theta_derivative_sum_cl));
-  logp += sum(from_matrix_cl(logp_cl));
+  double theta_derivative_sum
+      = sum(from_matrix_cl<Dynamic, 1>(theta_derivative_sum_cl));
+  logp += sum(from_matrix_cl<Dynamic, 1>(logp_cl));
   if (!std::isfinite(theta_derivative_sum)) {
     results(check_cl(function, "Vector of dependent variables", y_val,
                      "nonnegative"),
@@ -142,13 +142,14 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
     // transposition of a vector can be done without copying
     const matrix_cl<double> theta_derivative_transpose_cl(
         theta_derivative_cl.buffer(), 1, theta_derivative_cl.rows());
+    matrix_cl<double>& edge3_partials
+        = forward_as<matrix_cl<double>&>(ops_partials.edge3_.partials_);
     matrix_cl<double> edge3_partials_transpose_cl
         = theta_derivative_transpose_cl * x_val;
-    ops_partials.edge3_.partials_
-        = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
-                            edge3_partials_transpose_cl.cols(), 1);
+    edge3_partials = matrix_cl<double>(edge3_partials_transpose_cl.buffer(),
+                                       edge3_partials_transpose_cl.cols(), 1);
     if (beta_val.rows() != 0) {
-      ops_partials.edge3_.partials_.add_write_event(
+      edge3_partials.add_write_event(
           edge3_partials_transpose_cl.write_events().back());
     }
   }
diff --git a/stan/math/opencl/prim/poisson_log_lpmf.hpp b/stan/math/opencl/prim/poisson_log_lpmf.hpp
index 31ed3dfaecd..c2fe7d1a6a6 100644
--- a/stan/math/opencl/prim/poisson_log_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_lpmf.hpp
@@ -60,8 +60,8 @@ return_type_t<T_log_rate_cl> poisson_log_lpmf(const T_n_cl& n,
       = check_cl(function, "Log rate parameter", alpha_val, "not nan");
   auto alpha_not_nan = !isnan(alpha_val);
 
-  auto return_log_zero
-      = colwise_max(cast<char>(isinf(alpha_val) && (alpha_val > 0 || n != 0)));
+  auto return_log_zero = colwise_max(
+      constant(0, N, 1) + (isinf(alpha_val) && (alpha_val > 0 || n != 0)));
   auto exp_alpha = exp(alpha_val);
 
   auto logp1 = elt_multiply(n, alpha_val);
@@ -72,7 +72,7 @@ return_type_t<T_log_rate_cl> poisson_log_lpmf(const T_n_cl& n,
 
   auto deriv = n - exp_alpha;
 
-  matrix_cl<char> return_log_zero_cl;
+  matrix_cl<int> return_log_zero_cl;
   matrix_cl<double> logp_cl;
   matrix_cl<double> deriv_cl;
 
diff --git a/stan/math/opencl/prim/poisson_lpmf.hpp b/stan/math/opencl/prim/poisson_lpmf.hpp
index 2878e9d2849..e4b047f0f34 100644
--- a/stan/math/opencl/prim/poisson_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_lpmf.hpp
@@ -60,7 +60,7 @@ return_type_t<T_rate_cl> poisson_lpmf(const T_n_cl& n,
   auto lambda_nonnegative = 0.0 <= lambda_val;
 
   auto return_log_zero = colwise_max(
-      cast<char>(isinf(lambda_val) || (lambda_val == 0 && n != 0)));
+      constant(0, N, 1) + (isinf(lambda_val) || (lambda_val == 0 && n != 0)));
 
   auto logp1 = multiply_log(n, lambda_val);
   auto logp2 = static_select<include_summand<propto, T_rate_cl>::value>(
@@ -70,7 +70,7 @@ return_type_t<T_rate_cl> poisson_lpmf(const T_n_cl& n,
 
   auto deriv = elt_divide(n, lambda_val) - 1.0;
 
-  matrix_cl<char> return_log_zero_cl;
+  matrix_cl<int> return_log_zero_cl;
   matrix_cl<double> logp_cl;
   matrix_cl<double> deriv_cl;
 

From a9b1dbbf02f0f3b80ffbd42af057813c8084ebc1 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 16:12:49 +0200
Subject: [PATCH 5/9] Revert extra backports

---
 stan/math/opencl/zeros_strict_tri.hpp         |   4 +-
 .../prim/fun/offset_multiplier_constrain.hpp  | 264 ++----------------
 stan/math/prim/prob/bernoulli_lpmf.hpp        |   7 +-
 stan/math/prim/prob/beta_proportion_lpdf.hpp  |  20 +-
 stan/math/prim/prob/cauchy_lpdf.hpp           |  20 +-
 stan/math/prim/prob/chi_square_lpdf.hpp       |  13 +-
 .../prim/prob/double_exponential_lpdf.hpp     |  16 +-
 stan/math/prim/prob/exp_mod_normal_lpdf.hpp   |  20 +-
 stan/math/prim/prob/exponential_lccdf.hpp     |  24 +-
 stan/math/prim/prob/exponential_lpdf.hpp      |  13 +-
 stan/math/prim/prob/frechet_lpdf.hpp          |  15 +-
 stan/math/prim/prob/gamma_lpdf.hpp            |  15 +-
 stan/math/prim/prob/gumbel_lpdf.hpp           |  15 +-
 13 files changed, 154 insertions(+), 292 deletions(-)

diff --git a/stan/math/opencl/zeros_strict_tri.hpp b/stan/math/opencl/zeros_strict_tri.hpp
index 3ff0b92842a..abf7349339a 100644
--- a/stan/math/opencl/zeros_strict_tri.hpp
+++ b/stan/math/opencl/zeros_strict_tri.hpp
@@ -9,7 +9,7 @@
 #include <stan/math/opencl/kernels/fill_strict_tri.hpp>
 #include <stan/math/prim/meta.hpp>
 
-#include <CL/opencl.hpp>
+#include <CL/cl2.hpp>
 
 namespace stan {
 namespace math {
@@ -29,7 +29,7 @@ namespace math {
  */
 template <typename T>
 template <matrix_cl_view matrix_view>
-inline void matrix_cl<T>::zeros_strict_tri() try {
+inline void matrix_cl<T, require_arithmetic_t<T>>::zeros_strict_tri() try {
   if (matrix_view == matrix_cl_view::Entire) {
     invalid_argument(
         "zeros_strict_tri", "matrix_view",
diff --git a/stan/math/prim/fun/offset_multiplier_constrain.hpp b/stan/math/prim/fun/offset_multiplier_constrain.hpp
index 85a58ce4771..c724610deda 100644
--- a/stan/math/prim/fun/offset_multiplier_constrain.hpp
+++ b/stan/math/prim/fun/offset_multiplier_constrain.hpp
@@ -7,9 +7,6 @@
 #include <stan/math/prim/fun/identity_constrain.hpp>
 #include <stan/math/prim/fun/multiply_log.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <stan/math/prim/fun/sum.hpp>
-#include <stan/math/prim/fun/to_ref.hpp>
-#include <stan/math/prim/functor/apply.hpp>
 #include <cmath>
 
 namespace stan {
@@ -38,27 +35,19 @@ namespace math {
  * @throw std::domain_error if sigma <= 0
  * @throw std::domain_error if mu is not finite
  */
-template <typename T, typename M, typename S,
-          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
-              T, M, S>* = nullptr>
-inline auto offset_multiplier_constrain(const T& x, const M& mu,
-                                        const S& sigma) {
-  const auto& mu_ref = to_ref(mu);
-  const auto& sigma_ref = to_ref(sigma);
-  if (is_matrix<T>::value && is_matrix<M>::value) {
-    check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
-  }
-  if (is_matrix<T>::value && is_matrix<S>::value) {
-    check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
-  } else if (is_matrix<M>::value && is_matrix<S>::value) {
-    check_matching_dims("offset_multiplier_constrain", "mu", mu, "sigma",
-                        sigma);
+template <typename T, typename M, typename S>
+inline return_type_t<T, M, S> offset_multiplier_constrain(const T& x,
+                                                          const M& mu,
+                                                          const S& sigma) {
+  check_finite("offset_multiplier_constrain", "offset", mu);
+  if (sigma == 1) {
+    if (mu == 0) {
+      return identity_constrain(x);
+    }
+    return mu + x;
   }
-
-  check_finite("offset_multiplier_constrain", "offset", value_of_rec(mu_ref));
-  check_positive_finite("offset_multiplier_constrain", "multiplier",
-                        value_of_rec(sigma_ref));
-  return fma(sigma_ref, x, mu_ref);
+  check_positive_finite("offset_multiplier_constrain", "multiplier", sigma);
+  return fma(sigma, x, mu);
 }
 
 /**
@@ -87,223 +76,22 @@ inline auto offset_multiplier_constrain(const T& x, const M& mu,
  * @throw std::domain_error if sigma <= 0
  * @throw std::domain_error if mu is not finite
  */
-template <typename T, typename M, typename S,
-          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
-              T, M, S>* = nullptr>
-inline auto offset_multiplier_constrain(const T& x, const M& mu, const S& sigma,
-                                        return_type_t<T, M, S>& lp) {
-  const auto& mu_ref = to_ref(mu);
-  const auto& sigma_ref = to_ref(sigma);
-  if (is_matrix<T>::value && is_matrix<M>::value) {
-    check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
-  }
-  if (is_matrix<T>::value && is_matrix<S>::value) {
-    check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
-  } else if (is_matrix<M>::value && is_matrix<S>::value) {
-    check_matching_dims("offset_multiplier_constrain", "mu", mu, "sigma",
-                        sigma);
-  }
-
-  check_finite("offset_multiplier_constrain", "offset", value_of_rec(mu_ref));
-  check_positive_finite("offset_multiplier_constrain", "multiplier",
-                        value_of_rec(sigma_ref));
-  if (math::size(sigma_ref) == 1) {
-    lp += sum(multiply_log(math::size(x), sigma_ref));
-  } else {
-    lp += sum(log(sigma_ref));
-  }
-  return fma(sigma_ref, x, mu_ref);
-}
-
-/**
- * Overload for array of x and non-array mu and sigma
- */
-template <typename T, typename M, typename S,
-          require_all_not_std_vector_t<M, S>* = nullptr>
-inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
-                                        const S& sigma) {
-  std::vector<
-      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu, sigma))>>
-      ret;
-  ret.reserve(x.size());
-  const auto& mu_ref = to_ref(mu);
-  const auto& sigma_ref = to_ref(sigma);
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma_ref));
-  }
-  return ret;
-}
-
-/**
- * Overload for array of x and non-array mu and sigma with lp
- */
-template <typename T, typename M, typename S,
-          require_all_not_std_vector_t<M, S>* = nullptr>
-inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
-                                        const S& sigma,
-                                        return_type_t<T, M, S>& lp) {
-  std::vector<
-      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu, sigma, lp))>>
-      ret;
-  ret.reserve(x.size());
-  const auto& mu_ref = to_ref(mu);
-  const auto& sigma_ref = to_ref(sigma);
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma_ref, lp));
-  }
-  return ret;
-}
-
-/**
- * Overload for array of x and sigma and non-array mu
- */
-template <typename T, typename M, typename S,
-          require_not_std_vector_t<M>* = nullptr>
-inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
-                                        const std::vector<S>& sigma) {
-  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
-  std::vector<
-      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu, sigma[0]))>>
-      ret;
-  ret.reserve(x.size());
-  const auto& mu_ref = to_ref(mu);
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma[i]));
-  }
-  return ret;
-}
-
-/**
- * Overload for array of x and sigma and non-array mu with lp
- */
-template <typename T, typename M, typename S,
-          require_not_std_vector_t<M>* = nullptr>
-inline auto offset_multiplier_constrain(const std::vector<T>& x, const M& mu,
-                                        const std::vector<S>& sigma,
-                                        return_type_t<T, M, S>& lp) {
-  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
-  std::vector<plain_type_t<decltype(
-      offset_multiplier_constrain(x[0], mu, sigma[0], lp))>>
-      ret;
-  ret.reserve(x.size());
-  const auto& mu_ref = to_ref(mu);
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu_ref, sigma[i], lp));
-  }
-  return ret;
-}
-
-/**
- * Overload for array of x and mu and non-array sigma
- */
-template <typename T, typename M, typename S,
-          require_not_std_vector_t<S>* = nullptr>
-inline auto offset_multiplier_constrain(const std::vector<T>& x,
-                                        const std::vector<M>& mu,
-                                        const S& sigma) {
-  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
-  std::vector<
-      plain_type_t<decltype(offset_multiplier_constrain(x[0], mu[0], sigma))>>
-      ret;
-  ret.reserve(x.size());
-  const auto& sigma_ref = to_ref(sigma);
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma_ref));
-  }
-  return ret;
-}
-
-/**
- * Overload for array of x and mu and non-array sigma with lp
- */
-template <typename T, typename M, typename S,
-          require_not_std_vector_t<S>* = nullptr>
-inline auto offset_multiplier_constrain(const std::vector<T>& x,
-                                        const std::vector<M>& mu,
-                                        const S& sigma,
-                                        return_type_t<T, M, S>& lp) {
-  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
-  std::vector<plain_type_t<decltype(
-      offset_multiplier_constrain(x[0], mu[0], sigma, lp))>>
-      ret;
-  ret.reserve(x.size());
-  const auto& sigma_ref = to_ref(sigma);
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma_ref, lp));
-  }
-  return ret;
-}
-
-/**
- * Overload for array of x, mu, and sigma
- */
-template <typename T, typename M, typename S>
-inline auto offset_multiplier_constrain(const std::vector<T>& x,
-                                        const std::vector<M>& mu,
-                                        const std::vector<S>& sigma) {
-  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
-  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
-  std::vector<plain_type_t<decltype(
-      offset_multiplier_constrain(x[0], mu[0], sigma[0]))>>
-      ret;
-  ret.reserve(x.size());
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma[i]));
-  }
-  return ret;
-}
-
-/**
- * Overload for array of x, mu, and sigma with lp
- */
 template <typename T, typename M, typename S>
-inline auto offset_multiplier_constrain(const std::vector<T>& x,
-                                        const std::vector<M>& mu,
-                                        const std::vector<S>& sigma,
-                                        return_type_t<T, M, S>& lp) {
-  check_matching_dims("offset_multiplier_constrain", "x", x, "mu", mu);
-  check_matching_dims("offset_multiplier_constrain", "x", x, "sigma", sigma);
-  std::vector<plain_type_t<decltype(
-      offset_multiplier_constrain(x[0], mu[0], sigma[0], lp))>>
-      ret;
-  ret.reserve(x.size());
-  for (size_t i = 0; i < x.size(); ++i) {
-    ret.emplace_back(offset_multiplier_constrain(x[i], mu[i], sigma[i], lp));
-  }
-  return ret;
-}
-
-/**
- * Return the linearly transformed value for the specified unconstrained input
- * and specified offset and multiplier. If the `Jacobian` parameter is `true`,
- * the log density accumulator is incremented with the log absolute Jacobian
- * determinant of the transform.  All of the transforms are specified with their
- * Jacobians in the *Stan Reference Manual* chapter Constraint Transforms.
- *
- * @tparam Jacobian if `true`, increment log density accumulator with log
- * absolute Jacobian determinant of constraining transform
- * @tparam T A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
- * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
- * @tparam M A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
- * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
- * @tparam S A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
- * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
- * @param[in] x Unconstrained scalar input
- * @param[in] mu offset of constrained output
- * @param[in] sigma multiplier of constrained output
- * @param[in, out] lp log density accumulator
- * @return linear transformed value corresponding to inputs
- * @throw std::domain_error if sigma <= 0
- * @throw std::domain_error if mu is not finite
- */
-template <bool Jacobian, typename T, typename M, typename S>
-inline auto offset_multiplier_constrain(const T& x, const M& mu, const S& sigma,
-                                        return_type_t<T, M, S>& lp) {
-  if (Jacobian) {
-    return offset_multiplier_constrain(x, mu, sigma, lp);
-  } else {
-    return offset_multiplier_constrain(x, mu, sigma);
+inline return_type_t<T, M, S> offset_multiplier_constrain(const T& x,
+                                                          const M& mu,
+                                                          const S& sigma,
+                                                          T& lp) {
+  using std::log;
+  check_finite("offset_multiplier_constrain", "offset", mu);
+  if (sigma == 1) {
+    if (mu == 0) {
+      return identity_constrain(x);
+    }
+    return mu + x;
   }
+  check_positive_finite("offset_multiplier_constrain", "multiplier", sigma);
+  lp += multiply_log(math::size(x), sigma);
+  return fma(sigma, x, mu);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/prob/bernoulli_lpmf.hpp b/stan/math/prim/prob/bernoulli_lpmf.hpp
index 9b54f02d839..c288a3ddd48 100644
--- a/stan/math/prim/prob/bernoulli_lpmf.hpp
+++ b/stan/math/prim/prob/bernoulli_lpmf.hpp
@@ -42,8 +42,7 @@ return_type_t<T_prob> bernoulli_lpmf(const T_n& n, const T_prob& theta) {
   const T_n_ref n_ref = to_ref(n);
   const T_theta_ref theta_ref = to_ref(theta);
   check_bounded(function, "n", n_ref, 0, 1);
-  check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0,
-                1.0);
+  check_bounded(function, "Probability parameter", theta_ref, 0.0, 1.0);
 
   if (size_zero(n, theta)) {
     return 0.0;
@@ -84,8 +83,8 @@ return_type_t<T_prob> bernoulli_lpmf(const T_n& n, const T_prob& theta) {
       logp += (N - sum) * log1m_theta;
 
       if (!is_constant_all<T_prob>::value) {
-        ops_partials.edge1_.partials_[0] += sum * inv(theta_dbl);
-        ops_partials.edge1_.partials_[0] += (N - sum) * inv(theta_dbl - 1);
+        ops_partials.edge1_.partials_[0] += sum / theta_dbl;
+        ops_partials.edge1_.partials_[0] += (N - sum) / (theta_dbl - 1);
       }
     }
   } else {
diff --git a/stan/math/prim/prob/beta_proportion_lpdf.hpp b/stan/math/prim/prob/beta_proportion_lpdf.hpp
index 523fdc1c138..73aa67acab0 100644
--- a/stan/math/prim/prob/beta_proportion_lpdf.hpp
+++ b/stan/math/prim/prob/beta_proportion_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
@@ -52,6 +51,8 @@ return_type_t<T_y, T_loc, T_prec> beta_proportion_lpdf(const T_y& y,
                                                        const T_loc& mu,
                                                        const T_prec& kappa) {
   using T_partials_return = partials_return_t<T_y, T_loc, T_prec>;
+  using T_partials_return_kappa = return_type_t<T_prec>;
+  using T_partials_array = Eigen::Array<T_partials_return, Eigen::Dynamic, 1>;
   using std::log;
   using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
   using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
@@ -67,14 +68,23 @@ return_type_t<T_y, T_loc, T_prec> beta_proportion_lpdf(const T_y& y,
   T_mu_ref mu_ref = mu;
   T_kappa_ref kappa_ref = kappa;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
-  decltype(auto) kappa_val = to_ref(as_value_column_array_or_scalar(kappa_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
+  const auto& kappa_col = as_column_vector_or_scalar(kappa_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& mu_arr = as_array_or_scalar(mu_col);
+  const auto& kappa_arr
+      = promote_scalar<T_partials_return_kappa>(as_array_or_scalar(kappa_col));
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
+  ref_type_t<decltype(value_of(kappa_arr))> kappa_val = value_of(kappa_arr);
 
   check_positive(function, "Location parameter", mu_val);
   check_less(function, "Location parameter", mu_val, 1.0);
   check_positive_finite(function, "Precision parameter", kappa_val);
-  check_bounded(function, "Random variable", value_of(y_val), 0, 1);
+  check_bounded(function, "Random variable", y_val, 0, 1);
 
   if (!include_summand<propto, T_y, T_loc, T_prec>::value) {
     return 0;
diff --git a/stan/math/prim/prob/cauchy_lpdf.hpp b/stan/math/prim/prob/cauchy_lpdf.hpp
index 6e16581e786..302c3c22a9b 100644
--- a/stan/math/prim/prob/cauchy_lpdf.hpp
+++ b/stan/math/prim/prob/cauchy_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/log1p.hpp>
@@ -43,6 +42,7 @@ template <bool propto, typename T_y, typename T_loc, typename T_scale,
 return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
                                                const T_scale& sigma) {
   using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
+  using T_partials_array = Eigen::Array<T_partials_return, Eigen::Dynamic, 1>;
   using std::log;
   using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
   using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
@@ -65,9 +65,17 @@ return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
   operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref> ops_partials(
       y_ref, mu_ref, sigma_ref);
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
-  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
+  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& mu_arr = as_array_or_scalar(mu_col);
+  const auto& sigma_arr = as_array_or_scalar(sigma_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
+  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
   check_not_nan(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
   check_positive_finite(function, "Scale parameter", sigma_val);
@@ -93,8 +101,8 @@ return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
     const auto& y_minus_mu_squared
         = to_ref_if<!is_constant_all<T_scale>::value>(square(y_minus_mu));
     if (!is_constant_all<T_y, T_loc>::value) {
-      auto mu_deriv = to_ref_if<(!is_constant_all<T_y>::value
-                                 && !is_constant_all<T_loc>::value)>(
+      const auto& mu_deriv = to_ref_if<(!is_constant_all<T_y>::value
+                                        && !is_constant_all<T_loc>::value)>(
           2 * y_minus_mu / (sigma_squared + y_minus_mu_squared));
       if (!is_constant_all<T_y>::value) {
         if (is_vector<T_y>::value) {
diff --git a/stan/math/prim/prob/chi_square_lpdf.hpp b/stan/math/prim/prob/chi_square_lpdf.hpp
index 885d9a3c988..f96260950ec 100644
--- a/stan/math/prim/prob/chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/chi_square_lpdf.hpp
@@ -4,7 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
+#include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
@@ -54,10 +54,17 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {
                          "Degrees of freedom parameter", nu);
   T_y_ref y_ref = y;
   T_nu_ref nu_ref = nu;
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& nu_col = as_column_vector_or_scalar(nu_ref);
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) nu_val = to_ref(as_value_column_array_or_scalar(nu_ref));
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& nu_arr = as_array_or_scalar(nu_col);
 
+  ref_type_if_t<include_summand<propto, T_y>::value, decltype(value_of(y_arr))>
+      y_val = value_of(y_arr);
+  ref_type_if_t<include_summand<propto, T_dof>::value,
+                decltype(value_of(nu_arr))>
+      nu_val = value_of(nu_arr);
   check_nonnegative(function, "Random variable", y_val);
   check_positive_finite(function, "Degrees of freedom parameter", nu_val);
 
diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp
index 3db919362d2..832cafa50d8 100644
--- a/stan/math/prim/prob/double_exponential_lpdf.hpp
+++ b/stan/math/prim/prob/double_exponential_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/fabs.hpp>
 #include <stan/math/prim/fun/inv.hpp>
@@ -42,6 +41,7 @@ template <bool propto, typename T_y, typename T_loc, typename T_scale,
 return_type_t<T_y, T_loc, T_scale> double_exponential_lpdf(
     const T_y& y, const T_loc& mu, const T_scale& sigma) {
   using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
+  using T_partials_array = Eigen::Array<T_partials_return, Eigen::Dynamic, 1>;
   using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
   using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
   using T_sigma_ref = ref_type_if_t<!is_constant<T_scale>::value, T_scale>;
@@ -63,9 +63,17 @@ return_type_t<T_y, T_loc, T_scale> double_exponential_lpdf(
   operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref> ops_partials(
       y_ref, mu_ref, sigma_ref);
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
-  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
+  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& mu_arr = as_array_or_scalar(mu_col);
+  const auto& sigma_arr = as_array_or_scalar(sigma_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
+  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
 
   check_finite(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
diff --git a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
index ba3eddbcbde..0ee1f629f11 100644
--- a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/erfc.hpp>
 #include <stan/math/prim/fun/exp.hpp>
@@ -44,11 +43,20 @@ return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_lpdf(
   T_sigma_ref sigma_ref = sigma;
   T_lambda_ref lambda_ref = lambda;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
-  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
-  decltype(auto) lambda_val
-      = to_ref(as_value_column_array_or_scalar(lambda_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
+  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
+  const auto& lambda_col = as_column_vector_or_scalar(lambda_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& mu_arr = as_array_or_scalar(mu_col);
+  const auto& sigma_arr = as_array_or_scalar(sigma_col);
+  const auto& lambda_arr = as_array_or_scalar(lambda_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
+  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
+  ref_type_t<decltype(value_of(lambda_arr))> lambda_val = value_of(lambda_arr);
 
   check_not_nan(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);
diff --git a/stan/math/prim/prob/exponential_lccdf.hpp b/stan/math/prim/prob/exponential_lccdf.hpp
index c48c90ab9d3..f1920bbb2dc 100644
--- a/stan/math/prim/prob/exponential_lccdf.hpp
+++ b/stan/math/prim/prob/exponential_lccdf.hpp
@@ -3,8 +3,8 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
@@ -14,9 +14,7 @@
 namespace stan {
 namespace math {
 
-template <typename T_y, typename T_inv_scale,
-          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
-              T_y, T_inv_scale>* = nullptr>
+template <typename T_y, typename T_inv_scale>
 return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
                                                   const T_inv_scale& beta) {
   using T_partials_return = partials_return_t<T_y, T_inv_scale>;
@@ -28,8 +26,14 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
   T_y_ref y_ref = y;
   T_beta_ref beta_ref = beta;
 
-  auto y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  auto beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& beta_arr = as_array_or_scalar(beta_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
 
   check_nonnegative(function, "Random variable", y_val);
   check_positive_finite(function, "Inverse scale parameter", beta_val);
@@ -51,7 +55,9 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     } else if (is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = -forward_as<beta_val_array>(beta_val);
     } else {
-      ops_partials.edge1_.partials_[0] = -sum(beta_val);
+      forward_as<internal::broadcast_array<T_partials_return>>(
+          ops_partials.edge1_.partials_)
+          = -forward_as<beta_val_scalar>(beta_val);
     }
   }
   if (!is_constant_all<T_inv_scale>::value) {
@@ -63,7 +69,9 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     } else if (is_vector<T_y>::value) {
       ops_partials.edge2_.partials_ = -forward_as<y_val_array>(y_val);
     } else {
-      ops_partials.edge2_.partials_[0] = -sum(y_val);
+      forward_as<internal::broadcast_array<T_partials_return>>(
+          ops_partials.edge2_.partials_)
+          = -forward_as<y_val_scalar>(y_val);
     }
   }
   return ops_partials.build(ccdf_log);
diff --git a/stan/math/prim/prob/exponential_lpdf.hpp b/stan/math/prim/prob/exponential_lpdf.hpp
index 18f3bbfd0da..2cbb4e7e06a 100644
--- a/stan/math/prim/prob/exponential_lpdf.hpp
+++ b/stan/math/prim/prob/exponential_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
@@ -62,8 +61,14 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
   T_y_ref y_ref = y;
   T_beta_ref beta_ref = beta;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& beta_arr = as_array_or_scalar(beta_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
 
   check_nonnegative(function, "Random variable", y_val);
   check_positive_finite(function, "Inverse scale parameter", beta_val);
@@ -76,7 +81,7 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
 
   T_partials_return logp(0.0);
   if (include_summand<propto, T_inv_scale>::value) {
-    logp = sum(log(beta_val)) * max_size(y, beta) / math::size(beta);
+    logp = sum(log(beta_val)) * max_size(y, beta) / size(beta);
   }
   if (include_summand<propto, T_y, T_inv_scale>::value) {
     logp -= sum(beta_val * y_val);
diff --git a/stan/math/prim/prob/frechet_lpdf.hpp b/stan/math/prim/prob/frechet_lpdf.hpp
index b305f0a8f20..56bc37b2b5a 100644
--- a/stan/math/prim/prob/frechet_lpdf.hpp
+++ b/stan/math/prim/prob/frechet_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
@@ -43,9 +42,17 @@ return_type_t<T_y, T_shape, T_scale> frechet_lpdf(const T_y& y,
   T_sigma_ref sigma_ref = sigma;
   using std::pow;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
-  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
+  const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& alpha_arr = as_array_or_scalar(alpha_col);
+  const auto& sigma_arr = as_array_or_scalar(sigma_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
+  ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
 
   check_positive(function, "Random variable", y_val);
   check_positive_finite(function, "Shape parameter", alpha_val);
diff --git a/stan/math/prim/prob/gamma_lpdf.hpp b/stan/math/prim/prob/gamma_lpdf.hpp
index 50e7a6b65c8..05aa0a13d3e 100644
--- a/stan/math/prim/prob/gamma_lpdf.hpp
+++ b/stan/math/prim/prob/gamma_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
@@ -64,9 +63,17 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lpdf(const T_y& y,
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
-  decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
+  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& alpha_arr = as_array_or_scalar(alpha_col);
+  const auto& beta_arr = as_array_or_scalar(beta_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
+  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
 
   check_not_nan(function, "Random variable", y_val);
   check_positive_finite(function, "Shape parameter", alpha_val);
diff --git a/stan/math/prim/prob/gumbel_lpdf.hpp b/stan/math/prim/prob/gumbel_lpdf.hpp
index 22814b8982b..ad02f4ac088 100644
--- a/stan/math/prim/prob/gumbel_lpdf.hpp
+++ b/stan/math/prim/prob/gumbel_lpdf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
@@ -50,9 +49,17 @@ return_type_t<T_y, T_loc, T_scale> gumbel_lpdf(const T_y& y, const T_loc& mu,
   T_mu_ref mu_ref = mu;
   T_beta_ref beta_ref = beta;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
-  decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
+  const auto& y_col = as_column_vector_or_scalar(y_ref);
+  const auto& mu_col = as_column_vector_or_scalar(mu_ref);
+  const auto& beta_col = as_column_vector_or_scalar(beta_ref);
+
+  const auto& y_arr = as_array_or_scalar(y_col);
+  const auto& mu_arr = as_array_or_scalar(mu_col);
+  const auto& beta_arr = as_array_or_scalar(beta_col);
+
+  ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
+  ref_type_t<decltype(value_of(mu_arr))> mu_val = value_of(mu_arr);
+  ref_type_t<decltype(value_of(beta_arr))> beta_val = value_of(beta_arr);
 
   check_not_nan(function, "Random variable", y_val);
   check_finite(function, "Location parameter", mu_val);

From 33b6e7755ebae2878dcf27ce59755abcdcc7c7f6 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 16:15:10 +0200
Subject: [PATCH 6/9] Missed parentheses

---
 stan/math/opencl/matrix_cl.hpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/stan/math/opencl/matrix_cl.hpp b/stan/math/opencl/matrix_cl.hpp
index 561694cda6f..5248d719b93 100644
--- a/stan/math/opencl/matrix_cl.hpp
+++ b/stan/math/opencl/matrix_cl.hpp
@@ -272,7 +272,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   matrix_cl(const int rows, const int cols,
             matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(rows), cols_(cols), view_(partial_view) {
-    if this->size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     cl::Context& ctx = opencl_context.context();
@@ -310,7 +310,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(A.rows()), cols_(A.cols()), view_(partial_view) {
     using Mat_type = std::decay_t<ref_type_for_opencl_t<Mat>>;
-    if this->size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer_no_heap_if<
@@ -489,7 +489,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(const T* A) {
     cl::Event transfer_event;
-    if this->size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -509,7 +509,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(T* A) {
     cl::Event transfer_event;
-    if this->size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -548,7 +548,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
    */
   template <bool No_heap, typename U, std::enable_if_t<No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
-    if this->size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer(obj.data());
@@ -558,7 +558,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool No_heap, typename U, std::enable_if_t<!No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
     using U_val = std::decay_t<ref_type_for_opencl_t<U>>;
-    if this->size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     auto* obj_heap = new U_val(std::move(obj));

From cf46f37754bb1ad7933b4ccbcebea3def92e5f91 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 18:45:00 +0200
Subject: [PATCH 7/9] Mismatched parentheses

---
 stan/math/rev/functor/coupled_ode_system.hpp | 6 +++---
 stan/math/rev/functor/cvodes_integrator.hpp  | 8 ++++----
 stan/math/rev/functor/ode_adams.hpp          | 2 +-
 stan/math/rev/functor/ode_bdf.hpp            | 2 +-
 stan/math/rev/functor/reduce_sum.hpp         | 6 +++---
 5 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/stan/math/rev/functor/coupled_ode_system.hpp b/stan/math/rev/functor/coupled_ode_system.hpp
index 3ea68004a65..33c0596a0a7 100644
--- a/stan/math/rev/functor/coupled_ode_system.hpp
+++ b/stan/math/rev/functor/coupled_ode_system.hpp
@@ -126,7 +126,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
       y_vars.coeffRef(n) = z[n];
 
     Eigen::Matrix<var, Eigen::Dynamic, 1> f_y_t_vars
-        = math::apply(([&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
+        = math::apply([&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
                 local_args_tuple_);
 
     check_size_match("coupled_ode_system", "dy_dt", f_y_t_vars.size(), "states",
@@ -141,7 +141,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
       // memset was faster than Eigen setZero
       memset(args_adjoints_.data(), 0, sizeof(double) * num_args_vars);
 
-      math::apply((
+      math::apply(
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },
@@ -149,7 +149,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
 
       // The vars here do not live on the nested stack so must be zero'd
       // separately
-      math::apply(([&](auto&&... args) { zero_adjoints(args...); }, local_args_tuple_);
+      math::apply([&](auto&&... args) { zero_adjoints(args...); }, local_args_tuple_);
 
       // No need to zero adjoints after last sweep
       if (i + 1 < N_) {
diff --git a/stan/math/rev/functor/cvodes_integrator.hpp b/stan/math/rev/functor/cvodes_integrator.hpp
index a174205d238..83c395160c5 100644
--- a/stan/math/rev/functor/cvodes_integrator.hpp
+++ b/stan/math/rev/functor/cvodes_integrator.hpp
@@ -104,7 +104,7 @@ class cvodes_integrator {
     const Eigen::VectorXd y_vec = Eigen::Map<const Eigen::VectorXd>(y, N_);
 
     Eigen::VectorXd dy_dt_vec
-        = math::apply(([&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
+        = math::apply([&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
                 value_of_args_tuple_);
 
     check_size_match("cvodes_integrator", "dy_dt", dy_dt_vec.size(), "states",
@@ -122,7 +122,7 @@ class cvodes_integrator {
     Eigen::MatrixXd Jfy;
 
     auto f_wrapped = [&](const Eigen::Matrix<var, Eigen::Dynamic, 1>& y) {
-      return math::apply(([&](auto&&... args) { return f_(t, y, msgs_, args...); },
+      return math::apply([&](auto&&... args) { return f_(t, y, msgs_, args...); },
                    value_of_args_tuple_);
     };
 
@@ -212,7 +212,7 @@ class cvodes_integrator {
 
     // Code from: https://stackoverflow.com/a/17340003 . Should probably do
     // something better
-    math::apply((
+    math::apply(
         [&](auto&&... args) {
           std::vector<int> unused_temp{
               0, (check_finite(function_name, "ode parameters and data", args),
@@ -327,7 +327,7 @@ class cvodes_integrator {
           }
         }
 
-        y.emplace_back(math::apply((
+        y.emplace_back(math::apply(
             [&](auto&&... args) {
               return ode_store_sensitivities(f_, coupled_state_, y0_, t0_,
                                              ts_[n], msgs_, args...);
diff --git a/stan/math/rev/functor/ode_adams.hpp b/stan/math/rev/functor/ode_adams.hpp
index 204adab6bce..500cd2e4c92 100644
--- a/stan/math/rev/functor/ode_adams.hpp
+++ b/stan/math/rev/functor/ode_adams.hpp
@@ -55,7 +55,7 @@ ode_adams_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                    long int max_num_steps,  // NOLINT(runtime/int)
                    std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return math::apply((
+  return math::apply(
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_ADAMS, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/ode_bdf.hpp b/stan/math/rev/functor/ode_bdf.hpp
index 75c8b6bc197..8042f848b69 100644
--- a/stan/math/rev/functor/ode_bdf.hpp
+++ b/stan/math/rev/functor/ode_bdf.hpp
@@ -56,7 +56,7 @@ ode_bdf_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                  long int max_num_steps,  // NOLINT(runtime/int)
                  std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return math::apply((
+  return math::apply(
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_BDF, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/reduce_sum.hpp b/stan/math/rev/functor/reduce_sum.hpp
index bfecd4017e2..2ee3927618a 100644
--- a/stan/math/rev/functor/reduce_sum.hpp
+++ b/stan/math/rev/functor/reduce_sum.hpp
@@ -113,7 +113,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
         // scope. In this case no need for zeroing adjoints, since the
         // fresh copy has all adjoints set to zero.
         local_args_tuple_scope_.stack_.execute([&]() {
-          math::apply((
+          math::apply(
               [&](auto&&... args) {
                 local_args_tuple_scope_.args_tuple_holder_ = std::make_unique<
                     typename scoped_args_tuple::args_tuple_t>(
@@ -140,7 +140,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
       }
 
       // Perform calculation
-      var sub_sum_v = math::apply((
+      var sub_sum_v = math::apply(
           [&](auto&&... args) {
             return ReduceFunction()(local_sub_slice, r.begin(), r.end() - 1,
                                     &msgs_, args...);
@@ -158,7 +158,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
                           std::move(local_sub_slice));
 
       // Accumulate adjoints of shared_arguments
-      math::apply((
+      math::apply(
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },

From 89139c65bcc8f59d2800a262d220a93f93274edf Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Wed, 8 Feb 2023 19:58:26 +0200
Subject: [PATCH 8/9] LDLT path wrong

---
 stan/math/prim/mat/fun/LDLT_factor.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stan/math/prim/mat/fun/LDLT_factor.hpp b/stan/math/prim/mat/fun/LDLT_factor.hpp
index ed6f1121c85..f782f289e9b 100644
--- a/stan/math/prim/mat/fun/LDLT_factor.hpp
+++ b/stan/math/prim/mat/fun/LDLT_factor.hpp
@@ -1 +1 @@
-#include "../prim/fun/LDLT_factor.hpp"
+#include "../../fun/LDLT_factor.hpp"

From a80d1e25906150284f7d14b4d3ffec203ca8e407 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Thu, 9 Feb 2023 07:02:35 +0200
Subject: [PATCH 9/9] Revert another stray change

---
 stan/math/prim/prob/binomial_logit_lpmf.hpp | 24 ++++++++++++---------
 1 file changed, 14 insertions(+), 10 deletions(-)

diff --git a/stan/math/prim/prob/binomial_logit_lpmf.hpp b/stan/math/prim/prob/binomial_logit_lpmf.hpp
index 08e42e2873a..816d06a75a7 100644
--- a/stan/math/prim/prob/binomial_logit_lpmf.hpp
+++ b/stan/math/prim/prob/binomial_logit_lpmf.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/binomial_coefficient_log.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
@@ -34,9 +33,7 @@ namespace math {
  * @throw std::domain_error if N is negative or probability parameter is invalid
  * @throw std::invalid_argument if vector sizes do not match
  */
-template <bool propto, typename T_n, typename T_N, typename T_prob,
-          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
-              T_n, T_N, T_prob>* = nullptr>
+template <bool propto, typename T_n, typename T_N, typename T_prob>
 return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
                                           const T_prob& alpha) {
   using T_partials_return = partials_return_t<T_n, T_N, T_prob>;
@@ -52,11 +49,19 @@ return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
   T_N_ref N_ref = N;
   T_alpha_ref alpha_ref = alpha;
 
-  decltype(auto) n_val = to_ref(as_value_column_array_or_scalar(n_ref));
-  decltype(auto) N_val = to_ref(as_value_column_array_or_scalar(N_ref));
-  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
+  const auto& n_col = as_column_vector_or_scalar(n_ref);
+  const auto& N_col = as_column_vector_or_scalar(N_ref);
+  const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
 
-  check_bounded(function, "Successes variable", value_of(n_val), 0, N_val);
+  const auto& n_arr = as_array_or_scalar(n_col);
+  const auto& N_arr = as_array_or_scalar(N_col);
+  const auto& alpha_arr = as_array_or_scalar(alpha_col);
+
+  ref_type_t<decltype(value_of(n_arr))> n_val = value_of(n_arr);
+  ref_type_t<decltype(value_of(N_arr))> N_val = value_of(N_arr);
+  ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
+
+  check_bounded(function, "Successes variable", n_val, 0, N_val);
   check_nonnegative(function, "Population size parameter", N_val);
   check_finite(function, "Probability parameter", alpha_val);
 
@@ -90,8 +95,7 @@ return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
       T_partials_return sum_n = sum(n_val) * maximum_size / math::size(n);
       ops_partials.edge1_.partials_[0] = forward_as<T_partials_return>(
           sum_n * inv_logit_neg_alpha
-          - (sum(N_val) * maximum_size / math::size(N) - sum_n)
-                * inv_logit_alpha);
+          - (sum(N_val) * maximum_size / math::size(N) - sum_n) * inv_logit_alpha);
     }
   }