diff --git a/README.md b/README.md index 95236ae..a11b7db 100644 --- a/README.md +++ b/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 diff --git a/include/riccati/chebyshev.hpp b/include/riccati/chebyshev.hpp index 2c39d93..3a30a2e 100644 --- a/include/riccati/chebyshev.hpp +++ b/include/riccati/chebyshev.hpp @@ -331,9 +331,9 @@ template 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::Scalar>::Ones(r, r)); - auto R = to_arena(alloc, + auto R = eval(alloc, matrix_t::Scalar>::Ones(q, r)); for (std::size_t i = 1; i < static_cast(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; using vectorc_t = vector_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::Zero(n + 3, n + 1)); + auto D2ic = eval(alloc, matrix_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)); } diff --git a/include/riccati/evolve.hpp b/include/riccati/evolve.hpp index a3f260b..b97053d 100644 --- a/include/riccati/evolve.hpp +++ b/include/riccati/evolve.hpp @@ -58,10 +58,10 @@ inline auto osc_evolve(SolverInfo &&info, Scalar xi, Scalar xf, using complex_t = std::complex; using vectorc_t = vector_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(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; diff --git a/include/riccati/memory.hpp b/include/riccati/memory.hpp index 70c65b8..b6b25ce 100644 --- a/include/riccati/memory.hpp +++ b/include/riccati/memory.hpp @@ -2,6 +2,7 @@ #define INCLUDE_RICCATI_MEMORY_ARENA_ALLOC_HPP #include +#include #include #include #include diff --git a/include/riccati/step.hpp b/include/riccati/step.hpp index 59d0d64..afd4ab2 100644 --- a/include/riccati/step.hpp +++ b/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; 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)); diff --git a/include/riccati/stepsize.hpp b/include/riccati/stepsize.hpp index d771342..06b537f 100644 --- a/include/riccati/stepsize.hpp +++ b/include/riccati/stepsize.hpp @@ -63,14 +63,14 @@ template 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); } diff --git a/include/riccati/utils.hpp b/include/riccati/utils.hpp index 358161e..882ae27 100644 --- a/include/riccati/utils.hpp +++ b/include/riccati/utils.hpp @@ -74,7 +74,7 @@ inline auto eval(T&& x) { template inline Eigen::Matrix eval(Eigen::Matrix&& x) { - return x; + return std::move(x); } template @@ -89,7 +89,7 @@ inline const auto& eval(const Eigen::Matrix& x) { template inline Eigen::Array eval(Eigen::Array&& x) { - return x; + return std::move(x); } template