使用 Boost.odeint 和 Boost.SIMD 加速 ODE 模拟

• 下载 odeint_simd.zip - 3.6 KB

引言

常微分方程（ODEs）出现在许多应用中，求解它们是科学和工程中的一个基本问题。通常，我们必须依赖数值方法来逼近这些解，因为在存在非线性等问题时，无法找到解析解。由于这是一个常见的问题，存在各种库提供数值求解 ODE 的标准算法。在这里，我们将使用 Boost.odeint 库 [1,2,3]。Boost.odeint 以一种通用方式实现，使其高度灵活。这种灵活性使我们能够轻松地利用现代处理器的 SIMD 功能，并比标准实现获得高达三倍的显著速度提升。为了利用这种计算能力，我们将使用另一个工具：NT2 的 SIMD 库 [4,5,6]（提议成为 Boost.SIMD，以下简称 Boost.SIMD）。odeint 的通用设计允许我们以透明且简单的方式使用 Boost.SIMD 提供的 SIMD 数据结构。从标准实现到 SIMD 就绪代码的编程工作量几乎为零，这显示了像 Boost.odeint 和 Boost.SIMD 这样通用、设计精良的库的重要性和强大功能。

常微分方程的集合

常微分方程（ODEs）通常表示为以下形式：

d x(t) / d t = f( x(t) , t )

其中 x 被称为 ODE 的状态（通常是高维的），t 是因变量，函数 f 定义了 ODE。 ODE 最突出的例子是动力系统，其中 t 代表时间，x(t) 是由 f(x,t) 引起的动力学产生的轨迹。在这里，我们将面临一种情况，即拥有一个动态系统集合，我们希望对所有这些系统进行积分以获得每个系统的轨迹。在这种设置下，所有系统都遵循相同的动力学方程 f(x,t)，这将使我们能够使用 SIMD 指令。一种可能的示例是蒙特卡洛模拟，但也可以是参数研究，其中 f(x,t) 依赖于某些参数，这些参数对于集合的每个实现都取不同的值。

在下面的性能测试中，我们将模拟一个 Roessler 系统集合。集合中的每个系统都具有三维状态空间 (x,y,z)，并遵循以下动力学方程：

dx/dt = -y - z
dy/dt = x + ay
dz/dt = b + z(x-c)

常数 a, b 和 c，我们将其设置为一些典型值：a=0.2, b=1.0, c=9.0。我们的集合将由 N 个这样的系统组成，并且我们将为每个系统使用不同的随机初始条件。

实现

我们将从基于 Boost.odeint 的直接模拟实现开始。然后，我们将使用 Boost.SIMD 添加 SIMD 功能。我们将看到结合这两个库所需的努力非常少，以及在我们的设置中可以获得的巨大速度提升。

模拟集合

使用 Boost.odeint，实现上述模拟相当简单。首先，我们必须选择单个 Roessler 系统状态的表示。由于问题的维度很小、固定且在编译时已知，我们为此选择了一个 boost::array。然后，我们将 f(x,t) 函数以函数对象的形式提供：

static const size_t dim = 3;  // roessler is 3D
typedef boost::array<double, dim> state_type;

//---------------------------------------------------------------------------
struct roessler_system {
    const double m_a, m_b, m_c;

    roessler_system(const double a, const double b, const double c)
        : m_a(a), m_b(b), m_c(c)
    {}

    void operator()(const state_type &x, state_type &dxdt, const double t) const
    {
        dxdt[0] = -1.0*x[1] - x[2];
        dxdt[1] = x[0] + m_a * x[1];
        dxdt[2] = m_b + x[2] * (x[0] - m_c);
    }
};

有了这些要素，我们已经可以模拟单个 Roessler 系统。下面的示例使用步长 dt=0.1 的四阶龙格-库塔算法沿轨迹进行迭代：

const int steps = 50;
const double dt = 0.01;
state_type x = {0.0, -1.0, 2.5};  // some "random" initial conditions
boost::numeric::odeint::runge_kutta4<state_type> stepper;
for( int step = 0; step<steps; ++step )
    stepper.do_step(roessler_system(0.2, 1.0, 9.0), x, step*dt, dt);

但是，我们想要模拟整个 Roessler 系统集合。一种方法是模拟第一个系统，然后生成新的初始条件，然后模拟下一个系统，依此类推。但通常您需要在每个模拟步骤中都拥有所有系统的状态，例如用于统计分析。因此，我们引入一个向量，该向量包含所有系统的状态，并在每个步骤中依次迭代它们：

