在C++中使用OpenCL求解常微分方程

• 下载源代码 - 201 KB

引言

在本文中，我们将展示如何使 odeint 与 VexCL 协同工作。odeint 是一个用 C++ 数值求解常微分方程 (ODE) 的库。ODE 在许多科学领域都很重要，因此 odeint 有着广泛的应用。VexCL 是一个用于 OpenCL 的高级 C++ 库。其主要特点是表达式模板，这大大简化了数值问题代码的编写方式。通过使用 OpenCL，生成的代码可以在 GPU 上运行，也可以在多个核心上并行运行。这两个库都已在 CodeProject 上介绍过

• odeint v2 - 在 C++ 中求解常微分方程
• VexCL: 用于 OpenCL 的向量表达式模板库

注意：本文不介绍 odeint 和常微分方程。如果您不熟悉这些内容，请先阅读 odeint 文章！

注意：VexCL 需要 C++11 特性！因此您必须在启用 C++11 支持的情况下进行编译。

odeint 提供了一种机制，允许用户更改执行基本数值计算（加法、乘法等）的方式。这种机制由 state_type、algebra 和 operations 的组合组成。state_type 在此表示 ODE 的状态，通常是像 std::vector、std::array 这样的向量类型。例如：

std::array< double , 3 > x1 , x2;
// initialize x1, x2
odeint::range_algebra algebra;
double dt = 0.1;
algebra.for_each2( x1 , x2 , default_operations::scale_sum1( dt ) );
// computes x1 = dt * x2 for all elements of x1 and x2

代数负责迭代状态的所有元素，而操作负责基本操作。在上面的示例中，使用了 for_each2，这意味着迭代了两种状态类型。操作是 scale_sum1，它简单地计算 x1[i] = dt * x2[i]。odeint 提供了一组预定义的代数：

• range_algebra：默认代数，适用于 Boost.Ranges
• array_algebra：用于 boost::array 的专用代数
• fusion_algebra：用于编译时序列的代数，例如 boost::fusion::vector、std::tuple 等，请参见 Boost.Fusion
• vector_space_algebra：用于向量空间类型的代数，将所有计算直接重定向到操作。
• thrust_algebra：用于 Thrust 的代数

许多向量和矩阵类型库都为基本操作提供了表达式模板。例如 Boost.Ublas、MTL4 和 VexCL。此类库不需要自己的代数，但可以与 vector_space_algebra 和 default_operations 一起使用，后者直接在矩阵或向量类型上调用操作。在这种情况下，您只需调整 odeint 的大小调整机制。本文将介绍如何为 VexCL 完成此操作。Boost.Ublas 和 MTL4 的调整与 VexCL 的调整非常相似。

调整 VexCL

VexCL 引入了几种存在于 OpenCL 设备上的向量类型。主要的向量类型是 vex::vector，它是 std::vector 的经典模拟。vex::vector 也可以在多个设备上进行拆分。VexCL 的主要优点之一是它支持表达式模板。例如，可以编写如下代码：

vex::vector< double > x() , y() , z();
// initialize x, y, z
z = 0.125 * ( x * x + 2.0 * x * y + y * y );

此示例中的第二行延迟地创建一个表达式模板，该模板在调用赋值运算符时进行求值。优点是显然避免了临时变量，并且不会损失任何性能。

由于 VexCL 已经支持表达式模板，因此可以直接与 odeint 的 vector_space_algebra 一起使用。无需引入额外的代数或新操作。上一节中的示例可以写成：

