From 953a35d6092535d6cb6a68d649ebfe2e85363080 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Mon, 14 Jan 2019 21:51:54 +0100
Subject: [PATCH 01/20] make thread-specific AD tape pointer much faster (on
 Mac)

---
 make/libraries                              |  5 ++++
 stan/math/rev/core/autodiffstackstorage.hpp | 29 ++++++++++++---------
 2 files changed, 22 insertions(+), 12 deletions(-)

diff --git a/make/libraries b/make/libraries
index 0140776c9b0..6d08cce9f8e 100644
--- a/make/libraries
+++ b/make/libraries
@@ -110,6 +110,11 @@ clean-mpi:
 
 endif
 
+STACK_TEMPLATE_INSTANTIATION_CPP := $(shell find $(MATH)stan -type f -name 'chainablestack_inst.cpp')
+STACK_TEMPLATE_INSTANTIATION := $(STACK_TEMPLATE_INSTANTIATION_CPP:%.cpp=%.o)
+
+$(STACK_TEMPLATE_INSTANTIATION) : CXXFLAGS += -fPIC
+
 ############################################################
 # Google Test:
 #   Build the google test library.
diff --git a/stan/math/rev/core/autodiffstackstorage.hpp b/stan/math/rev/core/autodiffstackstorage.hpp
index 1b08e3cd165..0d7d3fafc43 100644
--- a/stan/math/rev/core/autodiffstackstorage.hpp
+++ b/stan/math/rev/core/autodiffstackstorage.hpp
@@ -61,26 +61,31 @@ struct AutodiffStackSingleton {
   explicit AutodiffStackSingleton(AutodiffStackSingleton_t const &) = delete;
   AutodiffStackSingleton &operator=(const AutodiffStackSingleton_t &) = delete;
 
-  static inline AutodiffStackStorage &instance() {
-#ifdef STAN_THREADS
-    thread_local static AutodiffStackStorage instance_;
-#endif
-    return instance_;
+  constexpr static inline AutodiffStackStorage &instance() {
+    return *instance_;
   }
 
-#ifndef STAN_THREADS
-
+  static inline AutodiffStackStorage* instantiate() {
+    if(instance_ == 0)
+      instance_ = new AutodiffStackStorage();
+    return instance_;
+  }
+  
  private:
-  static AutodiffStackStorage instance_;
+  static
+#ifdef STAN_THREADS
+  thread_local
 #endif
+  AutodiffStackStorage* instance_;
 };
 
-#ifndef STAN_THREADS
 template <typename ChainableT, typename ChainableAllocT>
-typename AutodiffStackSingleton<ChainableT,
-                                ChainableAllocT>::AutodiffStackStorage
-    AutodiffStackSingleton<ChainableT, ChainableAllocT>::instance_;
+#ifdef STAN_THREADS
+thread_local
 #endif
+typename AutodiffStackSingleton<ChainableT,
+                                ChainableAllocT>::AutodiffStackStorage*
+AutodiffStackSingleton<ChainableT, ChainableAllocT>::instance_ = 0;
 
 }  // namespace math
 }  // namespace stan

From 9179928f8f006d982a5f69ba93f2372dbce56429 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sebastian.weber@novartis.com>
Date: Tue, 15 Jan 2019 17:50:58 +0100
Subject: [PATCH 02/20] add makefile bits

---
 make/libraries                             |  6 +++++-
 stan/math/rev/core/chainablestack_inst.cpp | 15 +++++++++++++++
 2 files changed, 20 insertions(+), 1 deletion(-)
 create mode 100644 stan/math/rev/core/chainablestack_inst.cpp

diff --git a/make/libraries b/make/libraries
index 6d08cce9f8e..19211898422 100644
--- a/make/libraries
+++ b/make/libraries
@@ -115,6 +115,10 @@ STACK_TEMPLATE_INSTANTIATION := $(STACK_TEMPLATE_INSTANTIATION_CPP:%.cpp=%.o)
 
 $(STACK_TEMPLATE_INSTANTIATION) : CXXFLAGS += -fPIC
 
+clean-math:
+	@echo '  cleaning stan-math targets'
+	$(RM) $(STACK_TEMPLATE_INSTANTIATION)
+
 ############################################################
 # Google Test:
 #   Build the google test library.
@@ -128,4 +132,4 @@ $(GTEST)/src/gtest-all.o: INC += $(INC_GTEST)
 # Clean all libraries
 
 .PHONY: clean-libraries clean-sundials clean-mpi
-clean-libraries: clean-sundials clean-mpi
+clean-libraries: clean-sundials clean-mpi clean-math
diff --git a/stan/math/rev/core/chainablestack_inst.cpp b/stan/math/rev/core/chainablestack_inst.cpp
new file mode 100644
index 00000000000..210a66cd8f6
--- /dev/null
+++ b/stan/math/rev/core/chainablestack_inst.cpp
@@ -0,0 +1,15 @@
+#ifndef STAN_MATH_REV_CORE_CHAINABLESTACK_INST_CPP
+#define STAN_MATH_REV_CORE_CHAINABLESTACK_INST_CPP
+
+#include <stan/math/rev/core/chainablestack.hpp>
+
+namespace stan {
+namespace math {
+namespace internal {
+
+  const ChainableStack::AutodiffStackStorage* __chainable_stack = ChainableStack::instantiate();
+
+}
+}  // namespace math
+}  // namespace stan
+#endif

From 27e9ee9f0e0e07cff5c300ad877a8e8c0cae6139 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Tue, 15 Jan 2019 21:54:39 +0100
Subject: [PATCH 03/20] switch to new thread_local scheme

---
 make/libraries                                |  8 +++----
 make/tests                                    |  2 +-
 stan/math/rev/core/autodiffstackstorage.hpp   |  3 +++
 stan/math/rev/mat/functor/map_rect_reduce.hpp |  2 ++
 .../math/rev/mat/functor/gradient_test.cpp    | 21 ++++++++++---------
 5 files changed, 21 insertions(+), 15 deletions(-)

diff --git a/make/libraries b/make/libraries
index 19211898422..01f2dc32bf0 100644
--- a/make/libraries
+++ b/make/libraries
@@ -110,14 +110,14 @@ clean-mpi:
 
 endif
 
-STACK_TEMPLATE_INSTANTIATION_CPP := $(shell find $(MATH)stan -type f -name 'chainablestack_inst.cpp')
-STACK_TEMPLATE_INSTANTIATION := $(STACK_TEMPLATE_INSTANTIATION_CPP:%.cpp=%.o)
+STACK_INSTANTIATION_CPP := $(shell find $(MATH)stan -type f -name 'chainablestack_inst.cpp')
+STACK_INSTANTIATION := $(STACK_INSTANTIATION_CPP:%.cpp=%.o)
 
-$(STACK_TEMPLATE_INSTANTIATION) : CXXFLAGS += -fPIC
+$(STACK_INSTANTIATION) : CXXFLAGS += -fPIC
 
 clean-math:
 	@echo '  cleaning stan-math targets'
-	$(RM) $(STACK_TEMPLATE_INSTANTIATION)
+	$(RM) $(STACK_INSTANTIATION)
 
 ############################################################
 # Google Test:
diff --git a/make/tests b/make/tests
index ba1aeacfe5d..6126b3b89c3 100644
--- a/make/tests
+++ b/make/tests
@@ -9,7 +9,7 @@ test/% : CXXFLAGS += $(CXXFLAGS_GTEST)
 test/% : CPPFLAGS += $(CPPFLAGS_GTEST)
 test/% : INC += $(INC_GTEST)
 
-test/%$(EXE) : test/%.o $(GTEST)/src/gtest_main.cc $(GTEST)/src/gtest-all.o $(MPI_TARGETS)
+test/%$(EXE) : test/%.o $(GTEST)/src/gtest_main.cc $(GTEST)/src/gtest-all.o $(MPI_TARGETS) $(STACK_INSTANTIATION)
 	$(LINK.cpp) $^ $(LDLIBS) $(OUTPUT_OPTION)
 
 ##
diff --git a/stan/math/rev/core/autodiffstackstorage.hpp b/stan/math/rev/core/autodiffstackstorage.hpp
index 0d7d3fafc43..58d8140830d 100644
--- a/stan/math/rev/core/autodiffstackstorage.hpp
+++ b/stan/math/rev/core/autodiffstackstorage.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_REV_CORE_AUTODIFFSTACKSTORAGE_HPP
 
 #include <stan/math/memory/stack_alloc.hpp>