typedef std::vector<state_type> state_vec;
const int N = 1024*1024;  // number of systems in ensemble

state_vec state(N);

// generate random initial conditions...

boost::numeric::odeint::runge_kutta4<state_type> stepper;
roessler_system roessler(0.2, 1.0, 9.0);
double t = 0.0;
for( int step = 0; step < steps; ++step )
{
    for(size_t i = 0; i < N; ++i)
    {
        stepper.do_step(roessler, state[i], t, dt);
    }
    t += dt;
}

请注意，使用 state_vec state(N)，总共分配了 N*3*8 字节的连续内存，没有任何可能导致性能损失的碎片。此外，odeint 一次只处理一个系统，并使用 boost::array 作为其内部临时内存（由 stepper 的模板参数定义）。因此，代码中不再进行额外的内存分配，这使得代码具有缓存友好性，并为编译器提供了出色的优化可能性。最后，do_step 方法通过引用获取状态并在原地计算结果，这意味着没有执行不必要的复制操作。

上面省略的唯一部分是随机初始条件的初始化。我们将使用 C++11 引入的 <random> 库来生成初始条件并将其放入我们的 state_vec 结构中：

// random initial conditions
std::vector<double> x(N), y(N), z(N);
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution_xy(-8.0, 8.0);
std::uniform_real_distribution<double> distribution_z(0.0, 20.0);
auto rand_xy = std::bind(distribution_xy, std::ref(generator));
auto rand_z = std::bind(distribution_z, std::ref(generator));
std::generate(x.begin(), x.end(), rand_xy);
std::generate(y.begin(), y.end(), rand_xy);
std::generate(z.begin(), z.end(), rand_z);

for(size_t i=0; i<n; ++i)
{
    state[i][0] = x[i];
    state[i][1] = y[i];
    state[i][2] = z[i];
}

现在，我们已经集齐了所有部分，可以运行一些模拟性能测试了。所有性能测试均使用 Intel Compiler v15.0 在 3.8 GHz（睿频模式）的 Intel Xeon E5-2690 CPU 上进行。下图显示了该实现（蓝色条）在 MFlops/s（顶部图）下的性能，针对不同的系统大小。从图中可以看出，性能非常好，对于所有大小的系统，都能达到每周期 1 Flops 以上。仅使用 FPU，最多只能达到每周期 1 Flops。我们的代码性能甚至略好，意味着 Intel 编译器已经通过自动矢量化使用了一些 SIMD 指令。然而，在这个例子中，我们可以做得更好，如下所述。

使用 SIMD 指令

从上面的代码开始，我们将使用 Boost.SIMD 库来提高性能。使用 SIMD（单指令多数据）时，处理器能够在单个时钟周期内对多个值执行浮点运算。现代处理器支持 AVX 指令集，该指令集一次最多允许四次加法和乘法。因此，在理想情况下，我们可以在现代处理器上将性能提高四倍。

借助 Boost.SIMD，利用 CPU 的计算能力相当简单。本质上，我们只需更改计算的基本数据类型：

typedef boost::simd::pack<double> simd_pack;
typedef boost::array<simd_pack, dim> state_type;
typedef std::vector< state_type, boost::simd::allocator< state_type > > state_vec;

我们使用 boost::simd::pack 而不是普通的 doubles。这是 Boost.SIMD 提供的 SIMD 寄存器的抽象。Boost.SIMD 已经提供了所有基本的 SIMD pack 操作，这意味着我们甚至不需要更改 roessler 函数！请注意，SIMD 指令需要对齐的内存。为了确保我们的状态向量中有对齐的内存，我们使用 Boost.SIMD 的 simd::allocator。剩下要做的就是考虑到我们现在一次性处理了多个系统，因此 state 需要更少的条目。SIMD pack 的大小取决于处理器，并且可以通过 simd_pack::static_size 在编译时访问。有了这些，我们就可以使用 SIMD 指令实现上述模拟：

static const size_t pack_size = simd_pack::static_size;
state_vec state(N/pack_size);

// Generate random initial conditions...

for(int step = 0; step < steps; step++)
{
    for(size_t i = 0; i < N/pack_size; ++i)
    {
        stepper.do_step(sys, state[i], t, dt);
    }
    t += dt;
}

就这样！Boost.SIMD 负责所有 SIMD 增强的计算，而 Boost.odeint 使用 SIMD 数据类型执行 ODE 集成。唯一需要直接处理实现 SIMD 特性的地方是生成随机初始条件：

