File indexing completed on 2025-01-18 09:42:50
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_COMPUTE_COMPUTE_OPERATIONS_HPP_DEFINED
0019 #define BOOST_NUMERIC_ODEINT_EXTERNAL_COMPUTE_COMPUTE_OPERATIONS_HPP_DEFINED
0020
0021 #include <boost/preprocessor/repetition.hpp>
0022 #include <boost/compute.hpp>
0023
0024 namespace boost {
0025 namespace numeric {
0026 namespace odeint {
0027
0028 struct compute_operations {
0029
0030 #define BOOST_ODEINT_COMPUTE_TEMPL_FAC(z, n, unused) \
0031 , class Fac ## n = BOOST_PP_CAT(Fac, BOOST_PP_DEC(n))
0032
0033 #define BOOST_ODEINT_COMPUTE_MEMB_FAC(z, n, unused) \
0034 const Fac ## n m_alpha ## n;
0035
0036 #define BOOST_ODEINT_COMPUTE_PRM_FAC(z, n, unused) \
0037 BOOST_PP_COMMA_IF(n) const Fac ## n alpha ## n
0038
0039 #define BOOST_ODEINT_COMPUTE_INIT_FAC(z, n, unused) \
0040 BOOST_PP_COMMA_IF(n) m_alpha ## n (alpha ## n)
0041
0042 #define BOOST_ODEINT_COMPUTE_PRM_STATE(z, n, unused) \
0043 BOOST_PP_COMMA_IF(n) StateType ## n &s ## n
0044
0045 #define BOOST_ODEINT_COMPUTE_BEGIN_STATE(z, n, unused) \
0046 BOOST_PP_COMMA_IF( BOOST_PP_DEC(n) ) s ## n.begin()
0047
0048 #define BOOST_ODEINT_COMPUTE_END_STATE(z, n, unused) \
0049 BOOST_PP_COMMA_IF( BOOST_PP_DEC(n) ) s ## n.end()
0050
0051 #define BOOST_ODEINT_COMPUTE_LAMBDA(z, n, unused) \
0052 BOOST_PP_EXPR_IF(n, +) m_alpha ## n * bc::lambda::get< n >(bc::_1)
0053
0054 #define BOOST_ODEINT_COMPUTE_OPERATIONS(z, n, unused) \
0055 template< \
0056 class Fac0 = double \
0057 BOOST_PP_REPEAT_FROM_TO(1, n, BOOST_ODEINT_COMPUTE_TEMPL_FAC, ~) \
0058 > \
0059 struct scale_sum ## n { \
0060 BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_MEMB_FAC, ~) \
0061 scale_sum ## n( \
0062 BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_PRM_FAC, ~) \
0063 ) \
0064 : BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_INIT_FAC, ~) \
0065 { } \
0066 template< BOOST_PP_ENUM_PARAMS(BOOST_PP_INC(n), class StateType) > \
0067 void operator()( \
0068 BOOST_PP_REPEAT( \
0069 BOOST_PP_INC(n), \
0070 BOOST_ODEINT_COMPUTE_PRM_STATE, ~) \
0071 ) const \
0072 { \
0073 namespace bc = boost::compute; \
0074 bc::transform( \
0075 bc::make_zip_iterator( \
0076 boost::make_tuple( \
0077 BOOST_PP_REPEAT_FROM_TO( \
0078 1, BOOST_PP_INC(n), \
0079 BOOST_ODEINT_COMPUTE_BEGIN_STATE, ~) \
0080 ) \
0081 ), \
0082 bc::make_zip_iterator( \
0083 boost::make_tuple( \
0084 BOOST_PP_REPEAT_FROM_TO( \
0085 1, BOOST_PP_INC(n), \
0086 BOOST_ODEINT_COMPUTE_END_STATE, ~) \
0087 ) \
0088 ), \
0089 s0.begin(), \
0090 BOOST_PP_REPEAT(n, BOOST_ODEINT_COMPUTE_LAMBDA, ~) \
0091 ); \
0092 } \
0093 };
0094
0095 BOOST_PP_REPEAT_FROM_TO(2, 8, BOOST_ODEINT_COMPUTE_OPERATIONS, ~)
0096
0097 #undef BOOST_ODEINT_COMPUTE_TEMPL_FAC
0098 #undef BOOST_ODEINT_COMPUTE_MEMB_FAC
0099 #undef BOOST_ODEINT_COMPUTE_PRM_FAC
0100 #undef BOOST_ODEINT_COMPUTE_INIT_FAC
0101 #undef BOOST_ODEINT_COMPUTE_PRM_STATE
0102 #undef BOOST_ODEINT_COMPUTE_BEGIN_STATE
0103 #undef BOOST_ODEINT_COMPUTE_END_STATE
0104 #undef BOOST_ODEINT_COMPUTE_LAMBDA
0105 #undef BOOST_ODEINT_COMPUTE_OPERATIONS
0106
0107 template<class Fac1 = double, class Fac2 = Fac1>
0108 struct scale_sum_swap2 {
0109 const Fac1 m_alpha1;
0110 const Fac2 m_alpha2;
0111
0112 scale_sum_swap2(const Fac1 alpha1, const Fac2 alpha2)
0113 : m_alpha1(alpha1), m_alpha2(alpha2) { }
0114
0115 template<class State0, class State1, class State2>
0116 void operator()(State0 &s0, State1 &s1, State2 &s2) const {
0117 namespace bc = boost::compute;
0118
0119 bc::command_queue &queue = bc::system::default_queue();
0120 const bc::context &context = queue.get_context();
0121
0122 const char source[] = BOOST_COMPUTE_STRINGIZE_SOURCE(
0123 kernel void scale_sum_swap2(
0124 F1 a1, F2 a2,
0125 global T0 *x0, global T1 *x1, global T2 *x2,
0126 )
0127 {
0128 uint i = get_global_id(0);
0129 T0 tmp = x0[i];
0130 x0[i] = a1 * x1[i] + a2 * x2[i];
0131 x1[i] = tmp;
0132 }
0133 );
0134
0135 std::stringstream options;
0136 options
0137 << " -DT0=" << bc::type_name<typename State0::value_type>()
0138 << " -DT1=" << bc::type_name<typename State1::value_type>()
0139 << " -DT2=" << bc::type_name<typename State2::value_type>()
0140 << " -DF1=" << bc::type_name<Fac1>()
0141 << " -DF2=" << bc::type_name<Fac2>();
0142
0143 bc::program program =
0144 bc::program::build_with_source(source, context, options.str());
0145
0146 bc::kernel kernel(program, "scale_sum_swap2");
0147 kernel.set_arg(0, m_alpha1);
0148 kernel.set_arg(1, m_alpha2);
0149 kernel.set_arg(2, s0.get_buffer());
0150 kernel.set_arg(3, s1.get_buffer());
0151 kernel.set_arg(4, s2.get_buffer());
0152
0153 queue.enqueue_1d_range_kernel(kernel, 0, s0.size());
0154
0155 }
0156 };
0157
0158 template<class Fac1 = double>
0159 struct rel_error {
0160 const Fac1 m_eps_abs, m_eps_rel, m_a_x, m_a_dxdt;
0161
0162 rel_error(const Fac1 eps_abs, const Fac1 eps_rel, const Fac1 a_x, const Fac1 a_dxdt)
0163 : m_eps_abs(eps_abs), m_eps_rel(eps_rel), m_a_x(a_x), m_a_dxdt(a_dxdt) { }
0164
0165
0166 template <class State0, class State1, class State2>
0167 void operator()(State0 &s0, State1 &s1, State2 &s2) const {
0168 namespace bc = boost::compute;
0169 using bc::_1;
0170 using bc::lambda::get;
0171
0172 bc::for_each(
0173 bc::make_zip_iterator(
0174 boost::make_tuple(
0175 s0.begin(),
0176 s1.begin(),
0177 s2.begin()
0178 )
0179 ),
0180 bc::make_zip_iterator(
0181 boost::make_tuple(
0182 s0.end(),
0183 s1.end(),
0184 s2.end()
0185 )
0186 ),
0187 get<0>(_1) = abs( get<0>(_1) ) /
0188 (m_eps_abs + m_eps_rel * (m_a_x * abs(get<1>(_1) + m_a_dxdt * abs(get<2>(_1)))))
0189 );
0190 }
0191 };
0192 };
0193
0194 }
0195 }
0196 }
0197
0198 #endif