+#include <mutex>
 #include <vector>
 
 namespace stan {
@@ -66,6 +67,8 @@ struct AutodiffStackSingleton {
   }
 
   static inline AutodiffStackStorage* instantiate() {
+    static thread_local std::mutex init_mutex;
+    std::lock_guard<std::mutex> init_lock(init_mutex);
     if(instance_ == 0)
       instance_ = new AutodiffStackStorage();
     return instance_;
diff --git a/stan/math/rev/mat/functor/map_rect_reduce.hpp b/stan/math/rev/mat/functor/map_rect_reduce.hpp
index 4d9fe49d030..1b1ca9dce72 100644
--- a/stan/math/rev/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/rev/mat/functor/map_rect_reduce.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/mat/functor/map_rect_reduce.hpp>
 #include <stan/math/prim/mat/fun/typedefs.hpp>
+#include <stan/math/rev/core/chainablestack.hpp>
 #include <stan/math/rev/core/var.hpp>
 #include <stan/math/rev/mat/fun/typedefs.hpp>
 #include <stan/math/rev/mat/fun/to_var.hpp>
@@ -23,6 +24,7 @@ struct map_rect_reduce<F, var, var> {
     const size_type num_shared_params = shared_params.rows();
     const size_type num_job_specific_params = job_specific_params.rows();
     matrix_d out(1 + num_shared_params + num_job_specific_params, 0);
+    ChainableStack::instantiate();
 
     try {
       start_nested();
diff --git a/test/unit/math/rev/mat/functor/gradient_test.cpp b/test/unit/math/rev/mat/functor/gradient_test.cpp
index 0321e921785..28c50e63717 100644
--- a/test/unit/math/rev/mat/functor/gradient_test.cpp
+++ b/test/unit/math/rev/mat/functor/gradient_test.cpp
@@ -46,16 +46,17 @@ TEST(AgradAutoDiff, gradient_threaded) {
   EXPECT_FLOAT_EQ(x_ref(0) * x_ref(0) + 3 * 2 * x_ref(1), grad_fx_ref(1));
 
   auto thread_job = [&](double x1, double x2) {
-    double fx;
-    VectorXd x_local(2);
-    x_local << x1, x2;
-    VectorXd grad_fx;
-    stan::math::gradient(fun1(), x_local, fx, grad_fx);
-    VectorXd res(1 + grad_fx.size());
-    res(0) = fx;
-    res.tail(grad_fx.size()) = grad_fx;
-    return res;
-  };
+                      stan::math::ChainableStack::instantiate();
+                      double fx;
+                      VectorXd x_local(2);
+                      x_local << x1, x2;
+                      VectorXd grad_fx;
+                      stan::math::gradient(fun1(), x_local, fx, grad_fx);
+                      VectorXd res(1 + grad_fx.size());
+                      res(0) = fx;
+                      res.tail(grad_fx.size()) = grad_fx;
+                      return res;
+                    };
 
   // schedule a bunch of jobs which all do the same
   std::vector<std::future<VectorXd>> ad_futures_ref;

From 9cc04b000d1bb634d383f589be5846ccc198a176 Mon Sep 17 00:00:00 2001
From: Stan Jenkins <mc.stanislaw@gmail.com>
Date: Tue, 15 Jan 2019 21:15:18 +0000
Subject: [PATCH 04/20] [Jenkins] auto-formatting by clang-format version
 5.0.0-3~16.04.1 (tags/RELEASE_500/final)

---
 stan/math/rev/core/autodiffstackstorage.hpp   | 17 +++++++-------
 stan/math/rev/core/chainablestack_inst.cpp    |  4 ++--
 .../math/rev/mat/functor/gradient_test.cpp    | 22 +++++++++----------
 3 files changed, 22 insertions(+), 21 deletions(-)

diff --git a/stan/math/rev/core/autodiffstackstorage.hpp b/stan/math/rev/core/autodiffstackstorage.hpp
index 58d8140830d..29259e5dff2 100644
--- a/stan/math/rev/core/autodiffstackstorage.hpp
+++ b/stan/math/rev/core/autodiffstackstorage.hpp
@@ -66,29 +66,30 @@ struct AutodiffStackSingleton {
     return *instance_;
   }
 
-  static inline AutodiffStackStorage* instantiate() {
+  static inline AutodiffStackStorage *instantiate() {
     static thread_local std::mutex init_mutex;
     std::lock_guard<std::mutex> init_lock(init_mutex);
-    if(instance_ == 0)
+    if (instance_ == 0)
       instance_ = new AutodiffStackStorage();
     return instance_;
   }
-  
+
  private:
   static
 #ifdef STAN_THREADS
-  thread_local
+      thread_local
 #endif
-  AutodiffStackStorage* instance_;
+      AutodiffStackStorage *instance_;
 };
 
 template <typename ChainableT, typename ChainableAllocT>
 #ifdef STAN_THREADS
 thread_local
 #endif
-typename AutodiffStackSingleton<ChainableT,
-                                ChainableAllocT>::AutodiffStackStorage*
-AutodiffStackSingleton<ChainableT, ChainableAllocT>::instance_ = 0;
+    typename AutodiffStackSingleton<ChainableT,
+                                    ChainableAllocT>::AutodiffStackStorage
+        *AutodiffStackSingleton<ChainableT, ChainableAllocT>::instance_
+    = 0;
 
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/core/chainablestack_inst.cpp b/stan/math/rev/core/chainablestack_inst.cpp
index 210a66cd8f6..2a230804d90 100644
--- a/stan/math/rev/core/chainablestack_inst.cpp
+++ b/stan/math/rev/core/chainablestack_inst.cpp
@@ -7,8 +7,8 @@ namespace stan {
 namespace math {
 namespace internal {
 
-  const ChainableStack::AutodiffStackStorage* __chainable_stack = ChainableStack::instantiate();
-
+const ChainableStack::AutodiffStackStorage* __chainable_stack
+    = ChainableStack::instantiate();
 }
 }  // namespace math
 }  // namespace stan
diff --git a/test/unit/math/rev/mat/functor/gradient_test.cpp b/test/unit/math/rev/mat/functor/gradient_test.cpp
index 28c50e63717..14d3fabea2b 100644
--- a/test/unit/math/rev/mat/functor/gradient_test.cpp
+++ b/test/unit/math/rev/mat/functor/gradient_test.cpp
@@ -46,17 +46,17 @@ TEST(AgradAutoDiff, gradient_threaded) {
   EXPECT_FLOAT_EQ(x_ref(0) * x_ref(0) + 3 * 2 * x_ref(1), grad_fx_ref(1));
 
   auto thread_job = [&](double x1, double x2) {
-                      stan::math::ChainableStack::instantiate();
-                      double fx;
-                      VectorXd x_local(2);
-                      x_local << x1, x2;
-                      VectorXd grad_fx;
-                      stan::math::gradient(fun1(), x_local, fx, grad_fx);
-                      VectorXd res(1 + grad_fx.size());
-                      res(0) = fx;
-                      res.tail(grad_fx.size()) = grad_fx;
-                      return res;
-                    };
+    stan::math::ChainableStack::instantiate();
+    double fx;
+    VectorXd x_local(2);
+    x_local << x1, x2;
+    VectorXd grad_fx;
+    stan::math::gradient(fun1(), x_local, fx, grad_fx);
+    VectorXd res(1 + grad_fx.size());
+    res(0) = fx;
+    res.tail(grad_fx.size()) = grad_fx;
+    return res;
+  };
 
   // schedule a bunch of jobs which all do the same
   std::vector<std::future<VectorXd>> ad_futures_ref;

From fef78b5e640f97bb28bc46d7463d39aadc64c97d Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 19 Jan 2019 17:35:23 +0100
Subject: [PATCH 05/20] address PR comments from Sean: introduce init method,
 remove mutex there, add tests for init function

---
 stan/math/rev/core/autodiffstackstorage.hpp   | 11 +---
 stan/math/rev/core/chainablestack.hpp         | 12 ++++
 stan/math/rev/core/chainablestack_inst.cpp    |  3 +-
 stan/math/rev/mat/functor/map_rect_reduce.hpp |  2 +-
 test/unit/math/rev/core/init_test.cpp         | 60 +++++++++++++++++++
 .../math/rev/mat/functor/gradient_test.cpp    |  2 +-
 6 files changed, 76 insertions(+), 14 deletions(-)
 create mode 100644 test/unit/math/rev/core/init_test.cpp

diff --git a/stan/math/rev/core/autodiffstackstorage.hpp b/stan/math/rev/core/autodiffstackstorage.hpp
index 29259e5dff2..570026ce4f0 100644
--- a/stan/math/rev/core/autodiffstackstorage.hpp
+++ b/stan/math/rev/core/autodiffstackstorage.hpp
@@ -66,15 +66,6 @@ struct AutodiffStackSingleton {
     return *instance_;
   }
 
-  static inline AutodiffStackStorage *instantiate() {
-    static thread_local std::mutex init_mutex;
-    std::lock_guard<std::mutex> init_lock(init_mutex);
-    if (instance_ == 0)
-      instance_ = new AutodiffStackStorage();
-    return instance_;
-  }
-
- private:
   static
 #ifdef STAN_THREADS
       thread_local
@@ -89,7 +80,7 @@ thread_local
     typename AutodiffStackSingleton<ChainableT,
                                     ChainableAllocT>::AutodiffStackStorage
         *AutodiffStackSingleton<ChainableT, ChainableAllocT>::instance_
-    = 0;
+    = nullptr;
 
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/core/chainablestack.hpp b/stan/math/rev/core/chainablestack.hpp
index fe88c03cf29..f27dfc5d9bd 100644
--- a/stan/math/rev/core/chainablestack.hpp
+++ b/stan/math/rev/core/chainablestack.hpp
@@ -11,6 +11,18 @@ class chainable_alloc;
 
 typedef AutodiffStackSingleton<vari, chainable_alloc> ChainableStack;
 
+/**
+ * Instantiates an instance of the ChainableStack if not already
+ * initialized. This function must be called before any autodiff
+ * variable get's instantiated within any thread which performs
+ * reverse mode autodiff operations.
+ */
+static inline ChainableStack::AutodiffStackStorage* init() {
+  if (ChainableStack::instance_ == nullptr)
+    ChainableStack::instance_ = new ChainableStack::AutodiffStackStorage();
+  return ChainableStack::instance_;
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/core/chainablestack_inst.cpp b/stan/math/rev/core/chainablestack_inst.cpp
index 2a230804d90..70982c5a1cd 100644
--- a/stan/math/rev/core/chainablestack_inst.cpp
+++ b/stan/math/rev/core/chainablestack_inst.cpp
@@ -7,8 +7,7 @@ namespace stan {
 namespace math {
 namespace internal {
 
-const ChainableStack::AutodiffStackStorage* __chainable_stack
-    = ChainableStack::instantiate();
+const ChainableStack::AutodiffStackStorage* __chainable_stack = init();
 }
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/mat/functor/map_rect_reduce.hpp b/stan/math/rev/mat/functor/map_rect_reduce.hpp
index 1b1ca9dce72..176d8979abf 100644
--- a/stan/math/rev/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/rev/mat/functor/map_rect_reduce.hpp
@@ -24,7 +24,7 @@ struct map_rect_reduce<F, var, var> {
     const size_type num_shared_params = shared_params.rows();
     const size_type num_job_specific_params = job_specific_params.rows();
     matrix_d out(1 + num_shared_params + num_job_specific_params, 0);
-    ChainableStack::instantiate();
+    init();
 
     try {
       start_nested();
diff --git a/test/unit/math/rev/core/init_test.cpp b/test/unit/math/rev/core/init_test.cpp
new file mode 100644
index 00000000000..07bff9ac9ba
--- /dev/null
+++ b/test/unit/math/rev/core/init_test.cpp
@@ -0,0 +1,60 @@
+#include <gtest/gtest.h>
+#include <stan/math/rev/core.hpp>
+
+#include <thread>
+
+TEST(init, thread_initialize) {
+  // the main thread must be initialized by the time this code is
+  // reached
+  EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
+
+  stan::math::ChainableStack::AutodiffStackStorage& main_ad_stack
+      = stan::math::ChainableStack::instance();
+
+  EXPECT_TRUE(&main_ad_stack == stan::math::init());
+
+#ifdef STAN_THREADS
+  auto thread_tester = [&]() -> void {
+    EXPECT_TRUE(stan::math::ChainableStack::instance_ == nullptr);
+    stan::math::init();
+    EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
+    EXPECT_TRUE(stan::math::ChainableStack::instance_ != &main_ad_stack);
+  };
+#else
+  auto thread_tester = [&]() -> void {
+    EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
+    stan::math::init();
+    EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
+    EXPECT_TRUE(stan::math::ChainableStack::instance_ != &main_ad_stack);
+  };
+#endif
+  std::thread other_work(thread_tester);
+
+  other_work.join();
+}
+
+TEST(init, thread_instances) {
+  // place a var on the stack
+  stan::math::var a = 1;
+
+  stan::math::ChainableStack::AutodiffStackStorage& main_ad_stack
+      = stan::math::ChainableStack::instance();
+
+#ifdef STAN_THREADS
+  auto thread_tester = [&]() -> void {
+    stan::math::init();
+    EXPECT_TRUE(main_ad_stack.var_stack_.size()
+                > stan::math::ChainableStack::instance().var_stack_.size());
+  };
+#else
+  auto thread_tester = [&]() -> void {
+    stan::math::init();
+    EXPECT_TRUE(main_ad_stack.var_stack_.size()
+                == stan::math::ChainableStack::instance().var_stack_.size());
+  };
+#endif
+
+  std::thread other_work(thread_tester);
+
+  other_work.join();
+}
diff --git a/test/unit/math/rev/mat/functor/gradient_test.cpp b/test/unit/math/rev/mat/functor/gradient_test.cpp
index 14d3fabea2b..d4a2deadef4 100644
--- a/test/unit/math/rev/mat/functor/gradient_test.cpp
+++ b/test/unit/math/rev/mat/functor/gradient_test.cpp
@@ -46,7 +46,7 @@ TEST(AgradAutoDiff, gradient_threaded) {
   EXPECT_FLOAT_EQ(x_ref(0) * x_ref(0) + 3 * 2 * x_ref(1), grad_fx_ref(1));
 
   auto thread_job = [&](double x1, double x2) {
-    stan::math::ChainableStack::instantiate();
+    stan::math::init();
     double fx;
     VectorXd x_local(2);
     x_local << x1, x2;

From 7a3fffe3677146abf015b773cd9f63937520d256 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 19 Jan 2019 17:40:24 +0100
Subject: [PATCH 06/20] some more comments on test & makefile alignment

---
 make/libraries                        | 2 +-
 test/unit/math/rev/core/init_test.cpp | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/make/libraries b/make/libraries
index 01f2dc32bf0..85f38999a24 100644
--- a/make/libraries
+++ b/make/libraries
@@ -110,7 +110,7 @@ clean-mpi:
 
 endif
 
-STACK_INSTANTIATION_CPP := $(shell find $(MATH)stan -type f -name 'chainablestack_inst.cpp')
+STACK_INSTANTIATION_CPP := $(MATH)stan/math/rev/core/chainablestack_inst.cpp
 STACK_INSTANTIATION := $(STACK_INSTANTIATION_CPP:%.cpp=%.o)
 
 $(STACK_INSTANTIATION) : CXXFLAGS += -fPIC
diff --git a/test/unit/math/rev/core/init_test.cpp b/test/unit/math/rev/core/init_test.cpp
index 07bff9ac9ba..2887f7cc2df 100644
--- a/test/unit/math/rev/core/init_test.cpp
+++ b/test/unit/math/rev/core/init_test.cpp
@@ -25,7 +25,7 @@ TEST(init, thread_initialize) {
     EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
     stan::math::init();
     EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
-    EXPECT_TRUE(stan::math::ChainableStack::instance_ != &main_ad_stack);
+    EXPECT_TRUE(stan::math::ChainableStack::instance_ == &main_ad_stack);
   };
 #endif
   std::thread other_work(thread_tester);
@@ -34,7 +34,9 @@ TEST(init, thread_initialize) {
 }
 
 TEST(init, thread_instances) {
-  // place a var on the stack
+  // place a var on the stack such that a fresh stack in another
+  // thread will be different at initialization (if STAN_THREADS is
+  // set)
   stan::math::var a = 1;
 
   stan::math::ChainableStack::AutodiffStackStorage& main_ad_stack

From 710d387aaf4f2b6e0a9420bc7e9014ea3205fe08 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 19 Jan 2019 18:02:15 +0100
Subject: [PATCH 07/20] update doc on AutodiffStackSingleton

---
 stan/math/rev/core/autodiffstackstorage.hpp | 35 ++++++++++++++++-----
 1 file changed, 28 insertions(+), 7 deletions(-)

diff --git a/stan/math/rev/core/autodiffstackstorage.hpp b/stan/math/rev/core/autodiffstackstorage.hpp
index 570026ce4f0..2d8544cc644 100644
--- a/stan/math/rev/core/autodiffstackstorage.hpp
+++ b/stan/math/rev/core/autodiffstackstorage.hpp
@@ -9,21 +9,41 @@ namespace stan {
 namespace math {
 
 /**
- * Provides a thread_local singleton if needed. Read warnings below!
- * For performance reasons the singleton is a global static for the
- * case of no threading which is returned by a function. This design
- * should allow the compiler to apply necessary inlining to get
+ * Provides a thread_local storage (TLS) singleton if needed. Read
+ * warnings below!  For performance reasons the singleton is a global
+ * static pointer which is initialized to a constant expression (the
+ * nullptr). Doing allows for a faster access to the TLS pointer
+ * instance. If one were to directly initialize the pointer with a
+ * call to new, the compiler will insert additional code for each
+ * reference in the code to the TLS. The additional code checks the
+ * initialization status thus everytime one accesses the TLS, see [4]
+ * for details and an example which demonstrates this. However, this
+ * design requires that the pointer must be initialized with the
+ * init() function before any var's are
+ * instantiated. Otherwise a segfault will occur. The init() function
+ * must be called for any thread which wants to use reverse mode
+ * autodiff facilites.
+ *
+ * Furthermore, the declaration as a global (possibly thread local)
+ * pointer allows the compiler to apply necessary inlining to get
  * maximal performance. However, this design suffers from "the static
  * init order fiasco"[0].  Anywhere this is used, we must be
  * absolutely positive that it doesn't matter when the singleton will
  * get initialized relative to other static variables.  In exchange,
- * we get a more performant singleton pattern for the non-threading
- * case. In the threading case we use the defacto standard C++11
+ * we get a more performant singleton pattern.
+ *
+ * Formely, stan-math used in the threading case the defacto standard C++11
  * singleton pattern relying on a function wrapping a static local
  * variable. This standard pattern is expected to be well supported
  * by the major compilers (as its standard), but it does incur some
  * performance penalty.  There has been some discussion on this; see
- * [1] and [2] and the discussions those PRs link to as well.
+ * [1] and [2] and the discussions those PRs link to as well. The
+ * current design has a much reduced/almost no performance penalty
+ * since the access to the TLS is not wrapped in extra function
+ * calls. Moreover, the thread_local declaration below can be changed
+ * to the __thread keyword which is a compiler-specific extension of
+ * g++,clang++&intel and is around for much longer than C++11 such
+ * that these compilers should support this design just as robust.
  *
  * These are thread_local only if the user asks for it with
  * -DSTAN_THREADS. This is primarily because Apple clang compilers
@@ -38,6 +58,7 @@ namespace math {
  * [2] https://github.com/stan-dev/math/pull/826
  * [3]
  * http://discourse.mc-stan.org/t/potentially-dropping-support-for-older-versions-of-apples-version-of-clang/3780/
+ * [4] https://discourse.mc-stan.org/t/thread-performance-penalty/7306/8
  */
 template <typename ChainableT, typename ChainableAllocT>
 struct AutodiffStackSingleton {

From a81d7549835355a4c3271380696d58ea113a5514 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 19 Jan 2019 20:59:32 +0100
Subject: [PATCH 08/20] add init to all map_rect_reduce flavors

---
 stan/math/prim/mat/functor/map_rect_reduce.hpp | 5 +++++
 stan/math/rev/mat/functor/map_rect_reduce.hpp  | 2 ++
 2 files changed, 7 insertions(+)

diff --git a/stan/math/prim/mat/functor/map_rect_reduce.hpp b/stan/math/prim/mat/functor/map_rect_reduce.hpp
index 68486fecd12..90419edc8c3 100644
--- a/stan/math/prim/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/prim/mat/functor/map_rect_reduce.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_PRIM_MAT_FUNCTOR_MAP_RECT_REDUCE_HPP
 
 #include <stan/math/prim/mat/fun/typedefs.hpp>
+#include <stan/math/rev/core/chainablestack.hpp>
 
 #include <vector>
 
@@ -43,6 +44,10 @@ class map_rect_reduce<F, double, double> {
                       const std::vector<double>& x_r,
                       const std::vector<int>& x_i,
                       std::ostream* msgs = nullptr) const {
+    // It is rather unclear why the init of the AD stack must happen
+    // here. Otherwise we get segfaults with example stan programs
+    // which cannot be reproduced (yet) using tests.
+    init();
     return F()(shared_params, job_specific_params, x_r, x_i, msgs).transpose();
   }
 };
diff --git a/stan/math/rev/mat/functor/map_rect_reduce.hpp b/stan/math/rev/mat/functor/map_rect_reduce.hpp
index 176d8979abf..642d4b32eea 100644
--- a/stan/math/rev/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/rev/mat/functor/map_rect_reduce.hpp
@@ -66,6 +66,7 @@ struct map_rect_reduce<F, double, var> {
                       std::ostream* msgs = nullptr) const {
     const size_type num_job_specific_params = job_specific_params.rows();
     matrix_d out(1 + num_job_specific_params, 0);
+    init();
 
     try {
       start_nested();
@@ -102,6 +103,7 @@ struct map_rect_reduce<F, var, double> {
                       std::ostream* msgs = nullptr) const {
     const size_type num_shared_params = shared_params.rows();
     matrix_d out(1 + num_shared_params, 0);
+    init();
 
     try {
       start_nested();

From a3ee346ed905464dea1a5f4c5de545e6ba443f92 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 19 Jan 2019 21:09:52 +0100
Subject: [PATCH 09/20] add more notes

---
 stan/math/prim/mat/functor/map_rect_reduce.hpp              | 4 +++-
 test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp | 2 +-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/stan/math/prim/mat/functor/map_rect_reduce.hpp b/stan/math/prim/mat/functor/map_rect_reduce.hpp
index 90419edc8c3..f469596c4a8 100644
--- a/stan/math/prim/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/prim/mat/functor/map_rect_reduce.hpp
@@ -46,7 +46,9 @@ class map_rect_reduce<F, double, double> {
                       std::ostream* msgs = nullptr) const {
     // It is rather unclear why the init of the AD stack must happen
     // here. Otherwise we get segfaults with example stan programs
-    // which cannot be reproduced (yet) using tests.
+    // which cannot be reproduced (yet) using tests. The unit tests
+    // for prim and rev work just fine when the init is not given
+    // here!
     init();
     return F()(shared_params, job_specific_params, x_r, x_i, msgs).transpose();
   }
diff --git a/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp b/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp
index 4aebcdd56cc..46e3bc77315 100644
--- a/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp
+++ b/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp
@@ -20,7 +20,7 @@ struct map_rect : public ::testing::Test {
   std::vector<Eigen::VectorXd> job_params_d;
   std::vector<std::vector<double> > x_r;
   std::vector<std::vector<int> > x_i;
-  const std::size_t N = 10;
+  const std::size_t N = 100;
 
   virtual void SetUp() {
     shared_params_d.resize(2);

From 63cc01d8e2201f0bcc85bae5e61beaf25f57eb32 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sun, 20 Jan 2019 11:53:01 +0100
Subject: [PATCH 10/20] move all map_rect_reduce specialisations into rev in
 order to always initialize the AD tape

---
 .../math/prim/mat/functor/map_rect_reduce.hpp | 26 +++++--------------
 stan/math/rev/mat/functor/map_rect_reduce.hpp | 13 ++++++++++
 2 files changed, 20 insertions(+), 19 deletions(-)

diff --git a/stan/math/prim/mat/functor/map_rect_reduce.hpp b/stan/math/prim/mat/functor/map_rect_reduce.hpp
index f469596c4a8..fcb10ab1dec 100644
--- a/stan/math/prim/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/prim/mat/functor/map_rect_reduce.hpp
@@ -2,7 +2,6 @@
 #define STAN_MATH_PRIM_MAT_FUNCTOR_MAP_RECT_REDUCE_HPP
 
 #include <stan/math/prim/mat/fun/typedefs.hpp>
-#include <stan/math/rev/core/chainablestack.hpp>
 
 #include <vector>
 
@@ -29,6 +28,13 @@ namespace internal {
  *
  * No higher order output format is defined yet.
  *
+ * In the prim folder we only keep the basic defintion while all
+ * specialisations are in rev. This is to ensure that the AD tape of
+ * all thread can be initialized. This includes the variant which has
+ * no var arguments involved since nested AD operations can be
+ * performed (as for example is the case when calling the CVODES
+ * integrator which needs the Jacobian).
+ *
  * @tparam F user functor
  * @tparam T_shared_param type of shared parameters
  * @tparam T_job_param type of job specific parameters
@@ -36,24 +42,6 @@ namespace internal {
 template <typename F, typename T_shared_param, typename T_job_param>
 class map_rect_reduce {};
 
-template <typename F>
-class map_rect_reduce<F, double, double> {
- public:
-  matrix_d operator()(const vector_d& shared_params,
-                      const vector_d& job_specific_params,
-                      const std::vector<double>& x_r,
-                      const std::vector<int>& x_i,
-                      std::ostream* msgs = nullptr) const {
-    // It is rather unclear why the init of the AD stack must happen
-    // here. Otherwise we get segfaults with example stan programs
-    // which cannot be reproduced (yet) using tests. The unit tests
-    // for prim and rev work just fine when the init is not given
-    // here!
-    init();
-    return F()(shared_params, job_specific_params, x_r, x_i, msgs).transpose();
-  }
-};
-
 }  // namespace internal
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/mat/functor/map_rect_reduce.hpp b/stan/math/rev/mat/functor/map_rect_reduce.hpp
index 642d4b32eea..581b5884031 100644
--- a/stan/math/rev/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/rev/mat/functor/map_rect_reduce.hpp
@@ -14,6 +14,19 @@ namespace stan {
 namespace math {
 namespace internal {
 
+template <typename F>
+class map_rect_reduce<F, double, double> {
+ public:
+  matrix_d operator()(const vector_d& shared_params,
+                      const vector_d& job_specific_params,
+                      const std::vector<double>& x_r,
+                      const std::vector<int>& x_i,
+                      std::ostream* msgs = nullptr) const {
+    init();
+    return F()(shared_params, job_specific_params, x_r, x_i, msgs).transpose();
+  }
+};
+
 template <typename F>
 struct map_rect_reduce<F, var, var> {
   matrix_d operator()(const vector_d& shared_params,

From 2757f804e4ceca935587c9a8454f6c4e9a1e8e6e Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sun, 20 Jan 2019 12:05:24 +0100
Subject: [PATCH 11/20] ensure that 4 threads are used during
 map_rect_concurrent tests

---
 .../mat/functor/map_rect_concurrent_test.cpp  |  7 ++++--
 .../map_rect_concurrent_threads_test.cpp      | 13 ++--------
 .../prim/mat/functor/num_threads_test.cpp     | 25 +++++++------------
 .../math/prim/mat/functor/utils_threads.hpp   | 23 +++++++++++++++++
 .../mat/functor/map_rect_concurrent_test.cpp  |  4 ++-
 5 files changed, 42 insertions(+), 30 deletions(-)
 create mode 100644 test/unit/math/prim/mat/functor/utils_threads.hpp

diff --git a/test/unit/math/prim/mat/functor/map_rect_concurrent_test.cpp b/test/unit/math/prim/mat/functor/map_rect_concurrent_test.cpp
index 23e9e902636..bbb3fb94303 100644
--- a/test/unit/math/prim/mat/functor/map_rect_concurrent_test.cpp
+++ b/test/unit/math/prim/mat/functor/map_rect_concurrent_test.cpp
@@ -5,10 +5,12 @@
 #undef STAN_MPI
 #endif
 
-#include <stan/math/prim/mat.hpp>
 #include <gtest/gtest.h>
+#include <stan/math/prim/mat.hpp>
+#include <stan/math/rev/mat.hpp>
 
 #include <test/unit/math/prim/mat/functor/hard_work.hpp>
+#include <test/unit/math/prim/mat/functor/utils_threads.hpp>
 
 #include <iostream>
 #include <vector>
@@ -21,9 +23,10 @@ struct map_rect : public ::testing::Test {
   std::vector<Eigen::VectorXd> job_params_d;
   std::vector<std::vector<double> > x_r;
   std::vector<std::vector<int> > x_i;
-  const int N = 10;
+  const int N = 100;
 
   virtual void SetUp() {
+    set_n_threads(4);
     shared_params_d.resize(2);
     shared_params_d << 2, 0;
 
diff --git a/test/unit/math/prim/mat/functor/map_rect_concurrent_threads_test.cpp b/test/unit/math/prim/mat/functor/map_rect_concurrent_threads_test.cpp
index d474c5d0aa3..24a8bf668c1 100644
--- a/test/unit/math/prim/mat/functor/map_rect_concurrent_threads_test.cpp
+++ b/test/unit/math/prim/mat/functor/map_rect_concurrent_threads_test.cpp
@@ -8,26 +8,17 @@
 #undef STAN_MPI
 #endif
 
-#include <stdlib.h>
-
 #include <gtest/gtest.h>
 #include <stan/math/prim/mat.hpp>
+#include <stan/math/rev/mat.hpp>
 
 #include <test/unit/math/prim/mat/functor/hard_work.hpp>
+#include <test/unit/math/prim/mat/functor/utils_threads.hpp>
 
 #include <iostream>
 #include <vector>
 #include <string>
 
-// utility to set number of threads to use
-void set_n_threads(int num_threads) {
-  static char env_string[256];
-  std::string num_threads_str = std::to_string(num_threads);
-  snprintf(env_string, sizeof(env_string), "STAN_NUM_THREADS=%s",
-           num_threads_str.c_str());
-  putenv(env_string);
-}
-
 STAN_REGISTER_MAP_RECT(0, hard_work)
 STAN_REGISTER_MAP_RECT(1, hard_work)
 
diff --git a/test/unit/math/prim/mat/functor/num_threads_test.cpp b/test/unit/math/prim/mat/functor/num_threads_test.cpp
index 87023b4b42a..35aa29199ae 100644
--- a/test/unit/math/prim/mat/functor/num_threads_test.cpp
+++ b/test/unit/math/prim/mat/functor/num_threads_test.cpp
@@ -3,43 +3,36 @@
 #endif
 
 #include <stan/math/prim/mat/functor/map_rect_concurrent.hpp>
+
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
-
-#include <stdlib.h>
-
-// Can't easily use std::string as putenv require non-const char*
-void set_n_threads_var(const char* value) {
-  static char env_string[256];
-  snprintf(env_string, sizeof(env_string), "STAN_NUM_THREADS=%s", value);
-  putenv(env_string);
-}
+#include <test/unit/math/prim/mat/functor/utils_threads.hpp>
 
 TEST(num_threads, correct_values) {
-  set_n_threads_var("10");
+  set_n_threads("10");
   EXPECT_EQ(stan::math::internal::get_num_threads(100), 10);
 
-  set_n_threads_var("4");
+  set_n_threads("4");
   EXPECT_EQ(stan::math::internal::get_num_threads(3), 3);
 
-  set_n_threads_var("-1");
+  set_n_threads("-1");
   EXPECT_GE(stan::math::internal::get_num_threads(5), 1);
 }
 
 TEST(num_threads, incorrect_values) {
-  set_n_threads_var("abc");
+  set_n_threads("abc");
   EXPECT_THROW_MSG(stan::math::internal::get_num_threads(5),
                    std::invalid_argument, "positive number or -1");
 
-  set_n_threads_var("1c");
+  set_n_threads("1c");
   EXPECT_THROW_MSG(stan::math::internal::get_num_threads(5),
                    std::invalid_argument, "positive number or -1");
 
-  set_n_threads_var("-2");
+  set_n_threads("-2");
   EXPECT_THROW_MSG(stan::math::internal::get_num_threads(5),
                    std::invalid_argument, "must be positive or -1");
 
-  set_n_threads_var("0");
+  set_n_threads("0");
   EXPECT_THROW_MSG(stan::math::internal::get_num_threads(5),
                    std::invalid_argument, "must be positive or -1");
 }
diff --git a/test/unit/math/prim/mat/functor/utils_threads.hpp b/test/unit/math/prim/mat/functor/utils_threads.hpp
new file mode 100644
index 00000000000..9d92cf3a767
--- /dev/null
+++ b/test/unit/math/prim/mat/functor/utils_threads.hpp
@@ -0,0 +1,23 @@
+#ifndef TEST_UNIT_MATH_PRIM_MAT_FUNCTOR_UTILS_THREADS_HPP
+#define TEST_UNIT_MATH_PRIM_MAT_FUNCTOR_UTILS_THREADS_HPP
+
+#include <stdlib.h>
+#include <string>
+
+// utility to set number of threads to use
+void set_n_threads(int num_threads) {
+  static char env_string[256];
+  std::string num_threads_str = std::to_string(num_threads);
+  snprintf(env_string, sizeof(env_string), "STAN_NUM_THREADS=%s",
+           num_threads_str.c_str());
+  putenv(env_string);
+}
+
+// Can't easily use std::string as putenv require non-const char*
+void set_n_threads(const char* value) {
+  static char env_string[256];
+  snprintf(env_string, sizeof(env_string), "STAN_NUM_THREADS=%s", value);
+  putenv(env_string);
+}
+
+#endif
diff --git a/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp b/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp
index 46e3bc77315..f39ff8d5e0d 100644
--- a/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp
+++ b/test/unit/math/rev/mat/functor/map_rect_concurrent_test.cpp
@@ -5,10 +5,11 @@
 #undef STAN_MPI
 #endif
 
-#include <stan/math/rev/mat.hpp>
 #include <gtest/gtest.h>
+#include <stan/math/rev/mat.hpp>
 
 #include <test/unit/math/prim/mat/functor/hard_work.hpp>
+#include <test/unit/math/prim/mat/functor/utils_threads.hpp>
 
 #include <iostream>
 #include <vector>
@@ -23,6 +24,7 @@ struct map_rect : public ::testing::Test {
   const std::size_t N = 100;
 
   virtual void SetUp() {
+    set_n_threads(4);
     shared_params_d.resize(2);
     shared_params_d << 2, 0;
 

From 8b271252d332455e5887cda6cfb84d869fbc0bb0 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Mon, 21 Jan 2019 21:33:30 +0100
Subject: [PATCH 12/20] move map_rect_concurrent definition to rev branch and
 place init() call in lambda

---
 .../prim/mat/functor/map_rect_concurrent.hpp  |  68 +-----------
 .../math/prim/mat/functor/map_rect_reduce.hpp |  12 +++
 stan/math/rev/mat.hpp                         |   1 +
 .../rev/mat/functor/map_rect_concurrent.hpp   | 101 ++++++++++++++++++
 stan/math/rev/mat/functor/map_rect_reduce.hpp |  17 ---
 .../math/prim/mat/functor/map_rect_test.cpp   |   1 +
 6 files changed, 116 insertions(+), 84 deletions(-)
 create mode 100644 stan/math/rev/mat/functor/map_rect_concurrent.hpp

diff --git a/stan/math/prim/mat/functor/map_rect_concurrent.hpp b/stan/math/prim/mat/functor/map_rect_concurrent.hpp
index 0a739691cc8..24b65c64bb4 100644
--- a/stan/math/prim/mat/functor/map_rect_concurrent.hpp
+++ b/stan/math/prim/mat/functor/map_rect_concurrent.hpp
@@ -10,7 +10,6 @@
 
 #include <vector>
 #include <thread>
-#include <future>
 #include <cstdlib>
 
 namespace stan {
@@ -77,72 +76,7 @@ map_rect_concurrent(
     const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
         job_params,
     const std::vector<std::vector<double>>& x_r,
-    const std::vector<std::vector<int>>& x_i, std::ostream* msgs = nullptr) {
-  typedef map_rect_reduce<F, T_shared_param, T_job_param> ReduceF;
-  typedef map_rect_combine<F, T_shared_param, T_job_param> CombineF;
-
-  const int num_jobs = job_params.size();
-  const vector_d shared_params_dbl = value_of(shared_params);
-  std::vector<std::future<std::vector<matrix_d>>> futures;
-
-  auto execute_chunk = [&](int start, int size) -> std::vector<matrix_d> {
-    const int end = start + size;
-    std::vector<matrix_d> chunk_f_out;
-    chunk_f_out.reserve(size);
-    for (int i = start; i != end; i++)
-      chunk_f_out.push_back(ReduceF()(
-          shared_params_dbl, value_of(job_params[i]), x_r[i], x_i[i], msgs));
-    return chunk_f_out;
-  };
-
-  int num_threads = get_num_threads(num_jobs);
-  int num_jobs_per_thread = num_jobs / num_threads;
-  futures.emplace_back(
-      std::async(std::launch::deferred, execute_chunk, 0, num_jobs_per_thread));
-
-#ifdef STAN_THREADS
-  if (num_threads > 1) {
-    const int num_big_threads = num_jobs % num_threads;
-    const int first_big_thread = num_threads - num_big_threads;
-    for (int i = 1, job_start = num_jobs_per_thread, job_size = 0;
-         i < num_threads; ++i, job_start += job_size) {
-      job_size = i >= first_big_thread ? num_jobs_per_thread + 1
-                                       : num_jobs_per_thread;
-      futures.emplace_back(
-          std::async(std::launch::async, execute_chunk, job_start, job_size));
-    }
-  }
-#endif
-
-  // collect results
-  std::vector<int> world_f_out;
-  world_f_out.reserve(num_jobs);
-  matrix_d world_output(0, 0);
-
-  int offset = 0;
-  for (std::size_t i = 0; i < futures.size(); ++i) {
-    const std::vector<matrix_d>& chunk_result = futures[i].get();
-    if (i == 0)
-      world_output.resize(chunk_result[0].rows(),
-                          num_jobs * chunk_result[0].cols());
-
-    for (const auto& job_result : chunk_result) {
-      const int num_job_outputs = job_result.cols();
-      world_f_out.push_back(num_job_outputs);
-
-      if (world_output.cols() < offset + num_job_outputs)
-        world_output.conservativeResize(Eigen::NoChange,
-                                        2 * (offset + num_job_outputs));
-
-      world_output.block(0, offset, world_output.rows(), num_job_outputs)
-          = job_result;
-
-      offset += num_job_outputs;
-    }
-  }
-  CombineF combine(shared_params, job_params);
-  return combine(world_output, world_f_out);
-}
+    const std::vector<std::vector<int>>& x_i, std::ostream* msgs = nullptr);
 
 }  // namespace internal
 }  // namespace math
diff --git a/stan/math/prim/mat/functor/map_rect_reduce.hpp b/stan/math/prim/mat/functor/map_rect_reduce.hpp
index fcb10ab1dec..256241c3b16 100644
--- a/stan/math/prim/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/prim/mat/functor/map_rect_reduce.hpp
@@ -42,6 +42,18 @@ namespace internal {
 template <typename F, typename T_shared_param, typename T_job_param>
 class map_rect_reduce {};
 
+template <typename F>
+class map_rect_reduce<F, double, double> {
+ public:
+  matrix_d operator()(const vector_d& shared_params,
+                      const vector_d& job_specific_params,
+                      const std::vector<double>& x_r,
+                      const std::vector<int>& x_i,
+                      std::ostream* msgs = nullptr) const {
+    return F()(shared_params, job_specific_params, x_r, x_i, msgs).transpose();
+  }
+};
+
 }  // namespace internal
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/mat.hpp b/stan/math/rev/mat.hpp
index 9138a7b9827..5b39c0cf35c 100644
--- a/stan/math/rev/mat.hpp
+++ b/stan/math/rev/mat.hpp
@@ -68,6 +68,7 @@
 #include <stan/math/rev/mat/functor/integrate_ode_adams.hpp>
 #include <stan/math/rev/mat/functor/integrate_ode_bdf.hpp>
 #include <stan/math/rev/mat/functor/integrate_dae.hpp>
+#include <stan/math/rev/mat/functor/map_rect_concurrent.hpp>
 #include <stan/math/rev/mat/functor/map_rect_reduce.hpp>
 
 #endif
diff --git a/stan/math/rev/mat/functor/map_rect_concurrent.hpp b/stan/math/rev/mat/functor/map_rect_concurrent.hpp
new file mode 100644
index 00000000000..2e4e405523a
--- /dev/null
+++ b/stan/math/rev/mat/functor/map_rect_concurrent.hpp
@@ -0,0 +1,101 @@
+#ifndef STAN_MATH_REV_MAT_FUNCTOR_MAP_RECT_CONCURRENT_HPP
+#define STAN_MATH_REV_MAT_FUNCTOR_MAP_RECT_CONCURRENT_HPP
+
+#include <stan/math/prim/mat/fun/typedefs.hpp>
+
+#include <stan/math/prim/mat/functor/map_rect_reduce.hpp>
+#include <stan/math/prim/mat/functor/map_rect_concurrent.hpp>
+#include <stan/math/prim/scal/err/invalid_argument.hpp>
+#include <stan/math/rev/mat/functor/map_rect_reduce.hpp>
+#include <stan/math/rev/core/chainablestack.hpp>
+
+#include <vector>
+#include <thread>
+#include <future>
+
+namespace stan {
+namespace math {
+namespace internal {
+
+template <int call_id, typename F, typename T_shared_param,
+          typename T_job_param>
+Eigen::Matrix<typename stan::return_type<T_shared_param, T_job_param>::type,
+              Eigen::Dynamic, 1>
+map_rect_concurrent(
+    const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
+    const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
+        job_params,
+    const std::vector<std::vector<double>>& x_r,
+    const std::vector<std::vector<int>>& x_i, std::ostream* msgs) {
+  typedef map_rect_reduce<F, T_shared_param, T_job_param> ReduceF;
+  typedef map_rect_combine<F, T_shared_param, T_job_param> CombineF;
+
+  const int num_jobs = job_params.size();
+  const vector_d shared_params_dbl = value_of(shared_params);
+  std::vector<std::future<std::vector<matrix_d>>> futures;
+
+  auto execute_chunk = [&](int start, int size) -> std::vector<matrix_d> {
+    const int end = start + size;
+    init();
+    std::vector<matrix_d> chunk_f_out;
+    chunk_f_out.reserve(size);
+    for (int i = start; i != end; i++)
+      chunk_f_out.push_back(ReduceF()(
+          shared_params_dbl, value_of(job_params[i]), x_r[i], x_i[i], msgs));
+    return chunk_f_out;
+  };
+
+  int num_threads = get_num_threads(num_jobs);
+  int num_jobs_per_thread = num_jobs / num_threads;
+  futures.emplace_back(
+      std::async(std::launch::deferred, execute_chunk, 0, num_jobs_per_thread));
+
+#ifdef STAN_THREADS
+  if (num_threads > 1) {
+    const int num_big_threads = num_jobs % num_threads;
+    const int first_big_thread = num_threads - num_big_threads;
+    for (int i = 1, job_start = num_jobs_per_thread, job_size = 0;
+         i < num_threads; ++i, job_start += job_size) {
+      job_size = i >= first_big_thread ? num_jobs_per_thread + 1
+                                       : num_jobs_per_thread;
+      futures.emplace_back(
+          std::async(std::launch::async, execute_chunk, job_start, job_size));
+    }
+  }
+#endif
+
+  // collect results
+  std::vector<int> world_f_out;
+  world_f_out.reserve(num_jobs);
+  matrix_d world_output(0, 0);
+
+  int offset = 0;
+  for (std::size_t i = 0; i < futures.size(); ++i) {
+    const std::vector<matrix_d>& chunk_result = futures[i].get();
+    if (i == 0)
+      world_output.resize(chunk_result[0].rows(),
+                          num_jobs * chunk_result[0].cols());
+
+    for (const auto& job_result : chunk_result) {
+      const int num_job_outputs = job_result.cols();
+      world_f_out.push_back(num_job_outputs);
+
+      if (world_output.cols() < offset + num_job_outputs)
+        world_output.conservativeResize(Eigen::NoChange,
+                                        2 * (offset + num_job_outputs));
+
+      world_output.block(0, offset, world_output.rows(), num_job_outputs)
+          = job_result;
+
+      offset += num_job_outputs;
+    }
+  }
+  CombineF combine(shared_params, job_params);
+  return combine(world_output, world_f_out);
+}
+
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/mat/functor/map_rect_reduce.hpp b/stan/math/rev/mat/functor/map_rect_reduce.hpp
index 581b5884031..4d9fe49d030 100644
--- a/stan/math/rev/mat/functor/map_rect_reduce.hpp
+++ b/stan/math/rev/mat/functor/map_rect_reduce.hpp
@@ -3,7 +3,6 @@
 
 #include <stan/math/prim/mat/functor/map_rect_reduce.hpp>
 #include <stan/math/prim/mat/fun/typedefs.hpp>
-#include <stan/math/rev/core/chainablestack.hpp>
 #include <stan/math/rev/core/var.hpp>
 #include <stan/math/rev/mat/fun/typedefs.hpp>
 #include <stan/math/rev/mat/fun/to_var.hpp>
@@ -14,19 +13,6 @@ namespace stan {
 namespace math {
 namespace internal {
 
-template <typename F>
-class map_rect_reduce<F, double, double> {
- public:
-  matrix_d operator()(const vector_d& shared_params,
-                      const vector_d& job_specific_params,
-                      const std::vector<double>& x_r,
-                      const std::vector<int>& x_i,
-                      std::ostream* msgs = nullptr) const {
-    init();
-    return F()(shared_params, job_specific_params, x_r, x_i, msgs).transpose();
-  }
-};
-
 template <typename F>
 struct map_rect_reduce<F, var, var> {
   matrix_d operator()(const vector_d& shared_params,
@@ -37,7 +23,6 @@ struct map_rect_reduce<F, var, var> {
     const size_type num_shared_params = shared_params.rows();
     const size_type num_job_specific_params = job_specific_params.rows();
     matrix_d out(1 + num_shared_params + num_job_specific_params, 0);
-    init();
 
     try {
       start_nested();
@@ -79,7 +64,6 @@ struct map_rect_reduce<F, double, var> {
                       std::ostream* msgs = nullptr) const {
     const size_type num_job_specific_params = job_specific_params.rows();
     matrix_d out(1 + num_job_specific_params, 0);
-    init();
 
     try {
       start_nested();
@@ -116,7 +100,6 @@ struct map_rect_reduce<F, var, double> {
                       std::ostream* msgs = nullptr) const {
     const size_type num_shared_params = shared_params.rows();
     matrix_d out(1 + num_shared_params, 0);
-    init();
 
     try {
       start_nested();
diff --git a/test/unit/math/prim/mat/functor/map_rect_test.cpp b/test/unit/math/prim/mat/functor/map_rect_test.cpp
index 49696b0c058..93497f07c42 100644
--- a/test/unit/math/prim/mat/functor/map_rect_test.cpp
+++ b/test/unit/math/prim/mat/functor/map_rect_test.cpp
@@ -1,4 +1,5 @@
 #include <stan/math/prim/mat.hpp>
+#include <stan/math/rev/mat.hpp>
 #include <gtest/gtest.h>
 
 #include <test/unit/math/prim/mat/functor/hard_work.hpp>

From 4a521b9fa4accc40f414b6d62633405713b44b3a Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Mon, 21 Jan 2019 22:08:26 +0100
Subject: [PATCH 13/20] add missing include

---
 test/unit/math/prim/mat/functor/map_rect_mpi_test.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/unit/math/prim/mat/functor/map_rect_mpi_test.cpp b/test/unit/math/prim/mat/functor/map_rect_mpi_test.cpp
index d7d2ebc4e47..0f2e7aee549 100644
--- a/test/unit/math/prim/mat/functor/map_rect_mpi_test.cpp
+++ b/test/unit/math/prim/mat/functor/map_rect_mpi_test.cpp
@@ -3,6 +3,7 @@
 #ifdef STAN_MPI
 
 #include <stan/math/prim/mat.hpp>
+#include <stan/math/rev/mat.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
 

From 223df9559be02717231c22c4869a18afc30740d6 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sebastian.weber@novartis.com>
Date: Wed, 23 Jan 2019 09:37:37 +0100
Subject: [PATCH 14/20] address PR review comments

---
 stan/math/rev/core/autodiffstackstorage.hpp | 1 -
 stan/math/rev/core/chainablestack_inst.cpp  | 5 +++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/stan/math/rev/core/autodiffstackstorage.hpp b/stan/math/rev/core/autodiffstackstorage.hpp
index 2d8544cc644..173fb70511f 100644
--- a/stan/math/rev/core/autodiffstackstorage.hpp
+++ b/stan/math/rev/core/autodiffstackstorage.hpp
@@ -2,7 +2,6 @@
 #define STAN_MATH_REV_CORE_AUTODIFFSTACKSTORAGE_HPP
 
 #include <stan/math/memory/stack_alloc.hpp>
-#include <mutex>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/rev/core/chainablestack_inst.cpp b/stan/math/rev/core/chainablestack_inst.cpp
index 70982c5a1cd..041dc52029e 100644
--- a/stan/math/rev/core/chainablestack_inst.cpp
+++ b/stan/math/rev/core/chainablestack_inst.cpp
@@ -6,8 +6,9 @@
 namespace stan {
 namespace math {
 namespace internal {
-
-const ChainableStack::AutodiffStackStorage* __chainable_stack = init();
+  namespace {
+    const ChainableStack::AutodiffStackStorage* __chainable_stack = init();
+  }
 }
 }  // namespace math
 }  // namespace stan

From 61e937ccde595ea62ef59fe3d7b73ddd1c834cd0 Mon Sep 17 00:00:00 2001
From: Stan Jenkins <mc.stanislaw@gmail.com>
Date: Wed, 23 Jan 2019 08:40:33 +0000
Subject: [PATCH 15/20] [Jenkins] auto-formatting by clang-format version
 5.0.0-3~16.04.1 (tags/RELEASE_500/final)

---
 stan/math/rev/core/chainablestack_inst.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/stan/math/rev/core/chainablestack_inst.cpp b/stan/math/rev/core/chainablestack_inst.cpp
index 041dc52029e..e9a68e2eed6 100644
--- a/stan/math/rev/core/chainablestack_inst.cpp
+++ b/stan/math/rev/core/chainablestack_inst.cpp
@@ -6,10 +6,10 @@
 namespace stan {
 namespace math {
 namespace internal {
-  namespace {
-    const ChainableStack::AutodiffStackStorage* __chainable_stack = init();
-  }
+namespace {
+const ChainableStack::AutodiffStackStorage* __chainable_stack = init();
 }
+}  // namespace internal
 }  // namespace math
 }  // namespace stan
 #endif

From fed0097f4ce4e1789371341a329047ce74006bd9 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 26 Jan 2019 17:23:30 +0100
Subject: [PATCH 16/20] make AD tape initialized without an init call for
 non-threaded case

---
 stan/math/rev/core/autodiffstackstorage.hpp       | 10 ++++++++++
 stan/math/rev/core/chainablestack.hpp             |  5 +++++
 stan/math/rev/core/chainablestack_inst.cpp        |  7 ++++---
 stan/math/rev/mat/functor/map_rect_concurrent.hpp |  2 +-
 test/unit/math/rev/core/init_test.cpp             | 10 +++++-----
 5 files changed, 25 insertions(+), 9 deletions(-)

diff --git a/stan/math/rev/core/autodiffstackstorage.hpp b/stan/math/rev/core/autodiffstackstorage.hpp
index 173fb70511f..40565a174b4 100644
--- a/stan/math/rev/core/autodiffstackstorage.hpp
+++ b/stan/math/rev/core/autodiffstackstorage.hpp
@@ -86,6 +86,12 @@ struct AutodiffStackSingleton {
     return *instance_;
   }
 
+  static AutodiffStackStorage *init() {
+    if (instance_ == nullptr)
+      instance_ = new AutodiffStackStorage();
+    return instance_;
+  }
+
   static
 #ifdef STAN_THREADS
       thread_local
@@ -100,7 +106,11 @@ thread_local
     typename AutodiffStackSingleton<ChainableT,
                                     ChainableAllocT>::AutodiffStackStorage
         *AutodiffStackSingleton<ChainableT, ChainableAllocT>::instance_
+#ifdef STAN_THREADS
     = nullptr;
+#else
+    = AutodiffStackSingleton<ChainableT, ChainableAllocT>::init();
+#endif
 
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/core/chainablestack.hpp b/stan/math/rev/core/chainablestack.hpp
index f27dfc5d9bd..44983b1ca65 100644
--- a/stan/math/rev/core/chainablestack.hpp
+++ b/stan/math/rev/core/chainablestack.hpp
@@ -17,11 +17,16 @@ typedef AutodiffStackSingleton<vari, chainable_alloc> ChainableStack;
  * variable get's instantiated within any thread which performs
  * reverse mode autodiff operations.
  */
+/*
 static inline ChainableStack::AutodiffStackStorage* init() {
+#ifdef STAN_THREADS
+  return ChainableStack::init();
   if (ChainableStack::instance_ == nullptr)
     ChainableStack::instance_ = new ChainableStack::AutodiffStackStorage();
+#endif
   return ChainableStack::instance_;
 }
+*/
 
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/core/chainablestack_inst.cpp b/stan/math/rev/core/chainablestack_inst.cpp
index e9a68e2eed6..24b558a922c 100644
--- a/stan/math/rev/core/chainablestack_inst.cpp
+++ b/stan/math/rev/core/chainablestack_inst.cpp
@@ -6,9 +6,10 @@
 namespace stan {
 namespace math {
 namespace internal {
-namespace {
-const ChainableStack::AutodiffStackStorage* __chainable_stack = init();
-}
+
+const ChainableStack::AutodiffStackStorage* __chainable_stack
+    = ChainableStack::init();
+
 }  // namespace internal
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/rev/mat/functor/map_rect_concurrent.hpp b/stan/math/rev/mat/functor/map_rect_concurrent.hpp
index 2e4e405523a..2bd2e32ab5b 100644
--- a/stan/math/rev/mat/functor/map_rect_concurrent.hpp
+++ b/stan/math/rev/mat/functor/map_rect_concurrent.hpp
@@ -36,7 +36,7 @@ map_rect_concurrent(
 
   auto execute_chunk = [&](int start, int size) -> std::vector<matrix_d> {
     const int end = start + size;
-    init();
+    ChainableStack::init();
     std::vector<matrix_d> chunk_f_out;
     chunk_f_out.reserve(size);
     for (int i = start; i != end; i++)
diff --git a/test/unit/math/rev/core/init_test.cpp b/test/unit/math/rev/core/init_test.cpp
index 2887f7cc2df..c47ca0c060d 100644
--- a/test/unit/math/rev/core/init_test.cpp
+++ b/test/unit/math/rev/core/init_test.cpp
@@ -11,19 +11,19 @@ TEST(init, thread_initialize) {
   stan::math::ChainableStack::AutodiffStackStorage& main_ad_stack
       = stan::math::ChainableStack::instance();
 
-  EXPECT_TRUE(&main_ad_stack == stan::math::init());
+  EXPECT_TRUE(&main_ad_stack == stan::math::ChainableStack::init());
 
 #ifdef STAN_THREADS
   auto thread_tester = [&]() -> void {
     EXPECT_TRUE(stan::math::ChainableStack::instance_ == nullptr);
-    stan::math::init();
+    stan::math::ChainableStack::init();
     EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
     EXPECT_TRUE(stan::math::ChainableStack::instance_ != &main_ad_stack);
   };
 #else
   auto thread_tester = [&]() -> void {
     EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
-    stan::math::init();
+    stan::math::ChainableStack::init();
     EXPECT_TRUE(stan::math::ChainableStack::instance_ != nullptr);
     EXPECT_TRUE(stan::math::ChainableStack::instance_ == &main_ad_stack);
   };
@@ -44,13 +44,13 @@ TEST(init, thread_instances) {
 
 #ifdef STAN_THREADS
   auto thread_tester = [&]() -> void {
-    stan::math::init();
+    stan::math::ChainableStack::init();
     EXPECT_TRUE(main_ad_stack.var_stack_.size()
                 > stan::math::ChainableStack::instance().var_stack_.size());
   };
 #else
   auto thread_tester = [&]() -> void {
-    stan::math::init();
+    stan::math::ChainableStack::init();
     EXPECT_TRUE(main_ad_stack.var_stack_.size()
                 == stan::math::ChainableStack::instance().var_stack_.size());
   };

From 9e15996764989e07b8152e8ad500b8e322b54720 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 26 Jan 2019 17:38:10 +0100
Subject: [PATCH 17/20] align test with new init scheme

---
 test/unit/math/rev/mat/functor/gradient_test.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/unit/math/rev/mat/functor/gradient_test.cpp b/test/unit/math/rev/mat/functor/gradient_test.cpp
index d4a2deadef4..2e4a16c1d43 100644
--- a/test/unit/math/rev/mat/functor/gradient_test.cpp
+++ b/test/unit/math/rev/mat/functor/gradient_test.cpp
@@ -46,7 +46,7 @@ TEST(AgradAutoDiff, gradient_threaded) {
   EXPECT_FLOAT_EQ(x_ref(0) * x_ref(0) + 3 * 2 * x_ref(1), grad_fx_ref(1));
 
   auto thread_job = [&](double x1, double x2) {
-    stan::math::init();
+    stan::math::ChainableStack::init();
     double fx;
     VectorXd x_local(2);
     x_local << x1, x2;

From 52fffe95f52276d8a3a5090a889be8aae8d07727 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Fri, 15 Feb 2019 20:47:36 +0100
Subject: [PATCH 18/20] initialize AD tape within unnamed namespace

---
 stan/math/rev/core/chainablestack.hpp | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/stan/math/rev/core/chainablestack.hpp b/stan/math/rev/core/chainablestack.hpp
index 44983b1ca65..275840d7758 100644
--- a/stan/math/rev/core/chainablestack.hpp
+++ b/stan/math/rev/core/chainablestack.hpp
@@ -3,6 +3,8 @@
 
 #include <stan/math/rev/core/autodiffstackstorage.hpp>
 
+#include <iostream>
+
 namespace stan {
 namespace math {
 
@@ -28,6 +30,28 @@ static inline ChainableStack::AutodiffStackStorage* init() {
 }
 */
 
+namespace {
+
+struct ad_tape_initializer {
+  ad_tape_initializer() {
+    std::cout << "Initializer created" << std::endl;
+    tape_ = stan::math::ChainableStack::init();
+  }
+  ~ad_tape_initializer() {
+    // std::lock_guard<std::mutex> cout_lock(cout_mutex_init);
+    std::cout << "Initializer destructed" << std::endl;
+  }
+
+  typedef ChainableStack::AutodiffStackStorage* tape_ptr_t;
+  thread_local static tape_ptr_t tape_;
+};
+
+thread_local ad_tape_initializer::tape_ptr_t ad_tape_initializer::tape_
+    = ChainableStack::init();
+ad_tape_initializer global_initializer;
+
+}  // namespace
+
 }  // namespace math
 }  // namespace stan
 #endif

From 00c99a036d85ccd91688a47db3a979acb755d295 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 16 Feb 2019 07:42:37 +0100
Subject: [PATCH 19/20] remove cout which is unexpected for cmdstan and fully
 disable init calls

---
 stan/math/rev/core/chainablestack.hpp      | 4 ++--
 stan/math/rev/core/chainablestack_inst.cpp | 2 ++
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/stan/math/rev/core/chainablestack.hpp b/stan/math/rev/core/chainablestack.hpp
index 275840d7758..9407270b452 100644
--- a/stan/math/rev/core/chainablestack.hpp
+++ b/stan/math/rev/core/chainablestack.hpp
@@ -34,12 +34,12 @@ namespace {
 
 struct ad_tape_initializer {
   ad_tape_initializer() {
-    std::cout << "Initializer created" << std::endl;
+    // std::cout << "Initializer created" << std::endl;
     tape_ = stan::math::ChainableStack::init();
   }
   ~ad_tape_initializer() {
     // std::lock_guard<std::mutex> cout_lock(cout_mutex_init);
-    std::cout << "Initializer destructed" << std::endl;
+    // std::cout << "Initializer destructed" << std::endl;
   }
 
   typedef ChainableStack::AutodiffStackStorage* tape_ptr_t;
diff --git a/stan/math/rev/core/chainablestack_inst.cpp b/stan/math/rev/core/chainablestack_inst.cpp
index 24b558a922c..e78fdfdfa94 100644
--- a/stan/math/rev/core/chainablestack_inst.cpp
+++ b/stan/math/rev/core/chainablestack_inst.cpp
@@ -7,8 +7,10 @@ namespace stan {
 namespace math {
 namespace internal {
 
+/*
 const ChainableStack::AutodiffStackStorage* __chainable_stack
     = ChainableStack::init();
+*/
 
 }  // namespace internal
 }  // namespace math

From 1178f5f05e831b34eb53a837a38d9e6847837050 Mon Sep 17 00:00:00 2001
From: Sebastian Weber <sdw.post@waebers.de>
Date: Sat, 23 Feb 2019 14:30:45 +0100
Subject: [PATCH 20/20] add unit test which can serve as microbenchmark

---
 .../integrate_ode_rk45_stress_test.cpp        | 56 +++++++++++++++++++
 1 file changed, 56 insertions(+)
 create mode 100644 test/unit/math/rev/mat/functor/integrate_ode_rk45_stress_test.cpp

diff --git a/test/unit/math/rev/mat/functor/integrate_ode_rk45_stress_test.cpp b/test/unit/math/rev/mat/functor/integrate_ode_rk45_stress_test.cpp
new file mode 100644
index 00000000000..e90fed53e3d
--- /dev/null
+++ b/test/unit/math/rev/mat/functor/integrate_ode_rk45_stress_test.cpp
@@ -0,0 +1,56 @@
+#include <stan/math/rev/mat.hpp>
+#include <boost/random.hpp>
+#include <stan/math/prim/scal/prob/lognormal_rng.hpp>
+#include <gtest/gtest.h>
+// very small michaelis menten example
+#include <test/unit/math/rev/arr/functor/coupled_mm.hpp>
+#include <test/unit/util.hpp>
+#include <vector>
+
+// test which triggers the too much work exception from odeint
+TEST(StanOde_tooMuchWork_test, rk45_coupled_mm) {
+  coupled_mm_ode_fun f_;
+
+  boost::ecuyer1988 rng;
+
+  // initial value and parameters from model definition
+
+  double t0 = 0;
+
+  std::vector<double> ts_long;
+  ts_long.push_back(1E3);
+
+  std::vector<double> ts_short;
+  ts_short.push_back(1);
+
+  std::vector<double> data;
+
+  std::vector<int> data_int;
+
+  typedef stan::math::var scalar_type;
+
+  for (std::size_t i = 0; i < 10; i++) {
+    stan::math::start_nested();
+
+    std::vector<scalar_type> theta_v(4);
+
+    theta_v[0] = stan::math::lognormal_rng(1, 2, rng);
+    theta_v[1] = stan::math::lognormal_rng(-1, 2, rng);
+    theta_v[2] = stan::math::lognormal_rng(-1, 2, rng);
+    theta_v[3] = stan::math::lognormal_rng(-2, 2, rng);
+
+    std::vector<scalar_type> y0_v(2);
+    y0_v[0] = stan::math::lognormal_rng(5, 2, rng);
+    y0_v[1] = stan::math::lognormal_rng(-1, 2, rng);
+
+    std::vector<std::vector<scalar_type> > res
+        = stan::math::integrate_ode_rk45(f_, y0_v, t0, ts_long, theta_v, data,
+                                         data_int, 0, 1E-6, 1E-6, 1000000000);
+
+    stan::math::grad(res[0][0].vi_);
+
+    stan::math::recover_memory_nested();
+  }
+
+  stan::math::print_stack(std::cout);
+}