compute_operations.hpp 8.5 KB

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