Computations in single precision

Added in version 3.2.0.

In previous tutorials we saw how heyoka, in addition to the standard double precision, also supports computations in extended precision and arbitrary precision. Starting with version 3.2.0, heyoka supports also computations in single precision.

Single-precision computations can lead to substantial performance benefits when high accuracy is not required. In particular, single-precision batch mode can use a SIMD width twice larger than double precision, leading to an increase by a factor of 2 of the computational throughput. In scalar computations, the use of single precision reduces by half the memory usage with respect to double precision, which can help alleviating performance issues in large ODE systems. This can be particularly noticeable in applications such as neural ODEs.

In C++, single-precision values are usually represented via the standard floating-point type float. Correspondingly, and similarly to what explained in the extended precision tutorial, single-precision computations are activated by passing the float template parameter to functions and classes in the heyoka API.

A simple example

In order to verify that heyoka indeed is able to work in single precision, we will be monitoring the evolution of the energy constant in a low-precision numerical integration of the simple pendulum.

Let us begin as usual with the definition of the dynamical equations and the creation of the integrator object:

    // Create the symbolic variables x and v.
    auto [x, v] = make_vars("x", "v");

    // Create the integrator object
    // in single precision.
    auto ta = taylor_adaptive<float>{// Definition of the ODE system:
                                     // x' = v
                                     // v' = -9.8 * sin(x)
                                     {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                     // Initial conditions
                                     // for x and v.
                                     {-1.f, 0.f}};

In order to activate single precision, we created an integrator object of type taylor_adaptive<float> - that is, we specified float, instead of the usual double, as the (only) template parameter for the taylor_adaptive class template. Note that we specified a single-precision initial state via the use of the f suffix for the numerical constants. Note also that, when operating in single precision, all numerical values encapsulated in an integrator are represented in single precision - this includes not only the state vector, but also the time coordinate, the tolerance, the Taylor coefficients, etc. Similarly to double-precision integrators, the default value of the tolerance is the machine epsilon of float.

Next, we define a small helper function that will allow us to monitor the evolution of the energy constant throughout the integration:

    // Create a small helper to compute the energy constant
    // from the state vector.
    auto compute_energy = [](const auto &sv) {
        using std::cos;

        return (sv[1] * sv[1]) / 2 + 9.8 * (1 - cos(sv[0]));
    };

Before starting the integration, we compute and store the initial energy for later use:

    // Compute and store the intial energy.
    const auto orig_E = compute_energy(ta.get_state());

We can now begin a step-by-step integration. At the end of each step, we will be computing and printing to screen the relative energy error:

    // Integrate for a few timesteps.
    for (auto i = 0; i < 20; ++i) {
        using std::abs;

        ta.step();

        std::cout << "Relative energy error: " << abs((orig_E - compute_energy(ta.get_state())) / orig_E) << '\n';
    }

Relative energy error: 1.48183e-07
Relative energy error: 5.29227e-08
Relative energy error: 6.08611e-08
Relative energy error: 1.79937e-07
Relative energy error: 1.74645e-07
Relative energy error: 2.24921e-07
Relative energy error: 2.4609e-07
Relative energy error: 1.1643e-07
Relative energy error: 1.79937e-07
Relative energy error: 1.40245e-07
Relative energy error: 2.54029e-07
Relative energy error: 1.84899e-07
Relative energy error: 1.83245e-07
Relative energy error: 1.56122e-07
Relative energy error: 2.22275e-07
Relative energy error: 1.61414e-07
Relative energy error: 2.11691e-07
Relative energy error: 2.88428e-07
Relative energy error: 2.93721e-07
Relative energy error: 1.82583e-07

The console output indeed confirms that energy is conserved at the level of the epsilon of the single-precision format (that is, \(\sim 10^{-7}\)).

Other classes and functions

Besides the adaptive integrator, several other classes and functions in heyoka can be used in single precision.

The event classes, for instance, can be constructed in single precision by passing float as the template parameter (instead of double). Note that the precision of an event must match the precision of the integrator object in which the event is used, otherwise an error will be produced at compilation time.

Full code listing

#include <cmath>
#include <iostream>

#include <heyoka/heyoka.hpp>

using namespace heyoka;

int main()
{
    // Create the symbolic variables x and v.
    auto [x, v] = make_vars("x", "v");

    // Create the integrator object
    // in single precision.
    auto ta = taylor_adaptive<float>{// Definition of the ODE system:
                                     // x' = v
                                     // v' = -9.8 * sin(x)
                                     {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                     // Initial conditions
                                     // for x and v.
                                     {-1.f, 0.f}};

    // Create a small helper to compute the energy constant
    // from the state vector.
    auto compute_energy = [](const auto &sv) {
        using std::cos;

        return (sv[1] * sv[1]) / 2 + 9.8 * (1 - cos(sv[0]));
    };

    // Compute and store the intial energy.
    const auto orig_E = compute_energy(ta.get_state());

    // Integrate for a few timesteps.
    for (auto i = 0; i < 20; ++i) {
        using std::abs;

        ta.step();

        std::cout << "Relative energy error: " << abs((orig_E - compute_energy(ta.get_state())) / orig_E) << '\n';
    }
}