1 /* 2 [auto_generated] 3 boost/numeric/odeint/external/compute/compute_operations.hpp 4 5 [begin_description] 6 Operations of Boost.Compute zipped iterators. Is the counterpart of the compute_algebra. 7 [end_description] 8 9 Copyright 2009-2011 Karsten Ahnert 10 Copyright 2009-2011 Mario Mulansky 11 12 Distributed under the Boost Software License, Version 1.0. 13 (See accompanying file LICENSE_1_0.txt or 14 copy at http://www.boost.org/LICENSE_1_0.txt) 15 */ 16 17 18 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_COMPUTE_COMPUTE_OPERATIONS_HPP_DEFINED 19 #define BOOST_NUMERIC_ODEINT_EXTERNAL_COMPUTE_COMPUTE_OPERATIONS_HPP_DEFINED 20 21 #include <boost/preprocessor/repetition.hpp> 22 #include <boost/compute.hpp> 23 24 namespace boost { 25 namespace numeric { 26 namespace odeint { 27 28 struct compute_operations { 29 30 #define BOOST_ODEINT_COMPUTE_TEMPL_FAC(z, n, unused) \ 31 , class Fac ## n = BOOST_PP_CAT(Fac, BOOST_PP_DEC(n)) 32 33 #define BOOST_ODEINT_COMPUTE_MEMB_FAC(z, n, unused) \ 34 const Fac ## n m_alpha ## n; 35 36 #define BOOST_ODEINT_COMPUTE_PRM_FAC(z, n, unused) \ 37 BOOST_PP_COMMA_IF(n) const Fac ## n alpha ## n 38 39 #define BOOST_ODEINT_COMPUTE_INIT_FAC(z, n, unused) \ 40 BOOST_PP_COMMA_IF(n) m_alpha ## n (alpha ## n) 41 42 #define BOOST_ODEINT_COMPUTE_PRM_STATE(z, n, unused) \ 43 BOOST_PP_COMMA_IF(n) StateType ## n &s ## n 44 45 #define BOOST_ODEINT_COMPUTE_BEGIN_STATE(z, n, unused) \ 46 BOOST_PP_COMMA_IF( BOOST_PP_DEC(n) ) s ## n.begin() 47 48 #define BOOST_ODEINT_COMPUTE_END_STATE(z, n, unused) \ 49 BOOST_PP_COMMA_IF( BOOST_PP_DEC(n) ) s ## n.end() 50 51 #define BOOST_ODEINT_COMPUTE_LAMBDA(z, n, unused) \ 52 BOOST_PP_EXPR_IF(n, +) m_alpha ## n * bc::lambda::get< n >(bc::_1) 53 54 #define BOOST_ODEINT_COMPUTE_OPERATIONS(z, n, unused) \ 55 template< \ 56 class Fac0 = double \ 57 BOOST_PP_REPEAT_FROM_TO(1, n, BOOST_ODEINT_COMPUTE_TEMPL_FAC, ~) \ 58 > \ 59 struct scale_sum ## n { \ 60 BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_MEMB_FAC, ~) \ 61 scale_sum ## n( \ 62 BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_PRM_FAC, ~) \ 63 ) \ 64 : BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_INIT_FAC, ~) \ 65 { } \ 66 template< BOOST_PP_ENUM_PARAMS(BOOST_PP_INC(n), class StateType) > \ 67 void operator()( \ 68 BOOST_PP_REPEAT( \ 69 BOOST_PP_INC(n), \ 70 BOOST_ODEINT_COMPUTE_PRM_STATE, ~) \ 71 ) const \ 72 { \ 73 namespace bc = boost::compute; \ 74 bc::transform( \ 75 bc::make_zip_iterator( \ 76 boost::make_tuple( \ 77 BOOST_PP_REPEAT_FROM_TO( \ 78 1, BOOST_PP_INC(n), \ 79 BOOST_ODEINT_COMPUTE_BEGIN_STATE, ~) \ 80 ) \ 81 ), \ 82 bc::make_zip_iterator( \ 83 boost::make_tuple( \ 84 BOOST_PP_REPEAT_FROM_TO( \ 85 1, BOOST_PP_INC(n), \ 86 BOOST_ODEINT_COMPUTE_END_STATE, ~) \ 87 ) \ 88 ), \ 89 s0.begin(), \ 90 BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_LAMBDA, ~) \ 91 ); \ 92 } \ 93 }; 94 95 BOOST_PP_REPEAT_FROM_TO(2, 8, BOOST_ODEINT_COMPUTE_OPERATIONS, ~) 96 97 #undef BOOST_ODEINT_COMPUTE_TEMPL_FAC 98 #undef BOOST_ODEINT_COMPUTE_MEMB_FAC 99 #undef BOOST_ODEINT_COMPUTE_PRM_FAC 100 #undef BOOST_ODEINT_COMPUTE_INIT_FAC 101 #undef BOOST_ODEINT_COMPUTE_PRM_STATE 102 #undef BOOST_ODEINT_COMPUTE_BEGIN_STATE 103 #undef BOOST_ODEINT_COMPUTE_END_STATE 104 #undef BOOST_ODEINT_COMPUTE_LAMBDA 105 #undef BOOST_ODEINT_COMPUTE_OPERATIONS 106 107 template<class Fac1 = double, class Fac2 = Fac1> 108 struct scale_sum_swap2 { 109 const Fac1 m_alpha1; 110 const Fac2 m_alpha2; 111 scale_sum_swap2boost::numeric::odeint::compute_operations::scale_sum_swap2112 scale_sum_swap2(const Fac1 alpha1, const Fac2 alpha2) 113 : m_alpha1(alpha1), m_alpha2(alpha2) { } 114 115 template<class State0, class State1, class State2> operator ()boost::numeric::odeint::compute_operations::scale_sum_swap2116 void operator()(State0 &s0, State1 &s1, State2 &s2) const { 117 namespace bc = boost::compute; 118 119 bc::command_queue &queue = bc::system::default_queue(); 120 const bc::context &context = queue.get_context(); 121 122 const char source[] = BOOST_COMPUTE_STRINGIZE_SOURCE( 123 kernel void scale_sum_swap2( 124 F1 a1, F2 a2, 125 global T0 *x0, global T1 *x1, global T2 *x2, 126 ) 127 { 128 uint i = get_global_id(0); 129 T0 tmp = x0[i]; 130 x0[i] = a1 * x1[i] + a2 * x2[i]; 131 x1[i] = tmp; 132 } 133 ); 134 135 std::stringstream options; 136 options 137 << " -DT0=" << bc::type_name<typename State0::value_type>() 138 << " -DT1=" << bc::type_name<typename State1::value_type>() 139 << " -DT2=" << bc::type_name<typename State2::value_type>() 140 << " -DF1=" << bc::type_name<Fac1>() 141 << " -DF2=" << bc::type_name<Fac2>(); 142 143 bc::program program = 144 bc::program::build_with_source(source, context, options.str()); 145 146 bc::kernel kernel(program, "scale_sum_swap2"); 147 kernel.set_arg(0, m_alpha1); 148 kernel.set_arg(1, m_alpha2); 149 kernel.set_arg(2, s0.get_buffer()); 150 kernel.set_arg(3, s1.get_buffer()); 151 kernel.set_arg(4, s2.get_buffer()); 152 153 queue.enqueue_1d_range_kernel(kernel, 0, s0.size()); 154 155 } 156 }; 157 158 template<class Fac1 = double> 159 struct rel_error { 160 const Fac1 m_eps_abs, m_eps_rel, m_a_x, m_a_dxdt; 161 rel_errorboost::numeric::odeint::compute_operations::rel_error162 rel_error(const Fac1 eps_abs, const Fac1 eps_rel, const Fac1 a_x, const Fac1 a_dxdt) 163 : m_eps_abs(eps_abs), m_eps_rel(eps_rel), m_a_x(a_x), m_a_dxdt(a_dxdt) { } 164 165 166 template <class State0, class State1, class State2> operator ()boost::numeric::odeint::compute_operations::rel_error167 void operator()(State0 &s0, State1 &s1, State2 &s2) const { 168 namespace bc = boost::compute; 169 using bc::_1; 170 using bc::lambda::get; 171 172 bc::for_each( 173 bc::make_zip_iterator( 174 boost::make_tuple( 175 s0.begin(), 176 s1.begin(), 177 s2.begin() 178 ) 179 ), 180 bc::make_zip_iterator( 181 boost::make_tuple( 182 s0.end(), 183 s1.end(), 184 s2.end() 185 ) 186 ), 187 get<0>(_1) = abs( get<0>(_1) ) / 188 (m_eps_abs + m_eps_rel * (m_a_x * abs(get<1>(_1) + m_a_dxdt * abs(get<2>(_1))))) 189 ); 190 } 191 }; 192 }; 193 194 } // odeint 195 } // numeric 196 } // boost 197 198 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_COMPUTE_COMPUTE_OPERATIONS_HPP_DEFINED 199