Use connected lists for tape data structure

Instead of a single list use mutliple connected lists to store elements.
This allows to dynamically increase the size of tape without the need of
relocating elements.

benchmark/AlgorithmicComplexity.cpp

@@ -19,7 +19,7 @@ static void BM_NumericGausP(benchmark::State& state) {
   double p[] = {1, 2, 3, 4, 5};
   double dx[5] = {0, 0, 0, 0, 0};
   double dp[5] = {0, 0, 0, 0, 0};
-  clad::tape<clad::array_ref<double>> results = {};
+  clad::old_tape<clad::array_ref<double>> results = {};
   int dim = 5;
   results.emplace_back(dx, dim);
   results.emplace_back(dp, dim);

demos/CustomTypeNumDiff.cpp

@@ -134,7 +134,7 @@ int main() {
   // This is how we return the derivative with respect to all arguments.
   // The order of being placed in this tape should be the same as the order of
   // the arguments being passed to the function.
-  clad::tape<clad::array_ref<
+  clad::old_tape<clad::array_ref<
       double /*This should be the return value of the function you want to differentiate.*/>>
       grad = {};
   // Place the l-value reference of the variables in the tape.

include/clad/Differentiator/Differentiator.h

@@ -14,6 +14,7 @@
 #include "DynamicGraph.h"
 #include "FunctionTraits.h"
 #include "Matrix.h"
+#include "NewTape.h"
 #include "NumericalDiff.h"
 #include "Tape.h"
 
@@ -43,15 +44,16 @@ inline CUDA_HOST_DEVICE unsigned int GetLength(const char* code) {
 }
 
   /// Tape type used for storing values in reverse-mode AD inside loops.
-  template <typename T>
-  using tape = tape_impl<T>;
+template <typename T> using tape = new_tape_impl<T>;
 
-  /// Add value to the end of the tape, return the same value.
-  template <typename T, typename... ArgsT>
-  CUDA_HOST_DEVICE T push(tape<T>& to, ArgsT... val) {
-    to.emplace_back(std::forward<ArgsT>(val)...);
-    return to.back();
-  }
+template <typename T> using old_tape = tape_impl<T>;
+
+/// Add value to the end of the tape, return the same value.
+template <typename T, typename... ArgsT>
+CUDA_HOST_DEVICE T push(tape<T>& to, ArgsT... val) {
+  to.emplace_back(std::forward<ArgsT>(val)...);
+  return to.back();
+}
 
   /// Add value to the end of the tape, return the same value.
   /// A specialization for clad::array_ref types to use in reverse mode.

include/clad/Differentiator/NewTape.h

@@ -0,0 +1,133 @@
+#ifndef CLAD_DIFFERENTIATOR_NEWTAPE_H
+#define CLAD_DIFFERENTIATOR_NEWTAPE_H
+
+#include <cassert>
+#include <cstdio>
+#include <new>
+#include <type_traits>
+#include <utility>
+
+#include "clad/Differentiator/CladConfig.h"
+
+namespace clad {
+
+template <typename T> class Block {
+public:
+  T* data;
+  Block<T>* next;
+  Block<T>* prev;
+  std::size_t capacity;
+  using pointer = T*;
+  using iterator = pointer;
+
+  CUDA_HOST_DEVICE Block(std::size_t& capacity)
+      : next(nullptr), prev(nullptr), capacity(capacity),
+        data(allocateRawStorage(capacity)) {}
+
+  CUDA_HOST_DEVICE ~Block() {
+    destroy(block_begin(), block_end());
+    ::operator delete(
+        const_cast<void*>(static_cast<const volatile void*>(data)));
+  }
+  Block(const Block& other) = delete;
+  Block& operator=(const Block& other) = delete;
+
+  Block(Block&& other) = delete;
+  Block& operator=(const Block&& other) = delete;
+
+  CUDA_HOST_DEVICE iterator block_begin() { return data; }
+
+  CUDA_HOST_DEVICE iterator block_end() { return data + capacity; }
+
+  template <typename It> using value_type_of = decltype(*std::declval<It>());
+
+  CUDA_HOST_DEVICE T* allocateRawStorage(std::size_t capacity) {
+#ifdef __CUDACC__
+    return static_cast<T*>(::operator new(capacity * sizeof(T), std::nothrow));
+#else
+    return static_cast<T*>(::operator new(capacity * sizeof(T), std::nothrow));
+#endif
+  }
+
+  template <typename It>
+  static typename std::enable_if<
+      !std::is_trivially_destructible<value_type_of<It>>::value>::type
+  destroy(It B, It E) {
+    for (It I = E - 1; I >= B; --I)
+      I->~value_type_of<It>();
+  }
+
+  template <typename It>
+  static typename std::enable_if<
+      std::is_trivially_destructible<value_type_of<It>>::value>::type
+      CUDA_HOST_DEVICE
+      destroy(It B, It E) {}
+};
+
+template <typename T> class new_tape_impl {
+  Block<T>* m_cur_block = nullptr;
+  std::size_t m_capacity = 128;
+  std::size_t m_size = 0;
+  using pointer = T*;
+  using reference = T&;
+
+  using iterator = pointer;
+
+public:
+  new_tape_impl() = default;
+
+  ~new_tape_impl() {
+    while (m_cur_block) {
+      Block<T>* prev_block = m_cur_block->prev;
+      delete m_cur_block;
+      m_cur_block = prev_block;
+    }
+  }
+
+  new_tape_impl(new_tape_impl& other) = delete;
+  new_tape_impl operator=(new_tape_impl& other) = delete;
+
+  new_tape_impl(new_tape_impl&& other) = delete;
+  new_tape_impl& operator=(new_tape_impl&& other) = delete;
+
+  template <typename... ArgsT>
+
+  CUDA_HOST_DEVICE void emplace_back(ArgsT&&... args) {
+    if (!m_cur_block || m_size >= m_capacity) {
+      Block<T>* prev_block = m_cur_block;
+      m_cur_block = new Block<T>(m_capacity);
+      if (prev_block) {
+        prev_block->next = m_cur_block;
+        m_cur_block->prev = prev_block;
+      }
+      m_size = 0;
+    }
+    m_size += 1;
+    ::new (const_cast<void*>(static_cast<const volatile void*>(end())))
+        T(std::forward<ArgsT>(args)...);
+  }
+
+  [[nodiscard]] CUDA_HOST_DEVICE std::size_t size() const { return m_size; }
+
+  CUDA_HOST_DEVICE iterator end() { return m_cur_block->data + (m_size - 1); }
+
+  CUDA_HOST_DEVICE reference back() {
+    assert(m_size || m_cur_block->prev);
+    return *end();
+  }
+
+  CUDA_HOST_DEVICE void pop_back() {
+    assert(m_size || m_cur_block->prev);
+    end()->~T();
+    m_size -= 1;
+    if (m_size == 0) {
+      Block<T>* temp = m_cur_block;
+      m_cur_block = m_cur_block->prev;
+      temp->~Block<T>();
+      m_size = m_capacity;
+    }
+  }
+};
+} // namespace clad
+
+#endif // CLAD_DIFFERENTIATOR_NEWTAPE_H

test/NumericalDiff/PureCentralDiffCalls.C

@@ -75,7 +75,7 @@ int main() { // expected-no-diagnostics
   printf("Result is = %f\n", func1_res); // CHECK-EXEC: Result is = 2.000000
 
   // Gradients, derivative wrt all args
-  clad::tape<clad::array_ref<double>> grad = {};
+  clad::old_tape<clad::array_ref<double>> grad = {};
   grad.emplace_back(dx, 3);
   grad.emplace_back(&dy);
   grad.emplace_back(&dz);