123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233 |
- #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_OPERATIONS_HPP_INCLUDED
- #define BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_OPERATIONS_HPP_INCLUDED
- namespace boost {
- namespace numeric {
- namespace odeint {
- #include <thrust/tuple.h>
- #include <thrust/iterator/zip_iterator.h>
- struct thrust_operations
- {
- template< class Fac1 = double , class Fac2 = Fac1 >
- struct scale_sum2
- {
- const Fac1 m_alpha1;
- const Fac2 m_alpha2;
- scale_sum2( const Fac1 alpha1 , const Fac2 alpha2 )
- : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) + m_alpha2 * thrust::get<2>(t);
- }
- };
- template< class Fac1 = double , class Fac2 = Fac1 >
- struct scale_sum_swap2
- {
- const Fac1 m_alpha1;
- const Fac2 m_alpha2;
- scale_sum_swap2( const Fac1 alpha1 , const Fac2 alpha2 )
- : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- typename thrust::tuple_element<0,Tuple>::type tmp = thrust::get<0>(t);
- thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) + m_alpha2 * thrust::get<2>(t);
- thrust::get<1>(t) = tmp;
- }
- };
- template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 >
- struct scale_sum3
- {
- const Fac1 m_alpha1;
- const Fac2 m_alpha2;
- const Fac3 m_alpha3;
- scale_sum3( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 )
- : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) { }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) +
- m_alpha2 * thrust::get<2>(t) +
- m_alpha3 * thrust::get<3>(t);
- }
- };
- template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 >
- struct scale_sum4
- {
- const Fac1 m_alpha1;
- const Fac2 m_alpha2;
- const Fac3 m_alpha3;
- const Fac4 m_alpha4;
- scale_sum4( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 , const Fac4 alpha4 )
- : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ){ }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) +
- m_alpha2 * thrust::get<2>(t) +
- m_alpha3 * thrust::get<3>(t) +
- m_alpha4 * thrust::get<4>(t);
- }
- };
- template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 ,
- class Fac4 = Fac3 , class Fac5 = Fac4 >
- struct scale_sum5
- {
- const Fac1 m_alpha1;
- const Fac2 m_alpha2;
- const Fac3 m_alpha3;
- const Fac4 m_alpha4;
- const Fac5 m_alpha5;
- scale_sum5( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 ,
- const Fac4 alpha4 , const Fac5 alpha5 )
- : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) ,
- m_alpha4( alpha4 ) , m_alpha5( alpha5 ) { }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) +
- m_alpha2 * thrust::get<2>(t) +
- m_alpha3 * thrust::get<3>(t) +
- m_alpha4 * thrust::get<4>(t) +
- m_alpha5 * thrust::get<5>(t);
- }
- };
- template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 ,
- class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 >
- struct scale_sum6
- {
- const Fac1 m_alpha1;
- const Fac2 m_alpha2;
- const Fac3 m_alpha3;
- const Fac4 m_alpha4;
- const Fac5 m_alpha5;
- const Fac6 m_alpha6;
- scale_sum6( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 ,
- const Fac4 alpha4 , const Fac5 alpha5 , const Fac6 alpha6 )
- : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) ,
- m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) { }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) +
- m_alpha2 * thrust::get<2>(t) +
- m_alpha3 * thrust::get<3>(t) +
- m_alpha4 * thrust::get<4>(t) +
- m_alpha5 * thrust::get<5>(t) +
- m_alpha6 * thrust::get<6>(t);
- }
- };
- template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 ,
- class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 >
- struct scale_sum7
- {
- const Fac1 m_alpha1;
- const Fac2 m_alpha2;
- const Fac3 m_alpha3;
- const Fac4 m_alpha4;
- const Fac5 m_alpha5;
- const Fac6 m_alpha6;
- const Fac7 m_alpha7;
- scale_sum7( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 ,
- const Fac4 alpha4 , const Fac5 alpha5 , const Fac6 alpha6 , const Fac7 alpha7 )
- : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) ,
- m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) { }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) +
- m_alpha2 * thrust::get<2>(t) +
- m_alpha3 * thrust::get<3>(t) +
- m_alpha4 * thrust::get<4>(t) +
- m_alpha5 * thrust::get<5>(t) +
- m_alpha6 * thrust::get<6>(t) +
- m_alpha7 * thrust::get<7>(t) ;
- }
- };
- template< class Fac1 = double >
- struct rel_error
- {
- const Fac1 m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt;
- rel_error( const Fac1 eps_abs , const Fac1 eps_rel , const Fac1 a_x , const Fac1 a_dxdt )
- : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt ) { }
- template< class Tuple >
- __host__ __device__
- void operator()( Tuple t ) const
- {
- using std::abs;
- thrust::get< 0 >( t ) = abs( thrust::get< 0 >( t ) ) /
- ( m_eps_abs + m_eps_rel * ( m_a_x * abs( thrust::get< 1 >( t ) + m_a_dxdt * abs( thrust::get< 2 >( t ) ) ) ) );
- }
- typedef void result_type;
- };
- };
- }
- }
- }
- #endif
|