1 /*
2 Copyright 2011-2012 Karsten Ahnert
3 Copyright 2013 Mario Mulansky
4
5 Distributed under the Boost Software License, Version 1.0.
6 (See accompanying file LICENSE_1_0.txt or
7 copy at http://www.boost.org/LICENSE_1_0.txt)
8 */
9
10
11 /*
12 * Solves many relaxation equations dxdt = - a * x in parallel and for different values of a.
13 * The relaxation equations are completely uncoupled.
14 */
15
16 #include <thrust/device_vector.h>
17
18 #include <boost/ref.hpp>
19
20 #include <boost/numeric/odeint.hpp>
21 #include <boost/numeric/odeint/external/thrust/thrust.hpp>
22
23
24 using namespace std;
25 using namespace boost::numeric::odeint;
26
27 // change to float if your GPU does not support doubles
28 typedef double value_type;
29 typedef thrust::device_vector< value_type > state_type;
30 typedef runge_kutta4< state_type , value_type , state_type , value_type > stepper_type;
31
32 struct relaxation
33 {
34 struct relaxation_functor
35 {
36 template< class T >
37 __host__ __device__
operator ()relaxation::relaxation_functor38 void operator()( T t ) const
39 {
40 // unpack the parameter we want to vary and the Lorenz variables
41 value_type a = thrust::get< 1 >( t );
42 value_type x = thrust::get< 0 >( t );
43 thrust::get< 2 >( t ) = -a * x;
44 }
45 };
46
relaxationrelaxation47 relaxation( size_t N , const state_type &a )
48 : m_N( N ) , m_a( a ) { }
49
operator ()relaxation50 void operator()( const state_type &x , state_type &dxdt , value_type t ) const
51 {
52 thrust::for_each(
53 thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_a.begin() , dxdt.begin() ) ) ,
54 thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_a.end() , dxdt.end() ) ) ,
55 relaxation_functor() );
56 }
57
58 size_t m_N;
59 const state_type &m_a;
60 };
61
62 const size_t N = 1024 * 1024;
63 const value_type dt = 0.01;
64
main(int arc,char * argv[])65 int main( int arc , char* argv[] )
66 {
67 // initialize the relaxation constants a
68 vector< value_type > a_host( N );
69 for( size_t i=0 ; i<N ; ++i ) a_host[i] = drand48();
70 state_type a = a_host;
71
72 // initialize the intial state x
73 state_type x( N );
74 thrust::fill( x.begin() , x.end() , 1.0 );
75
76 // integrate
77 relaxation relax( N , a );
78 integrate_const( stepper_type() , boost::ref( relax ) , x , 0.0 , 10.0 , dt );
79
80 return 0;
81 }
82