Squashed 'third_party/boostorg/odeint/' content from commit 6ff2719
Change-Id: If4892e29c1a5e6cf3a7aa51486a2725c251b0c7d
git-subtree-dir: third_party/boostorg/odeint
git-subtree-split: 6ff2719b6907b86596c3d43e88c1bcfdf29df560
diff --git a/examples/thrust/relaxation.cu b/examples/thrust/relaxation.cu
new file mode 100644
index 0000000..f1d9f3a
--- /dev/null
+++ b/examples/thrust/relaxation.cu
@@ -0,0 +1,81 @@
+/*
+ Copyright 2011-2012 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+
+/*
+ * Solves many relaxation equations dxdt = - a * x in parallel and for different values of a.
+ * The relaxation equations are completely uncoupled.
+ */
+
+#include <thrust/device_vector.h>
+
+#include <boost/ref.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/thrust/thrust.hpp>
+
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+// change to float if your GPU does not support doubles
+typedef double value_type;
+typedef thrust::device_vector< value_type > state_type;
+typedef runge_kutta4< state_type , value_type , state_type , value_type > stepper_type;
+
+struct relaxation
+{
+ struct relaxation_functor
+ {
+ template< class T >
+ __host__ __device__
+ void operator()( T t ) const
+ {
+ // unpack the parameter we want to vary and the Lorenz variables
+ value_type a = thrust::get< 1 >( t );
+ value_type x = thrust::get< 0 >( t );
+ thrust::get< 2 >( t ) = -a * x;
+ }
+ };
+
+ relaxation( size_t N , const state_type &a )
+ : m_N( N ) , m_a( a ) { }
+
+ void operator()( const state_type &x , state_type &dxdt , value_type t ) const
+ {
+ thrust::for_each(
+ thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_a.begin() , dxdt.begin() ) ) ,
+ thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_a.end() , dxdt.end() ) ) ,
+ relaxation_functor() );
+ }
+
+ size_t m_N;
+ const state_type &m_a;
+};
+
+const size_t N = 1024 * 1024;
+const value_type dt = 0.01;
+
+int main( int arc , char* argv[] )
+{
+ // initialize the relaxation constants a
+ vector< value_type > a_host( N );
+ for( size_t i=0 ; i<N ; ++i ) a_host[i] = drand48();
+ state_type a = a_host;
+
+ // initialize the intial state x
+ state_type x( N );
+ thrust::fill( x.begin() , x.end() , 1.0 );
+
+ // integrate
+ relaxation relax( N , a );
+ integrate_const( stepper_type() , boost::ref( relax ) , x , 0.0 , 10.0 , dt );
+
+ return 0;
+}