// random initial conditions for SIMD
std::vector<double> x(N), y(N), z(N);
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution_xy(-8.0, 8.0);
std::uniform_real_distribution<double> distribution_z(0.0, 20.0);

auto rand_xy = std::bind(distribution_xy, std::ref(generator));
auto rand_z = std::bind(distribution_z, std::ref(generator));
std::generate(x.begin(), x.end(), rand_xy);
std::generate(y.begin(), y.end(), rand_xy);
std::generate(z.begin(), z.end(), rand_z);

state_vec state(N/pack_size);
for(size_t i=0; i<N/pack_size; ++i)
{
    for(size_t p=0; p<pack_size; ++p)
    {
        state[i][0][p] = x[i*pack_size+p];
        state[i][1][p] = y[i*pack_size+p];
        state[i][2][p] = z[i*pack_size+p];
    }
}

如您所见，我们使用 []-运算符显式填充 SIMD 类型。由于我们只能在编译时知道目标机器提供多少 SIMD 寄存器，因此我们使用一个循环遍历 pack 大小来填充寄存器。通过这一点，我们已经完全将 ODE 集合集成转换为利用 SIMD 指令。上面的图表显示了此实现的性能。结果令人印象深刻。在同一台 Xeon E5-2690 CPU 上，我们达到了近 16 GFlops/s，即每周期 4 Flops，而对于大约 32 KB（1024 个系统）的小系统大小，标准实现大约为每周期 1.2 Flops。增加系统大小会导致性能明显下降，因为如果系统变得太大，它将无法放入 L1 缓存，CPU 不时需要等待 L2（对于非常大的系统，最终是 L3）缓存加载，这会降低性能。

然而，尽管 Flops/s 是衡量实现效率的一个很好的指标，但真正的性能只能通过测量模拟的实际运行时间来获得。右侧的图表显示了 SIMD 实现相对于标准实现的速度提升，以模拟的实际运行时间改进来衡量。如您所见，对于小型系统，我们达到了三倍的性能增益，对于大型系统，达到了 2.4 倍。这些数字略低于 Flops/s 的改进，但仍然相当可观。请注意，Flops 是由性能计数器测量的，可能并不 100% 准确。因此，在判断实现性能时，应始终依赖实际运行时间测量。

总结和结论

本文展示了如何通过结合 Boost.odeint 和 Boost.SIMD，利用 C++ 的强大功能发挥现代 CPU 的全部潜力。我们能够以几乎零编程工作量实现高达三倍的速度提升。这表明像 Boost.odeint 和 Boost.SIMD 这样设计良好、通用的 C++ 库的重要性和适用性。Roessler 系统集合的示例源代码可以轻松地推广到其他情况，如参数研究或更复杂系统集合。Boost.SIMD 提供了对现代 CPU SIMD 功能的简单而良好的抽象，并且应该考虑用于任何性能关键的数值算法。但是，请注意，许多算法受限于缓存而不是 Flops。在这种情况下，在任何 SIMD 调优有意义之前，应专注于改进缓存使用。

致谢

我感谢 Joel Falcou 指出内存对齐要求并对本文提供进一步详细反馈。此外，我感谢 Denis Demidov 校对文本。

这项工作得到了欧盟委员会通过“神经工程变革技术”（NETT）玛丽·居里初始培训网络的支持，项目编号为 289146。

资源

源代码

模拟的源代码可供下载：odeint_simd.zip - 3.6 KB

此外，性能测试已成为 Boost.odeint 的一部分，可在 Github 上找到。

参考文献

[1] Karsten Ahnert and Mario Mulansky, odeint - Solving Ordinary Differential Equations in C++, AIP Conf. Proc. 1389 (ICNAAM 2011), pp. 1586-1589 (2011)

[2] odeint-v2 on Codeproject: https://codeproject.org.cn/Articles/268589/odeint-v-Solving-ordinary-differential-equations

[3] odeint source code on Github: https://github.com/headmyshoulder/odeint-v2

[4] Pierre Estérie, Joel Falcou, Mathias Gaunard, Jean-Thierry Lapresté, Boost.SIMD: Generic Programming for portable Simdization, Proceedings of the 2014 Workshop on programming models for SIMD/Vector processing (2014)

[5] NT² SIMD library: http://www.numscale.com/boost-simd/

[6] NT² library on Github: https://github.com/MetaScale/nt2