replace to_arena with eval

README.md

@@ -63,9 +63,9 @@ From the top level directory you can build and call the tests with the following
 
 ```bash
 # DEBUG build types enable 0g, ggdb3, and pretty printing helper functions in utils
-cmake -S . -DCMAKE_BUILD_TYPE=RELEASE  -DRICCATI_BUILD_TESTING=ON -DBUILD_TESTING=OFF
-cd build
-make tests_riccati_test && ./tests_riccati_test
+cmake -S . -B "build" -DCMAKE_BUILD_TYPE=RELEASE  -DBUILD_TESTING=ON
+cd build/tests
+make -j4 riccati_test && ctest
 ```
 
 ## Including

include/riccati/chebyshev.hpp

@@ -331,9 +331,9 @@ template <typename Vec1, typename Vec2, typename Allocator>
 inline auto interpolate(Vec1&& s, Vec2&& t, Allocator&& alloc) {
   const auto r = s.size();
   const auto q = t.size();
-  auto V = to_arena(alloc,
+  auto V = eval(alloc,
                     matrix_t<typename std::decay_t<Vec1>::Scalar>::Ones(r, r));
-  auto R = to_arena(alloc,
+  auto R = eval(alloc,
                     matrix_t<typename std::decay_t<Vec1>::Scalar>::Ones(q, r));
   for (std::size_t i = 1; i < static_cast<std::size_t>(r); ++i) {
     V.col(i).array() = V.col(i - 1).array() * s.array();
@@ -394,25 +394,25 @@ inline auto spectral_chebyshev(SolverInfo&& info, Scalar x0, Scalar h,
   using complex_t = std::complex<Scalar>;
   using vectorc_t = vector_t<complex_t>;
   auto x_scaled
-      = to_arena(alloc, riccati::scale(info.chebyshev_[niter].second, x0, h));
+      = eval(alloc, riccati::scale(info.chebyshev_[niter].second, x0, h));
   auto&& D = info.Dn(niter);
   auto ws = info.omega_fun_(x_scaled);
   auto gs = info.gamma_fun_(x_scaled);
-  auto D2 = to_arena(
+  auto D2 = eval(
       alloc, (4.0 / (h * h) * (D * D) + 4.0 / h * (gs.asDiagonal() * D)));
   D2 += (ws.array().square()).matrix().asDiagonal();
   const auto n = std::round(info.ns_[niter]);
-  auto D2ic = to_arena(alloc, matrix_t<complex_t>::Zero(n + 3, n + 1));
+  auto D2ic = eval(alloc, matrix_t<complex_t>::Zero(n + 3, n + 1));
   D2ic.topRows(n + 1) = D2;
   D2ic.row(n + 1) = 2.0 / h * D.row(D.rows() - 1);
-  auto ic = to_arena(alloc, vectorc_t::Zero(n + 1));
+  auto ic = eval(alloc, vectorc_t::Zero(n + 1));
   ic.coeffRef(n) = complex_t{1.0, 0.0};
   D2ic.row(n + 2) = ic;
-  auto rhs = to_arena(alloc, vectorc_t::Zero(n + 3));
+  auto rhs = eval(alloc, vectorc_t::Zero(n + 3));
   rhs.coeffRef(n + 1) = dy0;
   rhs.coeffRef(n + 2) = y0;
-  auto y1 = to_arena(alloc, D2ic.colPivHouseholderQr().solve(rhs));
-  auto dy1 = to_arena(alloc, 2.0 / h * (D * y1));
+  auto y1 = eval(alloc, D2ic.colPivHouseholderQr().solve(rhs));
+  auto dy1 = eval(alloc, 2.0 / h * (D * y1));
 
   return std::make_tuple(std::move(y1), std::move(dy1), std::move(x_scaled));
 }

include/riccati/evolve.hpp

@@ -58,10 +58,10 @@ inline auto osc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
   using complex_t = std::complex<Scalar>;
   using vectorc_t = vector_t<complex_t>;
   auto xi_scaled
-      = to_arena(alloc, scale(info.xn().array(), xi, init_stepsize).matrix());
+      = eval(alloc, scale(info.xn().array(), xi, init_stepsize).matrix());
   // Frequency and friction functions evaluated at n+1 Chebyshev nodes
-  auto omega_n = to_arena(alloc, info.omega_fun_(xi_scaled));
-  auto gamma_n = to_arena(alloc, info.gamma_fun_(xi_scaled));
+  auto omega_n = eval(alloc, info.omega_fun_(xi_scaled));
+  auto gamma_n = eval(alloc, info.gamma_fun_(xi_scaled));
   vectorc_t yeval;
   // TODO: Add this check to regular evolve
   if (sign * (xi + init_stepsize) > sign * xf) {
@@ -87,11 +87,11 @@ inline auto osc_evolve(SolverInfo &&info, Scalar xi, Scalar xf,
           = get_slice(x_eval, sign * xi, sign * (xi + init_stepsize));
       if (dense_size != 0) {
         auto x_eval_map = x_eval.segment(dense_start, dense_size);
-        auto x_eval_scaled = to_arena(
+        auto x_eval_scaled = eval(
             alloc,
             (2.0 / init_stepsize * (x_eval_map.array() - xi) - 1.0).matrix());
         auto Linterp = interpolate(info.xn(), x_eval_scaled, alloc);
-        auto fdense = to_arena(
+        auto fdense = eval(
             alloc, (Linterp * std::get<5>(osc_ret)).array().exp().matrix());
         yeval = std::get<6>(osc_ret).first * fdense
                 + std::get<6>(osc_ret).second * fdense.conjugate();
@@ -415,15 +415,15 @@ inline auto evolve(SolverInfo &&info, Scalar xi, Scalar xf,
         auto y_eval_map
             = Eigen::Map<vectorc_t>(yeval.data() + dense_start, dense_size);
         if (steptype) {
-          auto x_eval_scaled = to_arena(
+          auto x_eval_scaled = eval(
               alloc,
               (2.0 / h * (x_eval_map.array() - xcurrent) - 1.0).matrix());
           auto Linterp = interpolate(info.xn(), x_eval_scaled, alloc);
-          auto fdense = to_arena(alloc, (Linterp * un).array().exp().matrix());
+          auto fdense = eval(alloc, (Linterp * un).array().exp().matrix());
           y_eval_map
               = a_pair.first * fdense + a_pair.second * fdense.conjugate();
         } else {
-          auto xc_scaled = to_arena(
+          auto xc_scaled = eval(
               alloc, scale(info.chebyshev_[1].second, xcurrent, h).matrix());
           auto Linterp = interpolate(xc_scaled, x_eval_map, alloc);
           y_eval_map = Linterp * y_eval;

include/riccati/memory.hpp

@@ -2,6 +2,7 @@
 #define INCLUDE_RICCATI_MEMORY_ARENA_ALLOC_HPP
 
 #include <riccati/macros.hpp>
+#include <riccati/utils.hpp>
 #include <Eigen/Dense>
 #include <stdint.h>
 #include <cstdlib>

include/riccati/step.hpp

@@ -133,7 +133,7 @@ inline auto osc_step(SolverInfo &&info, OmegaVec &&omega_s, GammaVec &&gamma_s,
   using vectorc_t = vector_t<complex_t>;
   bool success = true;
   auto &&Dn = info.Dn();
-  auto y = to_arena(alloc, complex_t(0.0, 1.0) * omega_s);
+  auto y = eval(alloc, complex_t(0.0, 1.0) * omega_s);
   auto delta = [&](const auto &r, const auto &y) {
     return (-r.array() / (2.0 * (y.array() + gamma_s.array()))).matrix().eval();
   };
@@ -159,18 +159,18 @@ inline auto osc_step(SolverInfo &&info, OmegaVec &&omega_s, GammaVec &&gamma_s,
     prev_err = maxerr;
   }
   if (info.denseout_) {
-    auto u1 = to_arena(alloc, h / 2.0 * (info.integration_matrix_ * y));
-    auto f1 = to_arena(alloc, (u1).array().exp().matrix());
-    auto f2 = to_arena(alloc, f1.conjugate());
-    auto du2 = to_arena(alloc, y.conjugate());
+    auto u1 = eval(alloc, h / 2.0 * (info.integration_matrix_ * y));
+    auto f1 = eval(alloc, (u1).array().exp().matrix());
+    auto f2 = eval(alloc, f1.conjugate());
+    auto du2 = eval(alloc, y.conjugate());
     auto ap_top = (dy0 - y0 * du2(du2.size() - 1));
     auto ap_bottom = (y(y.size() - 1) - du2(du2.size() - 1));
     auto ap = ap_top / ap_bottom;
     auto am = (dy0 - y0 * y(y.size() - 1))
               / (du2(du2.size() - 1) - y(y.size() - 1));
-    auto y1 = to_arena(alloc, ap * f1 + am * f2);
+    auto y1 = eval(alloc, ap * f1 + am * f2);
     auto dy1
-        = to_arena(alloc, ap * y.cwiseProduct(f1) + am * du2.cwiseProduct(f2));
+        = eval(alloc, ap * y.cwiseProduct(f1) + am * du2.cwiseProduct(f2));
     Scalar phase = std::imag(f1(0));
     return std::make_tuple(success, y1(0), dy1(0), maxerr, phase, u1,
                            std::make_pair(ap, am));

include/riccati/stepsize.hpp

@@ -63,14 +63,14 @@ template <typename SolverInfo, typename FloatingPoint, typename Allocator>
 inline auto choose_osc_stepsize(SolverInfo&& info, FloatingPoint x0,
                                 FloatingPoint h, FloatingPoint epsilon_h,
                                 Allocator&& alloc) {
-  auto t = to_arena(alloc, riccati::scale(info.xp_interp(), x0, h));
-  auto s = to_arena(alloc, riccati::scale(info.xp(), x0, h));
+  auto t = eval(alloc, riccati::scale(info.xp_interp(), x0, h));
+  auto s = eval(alloc, riccati::scale(info.xp(), x0, h));
   // TODO: Use a memory arena for these
   auto ws = info.omega_fun_(s).eval();
   auto gs = info.gamma_fun_(s).eval();
-  auto omega_analytic = to_arena(alloc, info.omega_fun_(t));
+  auto omega_analytic = eval(alloc, info.omega_fun_(t));
   auto omega_estimate = info.L() * ws;
-  auto gamma_analytic = to_arena(alloc, info.gamma_fun_(t));
+  auto gamma_analytic = eval(alloc, info.gamma_fun_(t));
   auto gamma_estimate = info.L() * gs;
   FloatingPoint max_omega_err
       = (((omega_estimate - omega_analytic).array() / omega_analytic.array())
@@ -83,7 +83,7 @@ inline auto choose_osc_stepsize(SolverInfo&& info, FloatingPoint x0,
   FloatingPoint max_err = std::max(max_omega_err, max_gamma_err);
   if (max_err <= epsilon_h) {
     if (info.p_ != info.n_) {
-      auto xn_scaled = to_arena(alloc, riccati::scale(info.xn(), x0, h));
+      auto xn_scaled = eval(alloc, riccati::scale(info.xn(), x0, h));
       ws = info.omega_fun_(xn_scaled);
       gs = info.gamma_fun_(xn_scaled);
     }

include/riccati/utils.hpp

@@ -74,7 +74,7 @@ inline auto eval(T&& x) {
 
 template <typename T, Eigen::Index R, Eigen::Index C>
 inline Eigen::Matrix<T, R, C> eval(Eigen::Matrix<T, R, C>&& x) {
-  return x;
+  return std::move(x);
 }
 
 template <typename T, Eigen::Index R, Eigen::Index C>
@@ -89,7 +89,7 @@ inline const auto& eval(const Eigen::Matrix<T, R, C>& x) {
 
 template <typename T, Eigen::Index R, Eigen::Index C>
 inline Eigen::Array<T, R, C> eval(Eigen::Array<T, R, C>&& x) {
-  return x;
+  return std::move(x);
 }
 
 template <typename T, Eigen::Index R, Eigen::Index C>