vex::vector< double > x , y;
vector_space_algebra algebra;
double dt = 0.1;
algebra.for_each2( x , y , default_operations::scale_sum1( dt );

为了将 vex::vector 与 odeint 一起使用，只需调整 VexCL 的大小调整机制以适应 odeint。odeint 中的大小调整是必要的，因为许多求解器都需要临时状态类型。这些状态类型需要由 odeint 的大小调整机制构建和初始化。

odeint 的大小调整机制由三个类模板组成。这些类是 is_resizeable<>，它是一个简单的元函数，告诉 odeint 该类型是否真的可调整大小。第二个类是 same_size_impl<>，它有一个静态方法 same_size，接受两个 state_type 作为参数，并返回这两个类型是否具有相同的大小。第三个类是 resize_impl<>，它执行实际的大小调整。这些类都有默认实现，并且可以针对任何类型进行特化。

对于 VexCL，特化如下：

template< typename T >
struct is_resizeable< vex::vector< T > > : boost::true_type { };

template< typename T >
struct resize_impl< vex::vector< T > , vex::vector< T > >
{
    static void resize( vex::vector< T > &x1 , const vex::vector< T > &x2 )
    {
	x1.resize( x2.queue_list() , x2.size() );
    }
};

template< typename T >
struct same_size_impl< vex::vector< T > , vex::vector< T > >
{
    static bool same_size( const vex::vector< T > &x1 , const vex::vector< T > &x2 )
    {
	return x1.size() == x2.size();
    }
};

仅此而已。有了这些特化，就可以使用 VexCL 和 odeint 了。当然，这些特化已经在 odeint 中定义。您只需包含：

#include <boost/numeric/odeint/external/vexcl/vexcl_resizing.hpp>

VexCL 还有一个多向量，它打包了多个 vex::vector 实例，并允许同步操作所有这些实例。多向量的大小调整特化与 vex::vector 非常相似，并且也包含在上面的头文件中。

示例

只有当试图解决大型问题（可以并行解决许多子问题）时，GPU 的强大功能才能发挥作用。对于 ODEs，需要大约 10000 个耦合的 ODEs 才能从 GPU 获得比 CPU 更好的性能。本节将介绍两个典型的大型 ODEs 示例。

在第一个例子中，我们将使用洛伦兹系统，并研究其对其中一个参数的依赖性。洛伦兹系统是一个由三个耦合 ODE 组成的系统，在较大的参数范围内表现出混沌行为。ODE 为：

dx / dt = -sigma * ( x - y )
dy / dt = R * x - y - x * z
dz / dt = - b * z + x * y

我们将研究对参数的依赖性R. 因此，我们创建了一大组这样的系统（每个系统具有不同的参数 R），将它们全部打包到一个系统中，并在 GPU 上同时求解它们。洛伦兹系统是一个由三个耦合常微分方程组成的系统。如果我们要求解 N 个这样的系统，则总状态有 3*N 个条目。我们可以将每个分量分别打包到一个 VexCL 向量中。但由三个子向量组成的多向量更适合此问题。typedef 如下：

typedef vex::vector< double > vector_type;
typedef vex::multivector< double, 3 > state_type;

这里需要 vector_type 来存储参数R. state_type 的子向量可以通过以下方式访问：

state_type X;
// initialize X
auto &x = X(0);
auto &y = X(1);
auto &z = X(2);

因此，N 个 Lorenz 系统的所有 x 分量都在 X(0) 中，所有 y 分量都在 X(1) 中，所有 z 分量都在 X(2) 中。现在，我们实现系统函数。此函数表示 ODE，并由 odeint 用于求解 ODE。系统需要是一个具有三个参数的函数对象（一个函数子或一个普通 C 函数）。其签名为 void( const state_type&, state_type&, time_type )。第一个参数是输入参数，表示 ODE 的当前状态；第二个参数是输出参数，用于存储 ODE 的 RHS。第三个参数只是时间。如上所述，VexCL 支持数值计算的表达式模板。通过使用表达式模板，系统函数变得非常简单。

const double sigma = 10.0;
const double b = 8.0 / 3.0;

struct sys_func
{
    const vector_type &R;

    sys_func( const vector_type &_R ) : R( _R ) { }

    void operator()( const state_type &x , state_type &dxdt , double t ) const
    {
	dxdt(0) = -sigma * ( x(0) - x(1) );
	dxdt(1) = R * x(0) - x(1) - x(0) * x(2);
	dxdt(2) = - b * x(2) + x(0) * x(1);
    }
};

请注意，系统函数包含一个用于所有参数 R 的向量。系统函数中的每一行都计算向量的所有 N 个元素的整个集合的表达式。原则上这就是全部。我们现在可以实例化一个 odeint 求解器并求解 ODE。一个完整的 main 程序可能如下所示：

// setup the opencl context
vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
std::cout << ctx << std::endl;

// set up number of systems, time step and integration time
const size_t n = 1024 * 1024;
const double dt = 0.01;
const double t_max = 100.0;

// initialize R
double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
std::vector<double> r( n );
for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
vector_type R( ctx.queue() , r );

// initialize the state of the lorenz ensemble
state_type X(ctx.queue(), n);
X(0) = 10.0;
X(1) = 10.0;
X(2) = 10.0;

// instantiate a stepper
runge_kutta4<
    state_type , double , state_type , double ,
    odeint::vector_space_algebra , odeint::default_operations
    > stepper;

// solve the system
integrate_const( stepper , sys_func( R ) , X , 0.0 , t_max , dt );

如您所见，此处使用了 odeint 的 vector_space_algebra 和默认操作集。

作为第二个例子，我们选择了一个耦合相振荡器链。相振荡器是通常振荡器的一个非常简化的版本，其状态由一个 2π 周期变量描述。如果单个相振荡器未耦合，则其相位 φ 由线性增长 dφ / dt = ω 描述，其中 ω 是相速度。因此，只有当两个或更多振荡器耦合时，才能观察到有趣的行为。事实上，耦合相振荡器系统是涌现系统的一个突出例子，其中耦合系统表现出比其组成部分更复杂的行为。

我们在此分析的具体例子是

dφ(i) / dt = ω(i) + sin( dφ(i+1) - dφ(i) ) + sin( dφ(i) - dφ(i-1) )

请注意，dφ(i) 是时间的函数，这里的参数 i 表示链中的第 i 个相位。为了在 GPU 上高效地实现此类方程，Denis 在引入某种广义模板方面做得非常出色。我们问题的模板由以下代码生成：

extern const char oscillator_body[] = "return sin(X[-1] - X[0]) + sin(X[0] - X[1]);";
vex::StencilOperator< double, 3, 1, oscillator_body > S( queue_list() );

第一行简单地生成内核中完成的基本操作的 OpenCL 字符串。第二行实例化模板运算符。它可以通过施加 S(x) 应用于向量 x。相位振荡器链的完整系统函数如下：

extern const char oscillator_body[] = "return sin(X[-1] - X[0]) + sin(X[0] - X[1]);";

struct sys_func
{

    const state_type ω
    vex::StencilOperator< double, 3, 1, oscillator_body > S;

    sys_func( const state_type &_omega )
        : omega( _omega ) , S( _omega.queue_list() )
    { }

    void operator()( const state_type &x , state_type &dxdt , value_type t ) const
    {
        dxdt = omega + S( x );
    }
};

请再次注意，代码是多么紧凑。上面系统使用 Thrust 的等效版本要大得多。事实上，此处介绍的 Lorenz 系统实现有 78 行代码，而 Thrust 版本有 145 行。此外，Thrust 还需要单独的代数以及单独的操作类型。

性能对比

本节我们将 VexCL 的性能与 Thrust 进行比较。Thrust 是一个用于 CUDA 的高级库，它为 CUDA 设备上的向量提供类似 STL 的接口。此外，通过使用 OpenMP，它可以轻松地将计算放在 CPU 的一个或多个核心上。Thrust 的 Lorenz 示例在 odeint 的教程中描述。Thrust 不提供表达式模板，但它有一个高级迭代器系统，可以让您以简单的方式编写数值表达式。尽管如此，它比 VexCL 更复杂。

下图显示了 Lorenz 系统示例在 Thrust 和 VexCL 几种配置下的性能。具体而言，它显示了 VexCL 在一个 GPU、两个 GPU、三个 GPU 和一个 CPU 核心上的性能。此外，还显示了 Thrust 在 GPU（仅支持单个 GPU 计算）和一个 CPU 上的性能。明显可见，Thrust 在单个 GPU 上优于 VexCL。但是，如果使用多个 GPU（并且安装在您的计算机上），VexCL 会更快。CPU 版本也是如此。此外，可以清楚地看到，对于计算时间相对较小的系统尺寸，VexCL 具有恒定的运行时间。这是因为 OpenCL 在运行时为 GPU 编译内核。（右侧面板显示 VexCL 相对于 Thrust 在 GPU 上的性能。因此，如果曲线高于 1，则 VexCL 更快，否则更慢。）

相振荡器链的性能结果显示在下一图中。对于链的 Thrust 版本，系统函数已使用 Thrust 的迭代器系统实现。有趣的是，VexCL 在这里优于 Thrust。当然，对于小系统尺寸，存在恒定的开销，Thrust 在这里表现更好。但对于大型系统，VexCL 变得更快。这可能是由于 Thrust 中的迭代器有其代价，但目前尚不清楚具体是哪个版本更快。

结论

我们已经展示了如何将 VexCL 适应 odeint，以及如何利用它在求解大型常微分方程时提高性能。这里的大型意味着系统规模（耦合 ODE 的数量）应该在 10000-100000 的数量级，才能看到与普通 CPU 相比合理的性能提升。对于大型系统，这种提升可以达到约 20 倍。

VexCL 的性能也与 Thrust 进行了比较。在一个示例中，求解了一个大型 ODE 集。结果表明，在这种情况下，Thrust 比 VexCL 快约 10%。在第二个示例中，研究了一个耦合相振荡器系统。在这里，VexCL 快了 1.2 倍。尽管如此，VexCL 的表达式模板是该库的一大优势。它使用户能够在几分钟内解决复杂问题，而使用原生 CUDA 或 OpenCL 代码甚至 Thrust 进行开发则需要更多时间。

参考文献

• odeint v2 - 在 C++ 中求解常微分方程 - 一篇关于 odeint v2 的 CodeProject 文章。
• VexCL：用于 OpenCL 的向量表达式模板库 - 一篇关于 VexCL 的 CodeProject 文章
• 混沌与 Scholarpedia 上的洛伦兹系统
• 相振荡器在 Scholarpedia 上
• www.odeint.com - odeint 的项目页面
• VexCL 在 github 上的项目页面
• Thrust - 光速代码

历史

• 2012年7月26日 - 初始版本