odeint speaks OpenCL

In the previous post I wrote some words about VexCL, a library for OpenCL. This library can easily be used within odeint. I only needed to create an adapter for the odeint resizing mechanism and then it was done. VexCL supports expression template and therefore, it can be used with odeints vector_space_algebra and default_operations.Here is an example of a parameter study of the Lorenz system:

#include <iostream>
#include <vector>
 
#include <vexcl/vexcl.hpp>
 
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/vexcl/vexcl_resize.hpp>
 
typedef vex::vector< double >    vector_type;
typedef vex::multivector< double, 3 > state_type;
 
// system function for odeint
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
    {
        const double sigma = 10.0;
        const double b = 8.0 / 3.0;
	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);
    }
};
 
int main( int argc , char **argv )
{
    using namespace std;
    using namespace boost::numeric::odeint;
 
    // 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 = 100.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 , 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 );
 
    // output
    std::vector< double > res( 3 * n );
    vex::copy( X(0) , res );
    for( size_t i=0 ; i<n ; ++i )
        cout << r[i] << "\t" << res[i] << "\t" << "\n";
 
    return 0;
}

2 thoughts on “odeint speaks OpenCL

Leave a Reply

Your email address will not be published. Required fields are marked *