In the previous section the usage of odeint in combination with Thrust was shown. In this section we show how one can use OpenCL with odeint. The point of odeint is not to implement its own low-level data structures and algorithms, but to use high level libraries doing this task. Here, we will use the VexCL framework to use OpenCL. VexCL is a nice library for general computations and it uses heavily expression templates. With the help of VexCL it is possible to write very compact and expressive application.

	Note
	vexcl needs C++11 features! So you have to compile with C++11 support enabled.

To use VexCL one needs to include one additional header which includes the data-types and algorithms from vexcl and the adaption to odeint. Adaption to odeint means here only to adapt the resizing functionality of VexCL to odeint.

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

To demonstrate the use of VexCL we integrate an ensemble of Lorenz system. The example is very similar to the parameter study of the Lorenz system in the previous section except that we do not compute the Lyapunov exponents. Again, we vary the parameter R of the Lorenz system an solve a whole ensemble of Lorenz systems in parallel (each with a different parameter R). First, we define the state type and a vector type

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

The vector_type is used to represent the parameter R. The state_type is a multi-vector of three sub vectors and is used to represent. The first component of this multi-vector represent all x components of the Lorenz system, while the second all y components and the third all z components. The components of this vector can be obtained via

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

As already mentioned VexCL supports expression templates and we will use them to implement the system function for the Lorenz ensemble:

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);
    }
};

It's very easy, isn't it? These three little lines do all the computations for you. There is no need to write your own OpenCL kernels. VexCL does everything for you. Next we have to write the main application. We initialize the vector of parameters (R) and the initial state. Note that VexCL requires the vector_space_algebra, but that is automatically deduced and configured by odeint internally, so we only have to specify the state_type when instantiating the stepper and we are done:

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

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

// initialize R
double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
std::vector<double> x( n * 3 ) , 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;

// create a stepper
runge_kutta4< state_type > stepper;

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