Merge pull request #428 from bluescarni/pr/fixes

Compilation/testing fixes

benchmark/pendulum_jet.jl

@@ -0,0 +1,33 @@
+using TaylorIntegration
+
+H(x) = 0.5x[2]^2-cos(x[1])
+
+@taylorize function pendulum!(dx, x, p, t)
+    dx[1] = x[2]
+    dx[2] = -sin(x[1])
+end
+
+const varorder = 8
+using TaylorIntegration
+#ξ = set_variables("ξ", numvars=2, order=varorder)
+
+q0 = [1.3, 0.0]
+
+#q0TN = q0 .+ ξ
+q0TN = q0
+
+order = 20     #the order of the Taylor expansion wrt time
+abstol = 2e-16 #the absolute tolerance of the integration
+using Elliptic # we use Elliptic.jl to evaluate the elliptic integral K
+T = 4*Elliptic.K(sin(q0[1]/2)^2) #the libration period
+t0 = 0.0        #the initial time
+tmax = 6T       #the final time
+integstep = T/8 #the time interval between successive evaluations of the solution vector
+
+
+tv = t0:integstep:tmax # the times at which the solution will be evaluated
+xv = taylorinteg(pendulum!, q0TN, tv, order, abstol)
+
+using BenchmarkTools
+bV = @benchmark taylorinteg(pendulum!, q0TN, tv, order, abstol)
+println(bV)

doc/changelog.rst

@@ -12,6 +12,15 @@ New
 - Add relational and logical operators to the expression system
   (`#432 <https://github.com/bluescarni/heyoka/pull/432>`__).
 
+Fix
+~~~
+
+- Fix compilation on PPC64
+  (`#428 <https://github.com/bluescarni/heyoka/pull/428>`__).
+- Relax several tolerances in order to fix test failures reported
+  on FreeBSD
+  (`#428 <https://github.com/bluescarni/heyoka/pull/428>`__).
+
 5.0.0 (2024-06-13)
 ------------------
 

src/detail/tm_data.cpp

@@ -64,13 +64,25 @@ tm_data<T>::tm_data() = default;
 namespace
 {
 
+// LCOV_EXCL_START
+
 // Factorial implementations.
 template <typename F>
-F factorial(std::uint32_t n, long long)
+F factorial([[maybe_unused]] std::uint32_t n, long long)
 {
-    return boost::math::factorial<F>(boost::numeric_cast<unsigned>(n));
+#if defined(HEYOKA_ARCH_PPC)
+    if constexpr (std::same_as<F, long double>) {
+        throw std::invalid_argument("'long double' computations are not supported on this platform");
+    } else {
+#endif
+        return boost::math::factorial<F>(boost::numeric_cast<unsigned>(n));
+#if defined(HEYOKA_ARCH_PPC)
+    }
+#endif
 }
 
+// LCOV_EXCL_STOP
+
 #if defined(HEYOKA_HAVE_REAL128)
 
 template <>

src/llvm_state.cpp

@@ -150,6 +150,26 @@
 #include <heyoka/s11n.hpp>
 #include <heyoka/variable.hpp>
 
+// NOTE: logging here lhames' instructions on how to set up LLJIT
+// for parallel compilation of multiple modules.
+//
+//   auto J = LLJITBuilder()
+//              .setNumCompileThreads(<N>)
+//              .create();
+//   if (!J) { /* bail on error */ }
+//   (*J)->getIRTransformLayer().setTransform(
+//     [](ThreadSafeModule TSM, MaterializationResponsibility &R) -> Expected<ThreadSafeModule> {
+//       TSM.withModuleDo([](Module &M) {
+//         /* Apply your IR optimizations here */
+//       });
+//       return std::move(TSM);
+//     });
+//
+// Note that the optimisation passes in this approach are moved into the
+// transform layer. References:
+// https://discord.com/channels/636084430946959380/687692371038830597/1252428080648163328
+// https://discord.com/channels/636084430946959380/687692371038830597/1252118666187640892
+
 HEYOKA_BEGIN_NAMESPACE
 
 namespace detail

src/taylor_adaptive.cpp

@@ -450,13 +450,13 @@ void taylor_adaptive<T>::finalise_ctor_impl(sys_t vsys, std::vector<T> state,
 
     // Fetch the stepper.
     if (with_events) {
-        m_step_f = reinterpret_cast<i_data::step_f_e_t>(m_llvm.jit_lookup("step_e"));
+        m_step_f = reinterpret_cast<typename i_data::step_f_e_t>(m_llvm.jit_lookup("step_e"));
     } else {
-        m_step_f = reinterpret_cast<i_data::step_f_t>(m_llvm.jit_lookup("step"));
+        m_step_f = reinterpret_cast<typename i_data::step_f_t>(m_llvm.jit_lookup("step"));
     }
 
     // Fetch the function to compute the dense output.
-    m_d_out_f = reinterpret_cast<i_data::d_out_f_t>(m_llvm.jit_lookup("d_out_f"));
+    m_d_out_f = reinterpret_cast<typename i_data::d_out_f_t>(m_llvm.jit_lookup("d_out_f"));
 
     // Setup the vector for the Taylor coefficients.
     using su32_t = boost::safe_numerics::safe<std::uint32_t>;
@@ -545,9 +545,9 @@ taylor_adaptive<T>::taylor_adaptive(const taylor_adaptive &other)
       m_ed_data(other.m_ed_data ? std::make_unique<ed_data>(*other.m_ed_data) : nullptr)
 {
     if (m_ed_data) {
-        m_i_data->m_step_f = reinterpret_cast<i_data::step_f_e_t>(m_i_data->m_llvm.jit_lookup("step_e"));
+        m_i_data->m_step_f = reinterpret_cast<typename i_data::step_f_e_t>(m_i_data->m_llvm.jit_lookup("step_e"));
     } else {
-        m_i_data->m_step_f = reinterpret_cast<i_data::step_f_t>(m_i_data->m_llvm.jit_lookup("step"));
+        m_i_data->m_step_f = reinterpret_cast<typename i_data::step_f_t>(m_i_data->m_llvm.jit_lookup("step"));
     }
 }
 
@@ -613,9 +613,9 @@ void taylor_adaptive<T>::load_impl(Archive &ar, unsigned version)
 
         // Recover the function pointers.
         if (m_ed_data) {
-            m_i_data->m_step_f = reinterpret_cast<i_data::step_f_e_t>(m_i_data->m_llvm.jit_lookup("step_e"));
+            m_i_data->m_step_f = reinterpret_cast<typename i_data::step_f_e_t>(m_i_data->m_llvm.jit_lookup("step_e"));
         } else {
-            m_i_data->m_step_f = reinterpret_cast<i_data::step_f_t>(m_i_data->m_llvm.jit_lookup("step"));
+            m_i_data->m_step_f = reinterpret_cast<typename i_data::step_f_t>(m_i_data->m_llvm.jit_lookup("step"));
         }
         // LCOV_EXCL_START
     } catch (...) {

test/batch_event_detection.cpp

@@ -164,7 +164,7 @@ TEST_CASE("nte single step")
             REQUIRE(trig_times.size() == trig_times_batch[i].size());
             for (decltype(trig_times.size()) j = 0; j < trig_times.size(); ++j) {
                 REQUIRE(trig_times[j] == approximately(trig_times_batch[i][j], 1000.));
-                REQUIRE(v_vals[j] == approximately(v_vals_batch[i][j], 1000.));
+                REQUIRE(v_vals[j] == approximately(v_vals_batch[i][j], 10000.));
             }
         }
     };
@@ -248,7 +248,7 @@ TEST_CASE("te single step")
             REQUIRE(trig_times.size() == trig_times_batch[i].size());
             for (decltype(trig_times.size()) j = 0; j < trig_times.size(); ++j) {
                 REQUIRE(trig_times[j] == approximately(trig_times_batch[i][j], 1000.));
-                REQUIRE(v_vals[j] == approximately(v_vals_batch[i][j], 1000.));
+                REQUIRE(v_vals[j] == approximately(v_vals_batch[i][j], 10000.));
             }
         }
     };
@@ -1547,8 +1547,7 @@ TEST_CASE("te propagate_grid first step bug")
     }
 
     {
-        t_ev_t ev(
-            v, kw::callback = [](auto &, int, std::uint32_t) { return true; });
+        t_ev_t ev(v, kw::callback = [](auto &, int, std::uint32_t) { return true; });
 
         auto ta = taylor_adaptive_batch<double>{{prime(x) = v, prime(v) = -9.8 * sin(x)},
                                                 {0.05, 0.051, 0.052, 0.053, 0.025, 0.0251, 0.0252, 0.0253},

test/model_rotating.cpp

@@ -97,7 +97,7 @@ TEST_CASE("basic")
         double E = 0;
         cf(&E, ta.get_state().data(), nullptr, nullptr);
 
-        REQUIRE(E == approximately(E0));
+        REQUIRE(E == approximately(E0, 1000.));
 
         REQUIRE(0.5 * sum_sq({vx, vy, vz}) + model::rotating_potential(kw::omega = {.1, .2, .3})
                 == model::rotating_energy(kw::omega = {.1, .2, .3}));
@@ -138,7 +138,7 @@ TEST_CASE("basic")
         double E = 0;
         cf(&E, ta.get_state().data(), omega_vals.data(), nullptr);
 
-        REQUIRE(E == approximately(E0));
+        REQUIRE(E == approximately(E0, 1000.));
 
         REQUIRE(0.5 * sum_sq({vx, vy, vz}) + model::rotating_potential(kw::omega = {par[0], par[1], par[2]})
                 == model::rotating_energy(kw::omega = {par[0], par[1], par[2]}));

test/var_ode_sys.cpp

@@ -464,9 +464,9 @@ TEST_CASE("nonauto sys")
     ta_orig.propagate_until(3.);
 
     REQUIRE(ta_orig.get_state()[0]
-            == approximately(ta.get_state()[0] + ta.get_state()[4] * delta_par + ta.get_state()[5] * delta_tm));
+            == approximately(ta.get_state()[0] + ta.get_state()[4] * delta_par + ta.get_state()[5] * delta_tm, 1000.));
     REQUIRE(ta_orig.get_state()[1]
-            == approximately(ta.get_state()[1] + ta.get_state()[8] * delta_par + ta.get_state()[9] * delta_tm));
+            == approximately(ta.get_state()[1] + ta.get_state()[8] * delta_par + ta.get_state()[9] * delta_tm, 1000.));
 }
 
 // An order-2 test.