diff options
Diffstat (limited to 'boost/numeric')
45 files changed, 1100 insertions, 338 deletions
diff --git a/boost/numeric/interval/compare/set.hpp b/boost/numeric/interval/compare/set.hpp index aa4f1716b..44f4d6ee6 100644 --- a/boost/numeric/interval/compare/set.hpp +++ b/boost/numeric/interval/compare/set.hpp @@ -27,7 +27,7 @@ bool operator<(const interval<T, Policies1>& x, const interval<T, Policies2>& y) } template<class T, class Policies> inline -bool operator<(const interval<T, Policies>& x, const T& y) +bool operator<(const interval<T, Policies>& , const T& ) { throw comparison_error(); } @@ -39,7 +39,7 @@ bool operator<=(const interval<T, Policies1>& x, const interval<T, Policies2>& y } template<class T, class Policies> inline -bool operator<=(const interval<T, Policies>& x, const T& y) +bool operator<=(const interval<T, Policies>& , const T& ) { throw comparison_error(); } @@ -51,7 +51,7 @@ bool operator>(const interval<T, Policies1>& x, const interval<T, Policies2>& y) } template<class T, class Policies> inline -bool operator>(const interval<T, Policies>& x, const T& y) +bool operator>(const interval<T, Policies>& , const T& ) { throw comparison_error(); } @@ -63,7 +63,7 @@ bool operator>=(const interval<T, Policies1>& x, const interval<T, Policies2>& y } template<class T, class Policies> inline -bool operator>=(const interval<T, Policies>& x, const T& y) +bool operator>=(const interval<T, Policies>& , const T& ) { throw comparison_error(); } @@ -75,7 +75,7 @@ bool operator==(const interval<T, Policies1>& x, const interval<T, Policies2>& y } template<class T, class Policies> inline -bool operator==(const interval<T, Policies>& x, const T& y) +bool operator==(const interval<T, Policies>& , const T& ) { throw comparison_error(); } @@ -87,7 +87,7 @@ bool operator!=(const interval<T, Policies1>& x, const interval<T, Policies2>& y } template<class T, class Policies> inline -bool operator!=(const interval<T, Policies>& x, const T& y) +bool operator!=(const interval<T, Policies>& , const T& ) { throw comparison_error(); } diff --git a/boost/numeric/interval/detail/bugs.hpp b/boost/numeric/interval/detail/bugs.hpp index cc37988d7..b9b421a23 100644 --- a/boost/numeric/interval/detail/bugs.hpp +++ b/boost/numeric/interval/detail/bugs.hpp @@ -45,35 +45,4 @@ # define BOOST_NUMERIC_INTERVAL_using_ahyp(a) #endif -#if defined(__GNUC__) && (__GNUC__ <= 2) -// cf PR c++/1981 for a description of the bug -#include <algorithm> -#include <boost/config/no_tr1/cmath.hpp> -namespace boost { -namespace numeric { - using std::min; - using std::max; - using std::sqrt; - using std::exp; - using std::log; - using std::cos; - using std::tan; - using std::asin; - using std::acos; - using std::atan; - using std::ceil; - using std::floor; - using std::sinh; - using std::cosh; - using std::tanh; -# undef BOOST_NUMERIC_INTERVAL_using_max -# undef BOOST_NUMERIC_INTERVAL_using_math -# define BOOST_NUMERIC_INTERVAL_using_max(a) -# define BOOST_NUMERIC_INTERVAL_using_math(a) -# undef BOOST_NUMERIC_INTERVAL_using_ahyp -# define BOOST_NUMERIC_INTERVAL_using_ahyp(a) -} // namespace numeric -} // namespace boost -#endif - #endif // BOOST_NUMERIC_INTERVAL_DETAIL_BUGS diff --git a/boost/numeric/interval/detail/msvc_rounding_control.hpp b/boost/numeric/interval/detail/msvc_rounding_control.hpp index d36e80242..ccf1e8b47 100644 --- a/boost/numeric/interval/detail/msvc_rounding_control.hpp +++ b/boost/numeric/interval/detail/msvc_rounding_control.hpp @@ -88,7 +88,20 @@ struct x86_rounding static void get_rounding_mode(rounding_mode& mode) { mode = msvc2hard(_control87(0, 0)); } static void set_rounding_mode(const rounding_mode mode) - { _control87(hard2msvc(mode), _MCW_EM | _MCW_RC | _MCW_PC | _MCW_IC); } + { + _control87(hard2msvc(mode), + _MCW_EM | _MCW_RC +#if !defined(_M_AMD64) && !defined(_M_ARM) + // x64 ignores _MCW_PC and _MCW_IC, and the Debug CRT library actually + // asserts when these are passed to _control87. + // MSDN says on '_control87' that changing precision (_MCW_PC) or + // infinity (_MCW_IC) handling is not supported on the ARM and x64 + // architectures and that _control87 raises an assertion + // and the invalid parameter handler is invoked. + | _MCW_PC | _MCW_IC +#endif + ); + } static double to_int(const double& x) { return rint(x); } }; diff --git a/boost/numeric/interval/limits.hpp b/boost/numeric/interval/limits.hpp index d691ccee2..973681dbc 100644 --- a/boost/numeric/interval/limits.hpp +++ b/boost/numeric/interval/limits.hpp @@ -11,7 +11,6 @@ #ifndef BOOST_NUMERIC_INTERVAL_LIMITS_HPP #define BOOST_NUMERIC_INTERVAL_LIMITS_HPP -#ifndef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION #include <boost/config.hpp> #include <boost/limits.hpp> @@ -46,6 +45,5 @@ private: } // namespace std -#endif #endif // BOOST_NUMERIC_INTERVAL_LIMITS_HPP diff --git a/boost/numeric/interval/rounded_arith.hpp b/boost/numeric/interval/rounded_arith.hpp index 8b2d9a718..2520274a6 100644 --- a/boost/numeric/interval/rounded_arith.hpp +++ b/boost/numeric/interval/rounded_arith.hpp @@ -68,8 +68,8 @@ struct rounded_arith_std: Rounding { { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_DN(sqrt(x)); } T sqrt_up (const T& x) { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_UP(sqrt(x)); } - T int_down(const T& x) { this->downward(); return to_int(x); } - T int_up (const T& x) { this->upward(); return to_int(x); } + T int_down(const T& x) { this->downward(); return this->to_int(x); } + T int_up (const T& x) { this->upward(); return this->to_int(x); } # undef BOOST_DN # undef BOOST_NR # undef BOOST_UP @@ -105,8 +105,8 @@ struct rounded_arith_opp: Rounding { { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_DN(sqrt(x)); } T sqrt_up (const T& x) { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_UP(sqrt(x)); } - T int_down(const T& x) { return -to_int(-x); } - T int_up (const T& x) { return to_int(x); } + T int_down(const T& x) { return -this->to_int(-x); } + T int_up (const T& x) { return this->to_int(x); } # undef BOOST_DN # undef BOOST_NR # undef BOOST_UP diff --git a/boost/numeric/odeint/algebra/algebra_dispatcher.hpp b/boost/numeric/odeint/algebra/algebra_dispatcher.hpp index df5aabac9..88cf159e0 100644 --- a/boost/numeric/odeint/algebra/algebra_dispatcher.hpp +++ b/boost/numeric/odeint/algebra/algebra_dispatcher.hpp @@ -17,8 +17,9 @@ #ifndef BOOST_NUMERIC_ODEINT_ALGEBRA_ALGEBRA_DISPATCHER_HPP_INCLUDED #define BOOST_NUMERIC_ODEINT_ALGEBRA_ALGEBRA_DISPATCHER_HPP_INCLUDED -#include <complex> +#include <boost/numeric/odeint/config.hpp> +#include <complex> #include <boost/type_traits/is_floating_point.hpp> #include <boost/numeric/ublas/vector.hpp> @@ -38,14 +39,14 @@ namespace odeint { template< class StateType , class Enabler = void > struct algebra_dispatcher_sfinae { - // range_algebra is the standard algebra^ + // range_algebra is the standard algebra typedef range_algebra algebra_type; }; template< class StateType > struct algebra_dispatcher : algebra_dispatcher_sfinae< StateType > { }; -//specialize for array +// specialize for array template< class T , size_t N > struct algebra_dispatcher< boost::array< T , N > > { @@ -84,4 +85,26 @@ struct algebra_dispatcher< boost::numeric::ublas::matrix< T , L , A > > } } +#ifdef BOOST_NUMERIC_ODEINT_CXX11 + +// c++11 mode: specialization for std::array if available + +#include <array> + +namespace boost { +namespace numeric { +namespace odeint { + +// specialize for std::array +template< class T , size_t N > +struct algebra_dispatcher< std::array< T , N > > +{ + typedef array_algebra algebra_type; +}; + +} } } + +#endif + + #endif diff --git a/boost/numeric/odeint/algebra/array_algebra.hpp b/boost/numeric/odeint/algebra/array_algebra.hpp index c565d6dd7..471e866fe 100644 --- a/boost/numeric/odeint/algebra/array_algebra.hpp +++ b/boost/numeric/odeint/algebra/array_algebra.hpp @@ -3,7 +3,11 @@ boost/numeric/odeint/algebra/array_algebra.hpp [begin_description] - Algebra for boost::array. Highly specialized for odeint. Const arguments are introduce to work with odeint. + Algebra for Arrays. Highly specialized for odeint. Const arguments are + introduce to work with odeint. + The Array algebra can be used for Array structures with two template + parameters: + Array<T, N> [end_description] Copyright 2011-2013 Mario Mulansky @@ -29,234 +33,251 @@ namespace odeint { struct array_algebra { - template< typename T , size_t dim , class Op > - static void for_each1( boost::array< T , dim > &s1 , Op op ) + //template< typename T , size_t dim , class Op > + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each1( Array< T, dim > &s1, Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] ); } - - template< typename T1 , typename T2 , size_t dim , class Op > - static void for_each2( boost::array< T1 , dim > &s1 , - const boost::array< T2 , dim > &s2 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each2( Array< T, dim > &s1, const Array< T, dim > &s2, + Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] ); } - template< typename T , size_t dim , class Op > - static void for_each3( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each3( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] ); } /* different const signature - required for the scale_sum_swap2 operation */ - template< typename T , size_t dim , class Op > - static void for_each3( boost::array< T , dim > &s1 , - boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each3( Array< T , dim > &s1 , + Array< T , dim > &s2 , + const Array< T , dim > &s3 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] ); } - template< typename T , size_t dim , class Op > - static void for_each4( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each4( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] ); } - template< typename T , size_t dim , class Op > - static void for_each5( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each5( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] ); } - template< typename T , size_t dim , class Op > - static void for_each6( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each6( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] ); } - template< typename T , size_t dim , class Op > - static void for_each7( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each7( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] ); } - template< typename T , size_t dim , class Op > - static void for_each8( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each8( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] ); } - template< typename T , size_t dim , class Op > - static void for_each9( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , - const boost::array< T , dim > &s9 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each9( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , + const Array< T , dim > &s9 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] ); } - template< typename T , size_t dim , class Op > - static void for_each10( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , - const boost::array< T , dim > &s9 , - const boost::array< T , dim > &s10 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each10( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , + const Array< T , dim > &s9 , + const Array< T , dim > &s10 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] ); } - template< typename T , size_t dim , class Op > - static void for_each11( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , - const boost::array< T , dim > &s9 , - const boost::array< T , dim > &s10 , - const boost::array< T , dim > &s11 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each11( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , + const Array< T , dim > &s9 , + const Array< T , dim > &s10 , + const Array< T , dim > &s11 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] ); } - template< typename T , size_t dim , class Op > - static void for_each12( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , - const boost::array< T , dim > &s9 , - const boost::array< T , dim > &s10 , - const boost::array< T , dim > &s11 , - const boost::array< T , dim > &s12 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each12( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , + const Array< T , dim > &s9 , + const Array< T , dim > &s10 , + const Array< T , dim > &s11 , + const Array< T , dim > &s12 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] ); } - template< typename T , size_t dim , class Op > - static void for_each13( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , - const boost::array< T , dim > &s9 , - const boost::array< T , dim > &s10 , - const boost::array< T , dim > &s11 , - const boost::array< T , dim > &s12 , - const boost::array< T , dim > &s13 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each13( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , + const Array< T , dim > &s9 , + const Array< T , dim > &s10 , + const Array< T , dim > &s11 , + const Array< T , dim > &s12 , + const Array< T , dim > &s13 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] ); } - template< typename T , size_t dim , class Op > - static void for_each14( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , - const boost::array< T , dim > &s9 , - const boost::array< T , dim > &s10 , - const boost::array< T , dim > &s11 , - const boost::array< T , dim > &s12 , - const boost::array< T , dim > &s13 , - const boost::array< T , dim > &s14 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each14( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , + const Array< T , dim > &s9 , + const Array< T , dim > &s10 , + const Array< T , dim > &s11 , + const Array< T , dim > &s12 , + const Array< T , dim > &s13 , + const Array< T , dim > &s14 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] , s14[i] ); } - template< typename T , size_t dim , class Op > - static void for_each15( boost::array< T , dim > &s1 , - const boost::array< T , dim > &s2 , - const boost::array< T , dim > &s3 , - const boost::array< T , dim > &s4 , - const boost::array< T , dim > &s5 , - const boost::array< T , dim > &s6 , - const boost::array< T , dim > &s7 , - const boost::array< T , dim > &s8 , - const boost::array< T , dim > &s9 , - const boost::array< T , dim > &s10 , - const boost::array< T , dim > &s11 , - const boost::array< T , dim > &s12 , - const boost::array< T , dim > &s13 , - const boost::array< T , dim > &s14 , - const boost::array< T , dim > &s15 , Op op ) + template < template < typename, size_t > class Array, typename T, + size_t dim, class Op > + static void for_each15( Array< T , dim > &s1 , + const Array< T , dim > &s2 , + const Array< T , dim > &s3 , + const Array< T , dim > &s4 , + const Array< T , dim > &s5 , + const Array< T , dim > &s6 , + const Array< T , dim > &s7 , + const Array< T , dim > &s8 , + const Array< T , dim > &s9 , + const Array< T , dim > &s10 , + const Array< T , dim > &s11 , + const Array< T , dim > &s12 , + const Array< T , dim > &s13 , + const Array< T , dim > &s14 , + const Array< T , dim > &s15 , Op op ) { for( size_t i=0 ; i<dim ; ++i ) op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] , s14[i] , s15[i] ); } - template< typename T , size_t dim > - static typename norm_result_type< boost::array< T , dim > >::type norm_inf( const boost::array< T , dim > &s ) + template < template < typename, size_t > class Array, typename T, + size_t dim> + static typename norm_result_type< Array< T , dim > >::type norm_inf( const Array< T , dim > &s ) { BOOST_USING_STD_MAX(); using std::abs; - typedef typename norm_result_type< boost::array< T , dim > >::type result_type; + typedef typename norm_result_type< Array< T , dim > >::type result_type; result_type init = static_cast< result_type >( 0 ); for( size_t i=0 ; i<dim ; ++i ) init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type >(abs(s[i])) ); diff --git a/boost/numeric/odeint/algebra/multi_array_algebra.hpp b/boost/numeric/odeint/algebra/multi_array_algebra.hpp index 4d1fc16c3..0bc476e63 100644 --- a/boost/numeric/odeint/algebra/multi_array_algebra.hpp +++ b/boost/numeric/odeint/algebra/multi_array_algebra.hpp @@ -113,13 +113,13 @@ struct multi_array_algebra template< class S1 , class S2 , class S3 , class S4 , class S5 , class S6 ,class S7 , class S8 , class S9 , class S10 , class S11 , class S12 , class S13 , class S14 , class Op > static void for_each14( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , S5 &s5 , S6 &s6 , S7 &s7 , S8 &s8 , S9 &s9 , S10 &s10 , S11 &s11 , S12 &s12 , S13 &s13 , S14 &s14 , Op op ) { - detail::for_each14( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s14.data() , op ); + detail::for_each14( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , op ); } template< class S1 , class S2 , class S3 , class S4 , class S5 , class S6 ,class S7 , class S8 , class S9 , class S10 , class S11 , class S12 , class S13 , class S14 , class S15 , class Op > static void for_each15( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , S5 &s5 , S6 &s6 , S7 &s7 , S8 &s8 , S9 &s9 , S10 &s10 , S11 &s11 , S12 &s12 , S13 &s13 , S14 &s14 , S15 &s15 , Op op ) { - detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s14.data() , s15.data() , op ); + detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , s15.data() , op ); } template< typename S > diff --git a/boost/numeric/odeint/external/eigen/eigen.hpp b/boost/numeric/odeint/external/eigen/eigen.hpp new file mode 100644 index 000000000..620acae51 --- /dev/null +++ b/boost/numeric/odeint/external/eigen/eigen.hpp @@ -0,0 +1,27 @@ +/* + [auto_generated] + boost/numeric/odeint/external/eigen/eigen.hpp + + [begin_description] + tba. + [end_description] + + Copyright 2009-2012 Karsten Ahnert + Copyright 2009-2012 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) +*/ + + +#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED + + +#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp> +#include <boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp> +#include <boost/numeric/odeint/external/eigen/eigen_resize.hpp> + + +#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED diff --git a/boost/numeric/odeint/external/eigen/eigen_algebra.hpp b/boost/numeric/odeint/external/eigen/eigen_algebra.hpp index 6ca2026b9..b4ee5c3b1 100644 --- a/boost/numeric/odeint/external/eigen/eigen_algebra.hpp +++ b/boost/numeric/odeint/external/eigen/eigen_algebra.hpp @@ -71,6 +71,19 @@ operator/(const Eigen::MatrixBase<D1> &x1, const Eigen::MatrixBase<D2> &x2) { return x1.cwiseQuotient(x2); } + +template< typename D > +inline const +typename Eigen::CwiseUnaryOp< + typename Eigen::internal::scalar_abs_op< + typename Eigen::internal::traits< D >::Scalar > , + const D > +abs( const Eigen::MatrixBase< D > &m ) { + return m.cwiseAbs(); +} + + + } // end Eigen namespace diff --git a/boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp b/boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp new file mode 100644 index 000000000..6c8a3a276 --- /dev/null +++ b/boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp @@ -0,0 +1,49 @@ +/* + [auto_generated] + boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp + + [begin_description] + tba. + [end_description] + + Copyright 2009-2012 Karsten Ahnert + Copyright 2009-2012 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) +*/ + + +#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED + + +namespace boost { +namespace numeric { +namespace odeint { + + +template< class Derived > +struct algebra_dispatcher_sfinae< Derived , + typename boost::enable_if< typename boost::is_base_of< Eigen::MatrixBase< Derived > , Derived >::type >::type > +{ + typedef vector_space_algebra algebra_type; +}; + + +template < class Derived > +struct algebra_dispatcher_sfinae< Derived , + typename boost::enable_if< typename boost::is_base_of< Eigen::ArrayBase< Derived > , Derived >::type >::type > +{ + typedef vector_space_algebra algebra_type; +}; + + + +} // namespace odeint +} // namespace numeric +} // namespace boost + + +#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED diff --git a/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp b/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp index 5deba2cb5..096134635 100644 --- a/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp +++ b/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp @@ -24,6 +24,7 @@ #include <boost/numeric/odeint/external/thrust/thrust_algebra.hpp> #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> +// specializations for the standard thrust containers namespace boost { namespace numeric { @@ -48,5 +49,59 @@ struct algebra_dispatcher< thrust::device_vector< T , A > > } // namespace boost +// add support for thrust backend vectors, if available + +#include <thrust/version.h> + +#if THRUST_VERSION >= 100600 + +// specialization for thrust cpp vector +#include <thrust/system/cpp/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct algebra_dispatcher< thrust::cpp::vector< T , A > > + { + typedef thrust_algebra algebra_type; + }; +} } } + +// specialization for thrust omp vector +#ifdef _OPENMP +#include <thrust/system/omp/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct algebra_dispatcher< thrust::omp::vector< T , A > > + { + typedef thrust_algebra algebra_type; + }; +} } } +#endif // _OPENMP + +// specialization for thrust tbb vector +#ifdef TBB_VERSION_MAJOR +#include <thrust/system/tbb/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct algebra_dispatcher< thrust::tbb::vector< T , A > > + { + typedef thrust_algebra algebra_type; + }; +} } } +#endif // TBB_VERSION_MAJOR + +// specialization for thrust cuda vector +#ifdef __CUDACC__ +#include <thrust/system/cuda/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct algebra_dispatcher< thrust::cuda::vector< T , A > > + { + typedef thrust_algebra algebra_type; + }; +} } } +#endif // __CUDACC__ + +#endif // THRUST_VERSION >= 100600 + #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_ALGEBRA_DISPATCHER_HPP_DEFINED diff --git a/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp b/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp index 0105a4e34..e9b3a6423 100644 --- a/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp +++ b/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp @@ -6,8 +6,8 @@ operations_dispatcher specialization for thrust [end_description] - Copyright 2013 Karsten Ahnert - Copyright 2013 Mario Mulansky + Copyright 2013-2014 Karsten Ahnert + Copyright 2013-2014 Mario Mulansky Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or @@ -24,6 +24,7 @@ #include <boost/numeric/odeint/external/thrust/thrust_operations.hpp> #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> +// support for the standard thrust containers namespace boost { namespace numeric { @@ -47,6 +48,60 @@ struct operations_dispatcher< thrust::device_vector< T , A > > } // namespace numeric } // namespace boost +// add support for thrust backend vectors, if available -#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_ALGEBRA_DISPATCHER_HPP_DEFINED +#include <thrust/version.h> + +#if THRUST_VERSION >= 100600 + +// specialization for thrust cpp vector +#include <thrust/system/cpp/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct operations_dispatcher< thrust::cpp::vector< T , A > > + { + typedef thrust_operations operations_type; + }; +} } } + +// specialization for thrust omp vector +#ifdef _OPENMP +#include <thrust/system/omp/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct operations_dispatcher< thrust::omp::vector< T , A > > + { + typedef thrust_operations operations_type; + }; +} } } +#endif // _OPENMP + +// specialization for thrust tbb vector +#ifdef TBB_VERSION_MAJOR +#include <thrust/system/tbb/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct operations_dispatcher< thrust::tbb::vector< T , A > > + { + typedef thrust_operations operations_type; + }; +} } } +#endif // TBB_VERSION_MAJOR + +// specialization for thrust cuda vector +#ifdef __CUDACC__ +#include <thrust/system/cuda/vector.h> +namespace boost { namespace numeric { namespace odeint { + template< class T , class A > + struct operations_dispatcher< thrust::cuda::vector< T , A > > + { + typedef thrust_operations operations_type; + }; +} } } +#endif // __CUDACC__ + +#endif // THRUST_VERSION >= 100600 + + +#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_OPERATIONS_DISPATCHER_HPP_DEFINED diff --git a/boost/numeric/odeint/external/thrust/thrust_resize.hpp b/boost/numeric/odeint/external/thrust/thrust_resize.hpp index 8f1b181e2..2f7f7b23b 100644 --- a/boost/numeric/odeint/external/thrust/thrust_resize.hpp +++ b/boost/numeric/odeint/external/thrust/thrust_resize.hpp @@ -6,7 +6,7 @@ Enable resizing for thrusts device and host_vector. [end_description] - Copyright 2010-2012 Mario Mulansky + Copyright 2010-2014 Mario Mulansky Copyright 2010-2011 Karsten Ahnert Distributed under the Boost Software License, Version 1.0. @@ -18,102 +18,168 @@ #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED #define BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED +#include <boost/range.hpp> #include <thrust/device_vector.h> #include <thrust/host_vector.h> +#include <thrust/distance.h> +#include <boost/numeric/odeint/util/resize.hpp> +#include <boost/numeric/odeint/util/same_size.hpp> #include <boost/numeric/odeint/util/copy.hpp> namespace boost { namespace numeric { namespace odeint { -template< class T > -struct is_resizeable< thrust::device_vector< T > > -{ - struct type : public boost::true_type { }; - const static bool value = type::value; -}; - -template< class T > -struct same_size_impl< thrust::device_vector< T > , thrust::device_vector< T > > -{ - static bool same_size( const thrust::device_vector< T > &x , const thrust::device_vector< T > &y ) - { - return x.size() == y.size(); - } -}; - -template< class T > -struct resize_impl< thrust::device_vector< T > , thrust::device_vector< T > > -{ - static void resize( thrust::device_vector< T > &x , const thrust::device_vector< T > &y ) - { - x.resize( y.size() ); - } -}; - - -template< class T > -struct is_resizeable< thrust::host_vector< T > > -{ - struct type : public boost::true_type { }; - const static bool value = type::value; -}; - -template< class T > -struct same_size_impl< thrust::host_vector< T > , thrust::host_vector< T > > -{ - static bool same_size( const thrust::host_vector< T > &x , const thrust::host_vector< T > &y ) - { - return x.size() == y.size(); - } -}; - -template< class T > -struct resize_impl< thrust::host_vector< T > , thrust::host_vector< T > > -{ - static void resize( thrust::host_vector< T > &x , const thrust::host_vector< T > &y ) - { - x.resize( y.size() ); - } -}; - - - -template< class Container1, class Value > -struct copy_impl< Container1 , thrust::device_vector< Value > > -{ - static void copy( const Container1 &from , thrust::device_vector< Value > &to ) - { - thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) ); - } -}; - -template< class Value , class Container2 > -struct copy_impl< thrust::device_vector< Value > , Container2 > -{ - static void copy( const thrust::device_vector< Value > &from , Container2 &to ) - { - thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) ); - } -}; - -template< class Value > -struct copy_impl< thrust::device_vector< Value > , thrust::device_vector< Value > > -{ - static void copy( const thrust::device_vector< Value > &from , thrust::device_vector< Value > &to ) - { - thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) ); - } -}; - - +// some macros that define the necessary utilities + +#define ODEINT_THRUST_VECTOR_IS_RESIZEABLE( THRUST_VECTOR ) \ +template< class T , class A > \ +struct is_resizeable< THRUST_VECTOR<T,A> > \ +{ \ + struct type : public boost::true_type { }; \ + const static bool value = type::value; \ +}; \ + +#define ODEINT_TRHUST_VECTOR_RESIZE_IMPL( THRUST_VECTOR ) \ +template< class T, class A > \ +struct resize_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> > \ +{ \ + static void resize( THRUST_VECTOR<T,A> &x , \ + const THRUST_VECTOR<T,A> &y ) \ + { \ + x.resize( y.size() ); \ + } \ +}; \ +template< class T, class A, typename Range > \ +struct resize_impl< THRUST_VECTOR<T,A> , Range > \ +{ \ + static void resize( THRUST_VECTOR<T,A> &x , \ + const Range &y ) \ + { \ + x.resize( thrust::distance(boost::begin(y), \ + boost::end(y))); \ + } \ +}; \ + + +#define ODEINT_THRUST_SAME_SIZE_IMPL( THRUST_VECTOR ) \ +template< class T , class A > \ +struct same_size_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> > \ +{ \ + static bool same_size( const THRUST_VECTOR<T,A> &x , \ + const THRUST_VECTOR<T,A> &y ) \ + { \ + return x.size() == y.size(); \ + } \ +}; \ +template< class T , class A, typename Range > \ +struct same_size_impl< THRUST_VECTOR<T,A> , Range > \ +{ \ + static bool same_size( const THRUST_VECTOR<T,A> &x , \ + const Range &y ) \ + { \ + return x.size() == thrust::distance(boost::begin(y), \ + boost::end(y)); \ + } \ +}; \ + + +#define ODEINT_THRUST_COPY_IMPL( THRUST_VECTOR ) \ +template< class Container1 , class T , class A > \ +struct copy_impl< Container1 , THRUST_VECTOR<T,A> > \ +{ \ + static void copy( const Container1 &from , THRUST_VECTOR<T,A> &to ) \ + { \ + thrust::copy( boost::begin( from ) , boost::end( from ) , \ + boost::begin( to ) ); \ + } \ +}; \ + \ +template< class T , class A , class Container2 > \ +struct copy_impl< THRUST_VECTOR<T,A> , Container2 > \ +{ \ + static void copy( const THRUST_VECTOR<T,A> &from , Container2 &to ) \ + { \ + thrust::copy( boost::begin( from ) , boost::end( from ) , \ + boost::begin( to ) ); \ + } \ +}; \ + \ +template< class T , class A > \ +struct copy_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> > \ +{ \ + static void copy( const THRUST_VECTOR<T,A> &from , \ + THRUST_VECTOR<T,A> &to ) \ + { \ + thrust::copy( boost::begin( from ) , boost::end( from ) , \ + boost::begin( to ) ); \ + } \ +}; \ + +// add support for the standard thrust containers + +ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::device_vector ) +ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::device_vector ) +ODEINT_THRUST_SAME_SIZE_IMPL( thrust::device_vector ) +ODEINT_THRUST_COPY_IMPL( thrust::device_vector ) + + +ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::host_vector ) +ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::host_vector ) +ODEINT_THRUST_SAME_SIZE_IMPL( thrust::host_vector ) +ODEINT_THRUST_COPY_IMPL( thrust::host_vector ) } // odeint } // numeric } // boost +// add support for thrust backend vectors, if available + +#include <thrust/version.h> + +#if THRUST_VERSION >= 100600 + +#include <thrust/system/cpp/vector.h> +namespace boost { namespace numeric { namespace odeint { + ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::cpp::vector ) + ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::cpp::vector ) + ODEINT_THRUST_SAME_SIZE_IMPL( thrust::cpp::vector ) + ODEINT_THRUST_COPY_IMPL( thrust::cpp::vector ) +} } } + +#ifdef _OPENMP +#include <thrust/system/omp/vector.h> +namespace boost { namespace numeric { namespace odeint { + ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::omp::vector ) + ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::omp::vector ) + ODEINT_THRUST_SAME_SIZE_IMPL( thrust::omp::vector ) + ODEINT_THRUST_COPY_IMPL( thrust::omp::vector ) +} } } +#endif // _OPENMP + +#ifdef TBB_VERSION_MAJOR +#include <thrust/system/tbb/vector.h> +namespace boost { namespace numeric { namespace odeint { + ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::tbb::vector ) + ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::tbb::vector ) + ODEINT_THRUST_SAME_SIZE_IMPL( thrust::tbb::vector ) + ODEINT_THRUST_COPY_IMPL( thrust::tbb::vector ) +} } } +#endif // TBB_VERSION_MAJOR + +#ifdef __CUDACC__ +#include <thrust/system/cuda/vector.h> +namespace boost { namespace numeric { namespace odeint { + ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::cuda::vector ) + ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::cuda::vector ) + ODEINT_THRUST_SAME_SIZE_IMPL( thrust::cuda::vector ) + ODEINT_THRUST_COPY_IMPL( thrust::cuda::vector ) +} } } +#endif // __CUDACC__ + +#endif // THRUST_VERSION >= 100600 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED diff --git a/boost/numeric/odeint/integrate/integrate.hpp b/boost/numeric/odeint/integrate/integrate.hpp index 6cb96c0f1..446656b58 100644 --- a/boost/numeric/odeint/integrate/integrate.hpp +++ b/boost/numeric/odeint/integrate/integrate.hpp @@ -48,6 +48,15 @@ integrate( System system , State &start_state , Time start_time , Time end_time return integrate_adaptive( stepper_type() , system , start_state , start_time , end_time , dt , observer ); } +template< class Value , class System , class State , class Time , class Observer > +size_t +integrate( System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer ) +{ + typedef controlled_runge_kutta< runge_kutta_dopri5< State , Value , State , Time > > stepper_type; + return integrate_adaptive( stepper_type() , system , start_state , start_time , end_time , dt , observer ); +} + + /* @@ -59,6 +68,13 @@ size_t integrate( System system , State &start_state , Time start_time , Time en return integrate( system , start_state , start_time , end_time , dt , null_observer() ); } +template< class Value , class System , class State , class Time > +size_t integrate( System system , State &start_state , Time start_time , Time end_time , Time dt ) +{ + return integrate< Value >( system , start_state , start_time , end_time , dt , null_observer() ); +} + + /** * \fn integrate( System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer ) @@ -70,6 +86,9 @@ size_t integrate( System system , State &start_state , Time start_time , Time en * integration with step size control, thus dt changes during the integration. * This method uses standard error bounds of 1E-6. * After each step, the observer is called. + * + * \attention A second version of this function template exists which explicitly + * expects the value type as template parameter, i.e. integrate< double >( sys , x , t0 , t1 , dt , obs ); * * \param system The system function to solve, hence the r.h.s. of the * ordinary differential equation. @@ -92,6 +111,9 @@ size_t integrate( System system , State &start_state , Time start_time , Time en * integration with step size control, thus dt changes during the integration. * This method uses standard error bounds of 1E-6. * No observer is called. + * + * \attention A second version of this function template exists which explicitly + * expects the value type as template parameter, i.e. integrate< double >( sys , x , t0 , t1 , dt ); * * \param system The system function to solve, hence the r.h.s. of the * ordinary differential equation. diff --git a/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp b/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp index 643b53937..e5b4b3cf4 100644 --- a/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp +++ b/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp @@ -205,7 +205,7 @@ namespace odeint { * \param sys The system function (ODE) to solve. * \param s The initial state. */ - adaptive_iterator_impl( stepper_type stepper , system_type sys , state_type &s ) + adaptive_iterator_impl( stepper_type stepper , system_type sys , state_type& /* s */ ) : base_type( stepper , sys ) { } protected: diff --git a/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp b/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp index 449acfa79..e23474c7a 100644 --- a/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp +++ b/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp @@ -77,7 +77,7 @@ namespace odeint { * \param sys The system function (ODE) to solve. * \param s The initial state. const_step_iterator stores a reference of s and changes its value during the iteration. */ - const_step_iterator_impl( stepper_type stepper , system_type sys , state_type &s ) + const_step_iterator_impl( stepper_type stepper , system_type sys , state_type& /* s */ ) : base_type( stepper , sys ) { } protected: diff --git a/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp b/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp index 147799040..f90f8875c 100644 --- a/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp +++ b/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp @@ -48,7 +48,7 @@ Time integrate_n_steps( Observer observer , stepper_tag ) { // ToDo: is there a better way to extract the final time? - Time t; + Time t = start_time; // Assignment is only here to avoid warnings. boost::for_each( make_n_step_time_range( stepper , system , start_state , start_time , dt , num_of_steps ) , obs_caller_time< Observer , Time >( t , observer ) ); @@ -91,7 +91,7 @@ Time integrate_n_steps( Observer observer , dense_output_stepper_tag ) { // ToDo: is there a better way to extract the final time? - Time t; + Time t = start_time; // Assignment is only here to avoid warnings. boost::for_each( make_n_step_time_range( stepper , system , start_state , start_time , dt , num_of_steps ) , obs_caller_time< Observer , Time >( t , observer ) ); diff --git a/boost/numeric/odeint/stepper/adams_bashforth.hpp b/boost/numeric/odeint/stepper/adams_bashforth.hpp index 59ed1c902..5ff1e8358 100644 --- a/boost/numeric/odeint/stepper/adams_bashforth.hpp +++ b/boost/numeric/odeint/stepper/adams_bashforth.hpp @@ -37,6 +37,7 @@ #include <boost/numeric/odeint/stepper/stepper_categories.hpp> #include <boost/numeric/odeint/stepper/runge_kutta4.hpp> +#include <boost/numeric/odeint/stepper/extrapolation_stepper.hpp> #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp> @@ -44,11 +45,28 @@ #include <boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp> #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp> +#include <boost/mpl/arithmetic.hpp> +#include <boost/mpl/min_max.hpp> +#include <boost/mpl/equal_to.hpp> + +namespace mpl = boost::mpl; + namespace boost { namespace numeric { namespace odeint { + using mpl::int_; + + /* if N >= 4, returns the smallest even number > N, otherwise returns 4 */ + template < int N > + struct order_helper + : mpl::max< typename mpl::eval_if< + mpl::equal_to< mpl::modulus< int_< N >, int_< 2 > >, + int_< 0 > >, + int_< N >, int_< N + 1 > >::type, + int_< 4 > >::type + { }; template< size_t Steps , @@ -59,7 +77,9 @@ class Time = Value , class Algebra = typename algebra_dispatcher< State >::algebra_type , class Operations = typename operations_dispatcher< State >::operations_type , class Resizer = initially_resizer , -class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer > +class InitializingStepper = extrapolation_stepper< order_helper<Steps>::value, + State, Value, Deriv, Time, + Algebra, Operations, Resizer > > class adams_bashforth : public algebra_stepper_base< Algebra , Operations > { @@ -181,7 +201,8 @@ public : { if( i != 0 ) m_step_storage.rotate(); sys( x , m_step_storage[0].m_v , t ); - stepper.do_step( system , x , m_step_storage[0].m_v , t , dt ); + stepper.do_step_dxdt_impl( system, x, m_step_storage[0].m_v, t, + dt ); t += dt; } m_steps_initialized = steps; @@ -222,7 +243,8 @@ private: { if( m_steps_initialized != 0 ) m_step_storage.rotate(); sys( in , m_step_storage[0].m_v , t ); - m_initializing_stepper.do_step( system , in , m_step_storage[0].m_v , t , out , dt ); + m_initializing_stepper.do_step_dxdt_impl( + system, in, m_step_storage[0].m_v, t, out, dt ); ++m_steps_initialized; } else diff --git a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp index 19b023999..2f7cc4c6f 100644 --- a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp +++ b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp @@ -49,7 +49,8 @@ class Deriv = State , class Time = Value , class Algebra = typename algebra_dispatcher< State >::algebra_type , class Operations = typename operations_dispatcher< State >::operations_type , -class Resizer = initially_resizer +class Resizer = initially_resizer, +class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer > > class adams_bashforth_moulton { @@ -71,12 +72,13 @@ public : typedef Operations operations_type; typedef Resizer resizer_type; typedef stepper_tag stepper_category; + typedef InitializingStepper initializing_stepper_type; static const size_t steps = Steps; #ifndef DOXYGEN_SKIP - typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_bashforth_type; + typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type, initializing_stepper_type > adams_bashforth_type; typedef adams_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_moulton_type; - typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > stepper_type; + typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type , initializing_stepper_type> stepper_type; #endif //DOXYGEN_SKIP typedef unsigned short order_type; static const order_type order_value = steps; @@ -158,7 +160,7 @@ private: { m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) ); m_adams_bashforth.do_step( system , x , t , m_x.m_v , dt ); - m_adams_moulton.do_step( system , x , m_x.m_v , t , x , dt , m_adams_bashforth.step_storage() ); + m_adams_moulton.do_step( system , x , m_x.m_v , t+dt , x , dt , m_adams_bashforth.step_storage() ); } else { diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp index 3ba11a4fe..08009dc16 100644 --- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp +++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp @@ -152,6 +152,22 @@ public: /* + * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template < class System, class StateInOut, class DerivIn > + void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt, + time_type t, time_type dt ) + { + this->stepper().do_step_impl( system , x , dxdt , t , x , dt ); + } + + + + /* * Version 3 : do_step( sys , in , t , out , dt ) * * this version does not solve the forwarding problem, boost.range can not be used @@ -181,10 +197,21 @@ public: { this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); } - - - - + + /* + * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_dxdt_impl( System system, const StateIn &in, + const DerivIn &dxdt, time_type t, StateOut &out, + time_type dt ) + { + this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); + } /* * Version 5 :do_step( sys , x , t , dt , xerr ) diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp index a055c7fa5..b1d751a0c 100644 --- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp +++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp @@ -151,11 +151,27 @@ public: /* + * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template< class System , class StateInOut , class DerivInOut > + void do_step_dxdt_impl( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt ) + { + m_first_call = true; + this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt ); + } + + /* * version 3 : do_step( sys , in , t , out , dt ) * - * this version does not solve the forwarding problem, boost.range can not be used + * this version does not solve the forwarding problem, boost.range can not + * be used. * - * the disable is needed to avoid ambiguous overloads if state_type = time_type + * the disable is needed to avoid ambiguous overloads if + * state_type = time_type */ template< class System , class StateIn , class StateOut > typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type @@ -174,12 +190,14 @@ public: * * this version does not solve the forwarding problem, boost.range can not be used */ - template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut > - void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t , - StateOut &out , DerivOut &dxdt_out , time_type dt ) + template< class System, class StateIn, class DerivIn, class StateOut, + class DerivOut > + void do_step( System system, const StateIn &in, const DerivIn &dxdt_in, + time_type t, StateOut &out, DerivOut &dxdt_out, time_type dt ) { m_first_call = true; - this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt ); + this->stepper().do_step_impl( system, in, dxdt_in, t, out, dxdt_out, + dt ); } diff --git a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp index 40aab8039..d81c8c7a7 100644 --- a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp +++ b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp @@ -139,6 +139,21 @@ public: /* + * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version that is not disabled. Meant for internal use only. + */ + template < class System, class StateInOut, class DerivIn > + void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt, + time_type t, time_type dt ) + { + this->stepper().do_step_impl( system , x , dxdt , t , x , dt ); + } + + + /* * Version 3 : do_step( sys , in , t , out , dt ) * * this version does not solve the forwarding problem, boost.range can not be used @@ -164,6 +179,22 @@ public: this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); } + + /* + * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt ) + * + * this version is needed when this stepper is used for initializing + * multistep stepper like adams-bashforth. Hence we provide an explicitely + * named version. Meant for internal use only. + */ + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_dxdt_impl( System system, const StateIn &in, + const DerivIn &dxdt, time_type t, StateOut &out, + time_type dt ) + { + this->stepper().do_step_impl( system , in , dxdt , t , out , dt ); + } + template< class StateIn > void adjust_size( const StateIn &x ) { diff --git a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp index ea3523e5b..eb09aefc8 100644 --- a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp +++ b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp @@ -182,7 +182,7 @@ private: // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q) template< class System , class StateIn , class StateOut > - void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , boost::mpl::true_ ) + void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , boost::mpl::true_ ) { typedef typename odeint::unwrap_reference< System >::type system_type; typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type; diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp index b604f86a2..e71008a99 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp @@ -194,7 +194,6 @@ public: static const value_type val1( 1.0 ); - typename odeint::unwrap_reference< System >::type &sys = system; if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) ) { reset(); // system resized -> reset @@ -219,12 +218,12 @@ public: m_midpoint.set_steps( m_interval_sequence[k] ); if( k == 0 ) { - m_midpoint.do_step( sys , in , dxdt , t , out , dt ); + m_midpoint.do_step( system , in , dxdt , t , out , dt ); /* the first step, nothing more to do */ } else { - m_midpoint.do_step( sys , in , dxdt , t , m_table[k-1].m_v , dt ); + m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt ); extrapolate( k , m_table , m_coeff , out ); // get error estimate m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v , @@ -341,7 +340,7 @@ public: resize_m_dxdt( x ); resize_m_xnew( x ); resize_impl( x ); - m_midpoint.adjust_size(); + m_midpoint.adjust_size( x ); } diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp index c352c9605..0036d1d2a 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp @@ -155,8 +155,6 @@ public: static const value_type val1( 1.0 ); - typename odeint::unwrap_reference< System >::type &sys = system; - bool reject( true ); time_vector h_opt( m_k_max+1 ); @@ -261,6 +259,7 @@ public: { //calculate dxdt for next step and dense output + typename odeint::unwrap_reference< System >::type &sys = system; sys( out , dxdt_new , t+dt ); //prepare dense output @@ -367,7 +366,7 @@ public: void adjust_size( const StateIn &x ) { resize_impl( x ); - m_midpoint.adjust_size(); + m_midpoint.adjust_size( x ); } diff --git a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp index 168bc6920..cc279abbe 100644 --- a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp +++ b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp @@ -117,7 +117,7 @@ struct adams_bashforth_call_algebra< 7 , Algebra , Operations > { //BOOST_ASSERT( false ); // not implemented typedef typename Coefficients::value_type value_type; - Algebra::for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , + algebra.for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , typename Operations::template scale_sum8< value_type , Time , Time , Time , Time , Time , Time >( 1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] ) ); } @@ -132,7 +132,7 @@ struct adams_bashforth_call_algebra< 8 , Algebra , Operations > { //BOOST_ASSERT( false ); // not implemented typedef typename Coefficients::value_type value_type; - Algebra::for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v , + algebra.for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v , typename Operations::template scale_sum9< value_type , Time , Time , Time , Time , Time , Time , Time >( 1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] , dt * coef[7] ) ); } diff --git a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp index f8a150551..b6f5f2a40 100644 --- a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp +++ b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp @@ -32,7 +32,7 @@ template< class Algebra , class Operations > struct adams_moulton_call_algebra< 1 , Algebra , Operations > { template< class StateIn , class StateOut , class DerivIn , class StepStorage , class Coefficients , class Time > - void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage &steps , const Coefficients &coef , Time dt ) const + void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage& /* steps */ , const Coefficients &coef , Time dt ) const { typedef typename Coefficients::value_type value_type; algebra.for_each3( out , in , dxdt , typename Operations::template scale_sum2< value_type , Time >( 1.0 , dt * coef[0] ) ); diff --git a/boost/numeric/odeint/stepper/euler.hpp b/boost/numeric/odeint/stepper/euler.hpp index 443f94276..1c7c126b9 100644 --- a/boost/numeric/odeint/stepper/euler.hpp +++ b/boost/numeric/odeint/stepper/euler.hpp @@ -76,7 +76,7 @@ public : { } template< class System , class StateIn , class DerivIn , class StateOut > - void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt ) + void do_step_impl( System /* system */ , const StateIn &in , const DerivIn &dxdt , time_type /* t */ , StateOut &out , time_type dt ) { stepper_base_type::m_algebra.for_each3( out , in , dxdt , typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt ) ); @@ -84,7 +84,7 @@ public : } template< class StateOut , class StateIn1 , class StateIn2 > - void calc_state( StateOut &x , time_type t , const StateIn1 &old_state , time_type t_old , const StateIn2 ¤t_state , time_type t_new ) const + void calc_state( StateOut &x , time_type t , const StateIn1 &old_state , time_type t_old , const StateIn2 & /*current_state*/ , time_type /* t_new */ ) const { const time_type delta = t - t_old; stepper_base_type::m_algebra.for_each3( x , old_state , stepper_base_type::m_dxdt.m_v , diff --git a/boost/numeric/odeint/stepper/extrapolation_stepper.hpp b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp new file mode 100644 index 000000000..3b688e958 --- /dev/null +++ b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp @@ -0,0 +1,293 @@ +/* + [auto_generated] + boost/numeric/odeint/stepper/extrapolation_stepper.hpp + + [begin_description] + extrapolation stepper + [end_description] + + Copyright 2009-2015 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) +*/ + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED + +#include <iostream> + +#include <algorithm> + +#include <boost/config.hpp> // for min/max guidelines +#include <boost/static_assert.hpp> + +#include <boost/numeric/odeint/util/bind.hpp> +#include <boost/numeric/odeint/util/unwrap_reference.hpp> + +#include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp> +#include <boost/numeric/odeint/stepper/modified_midpoint.hpp> +#include <boost/numeric/odeint/stepper/controlled_step_result.hpp> +#include <boost/numeric/odeint/algebra/range_algebra.hpp> +#include <boost/numeric/odeint/algebra/default_operations.hpp> +#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp> +#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp> + +#include <boost/numeric/odeint/util/state_wrapper.hpp> +#include <boost/numeric/odeint/util/is_resizeable.hpp> +#include <boost/numeric/odeint/util/resizer.hpp> +#include <boost/numeric/odeint/util/unit_helper.hpp> +#include <boost/numeric/odeint/util/detail/less_with_sign.hpp> + +namespace boost +{ +namespace numeric +{ +namespace odeint +{ + +template < unsigned short Order, class State, class Value = double, + class Deriv = State, class Time = Value, + class Algebra = typename algebra_dispatcher< State >::algebra_type, + class Operations = + typename operations_dispatcher< State >::operations_type, + class Resizer = initially_resizer > +#ifndef DOXYGEN_SKIP +class extrapolation_stepper + : public explicit_error_stepper_base< + extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra, + Operations, Resizer >, + Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, + Operations, Resizer > +#else +class extrapolation_stepper : public explicit_error_stepper_base +#endif +{ + + private: + // check for Order being odd + BOOST_STATIC_ASSERT_MSG( + ( ( Order % 2 ) == 0 ) && ( Order > 2 ), + "extrapolation_stepper requires even Order larger than 2" ); + + public: +#ifndef DOXYGEN_SKIP + typedef explicit_error_stepper_base< + extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra, + Operations, Resizer >, + Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, Operations, + Resizer > stepper_base_type; +#else + typedef explicit_error_stepper_base< extrapolation_stepper< ... >, ... > + stepper_base_type; +#endif + + typedef typename stepper_base_type::state_type state_type; + typedef typename stepper_base_type::value_type value_type; + typedef typename stepper_base_type::deriv_type deriv_type; + typedef typename stepper_base_type::time_type time_type; + typedef typename stepper_base_type::algebra_type algebra_type; + typedef typename stepper_base_type::operations_type operations_type; + typedef typename stepper_base_type::resizer_type resizer_type; + +#ifndef DOXYGEN_SKIP + typedef typename stepper_base_type::stepper_type stepper_type; + typedef typename stepper_base_type::wrapped_state_type wrapped_state_type; + typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type; + + typedef std::vector< value_type > value_vector; + typedef std::vector< value_vector > value_matrix; + typedef std::vector< size_t > int_vector; + typedef std::vector< wrapped_state_type > state_table_type; + typedef modified_midpoint< state_type, value_type, deriv_type, time_type, + algebra_type, operations_type, + resizer_type > midpoint_stepper_type; + +#endif // DOXYGEN_SKIP + + typedef unsigned short order_type; + static const order_type order_value = stepper_base_type::order_value; + static const order_type stepper_order_value = + stepper_base_type::stepper_order_value; + static const order_type error_order_value = + stepper_base_type::error_order_value; + + const static size_t m_k_max = ( order_value - 2 ) / 2; + + extrapolation_stepper( const algebra_type &algebra = algebra_type() ) + : stepper_base_type( algebra ), m_interval_sequence( m_k_max + 1 ), + m_coeff( m_k_max + 1 ), m_table( m_k_max ) + { + for ( unsigned short i = 0; i < m_k_max + 1; i++ ) + { + m_interval_sequence[i] = 2 * ( i + 1 ); + m_coeff[i].resize( i ); + for ( size_t k = 0; k < i; ++k ) + { + const value_type r = + static_cast< value_type >( m_interval_sequence[i] ) / + static_cast< value_type >( m_interval_sequence[k] ); + m_coeff[i][k] = + static_cast< value_type >( 1 ) / + ( r * r - static_cast< value_type >( + 1 ) ); // coefficients for extrapolation + } + } + } + + template < class System, class StateIn, class DerivIn, class StateOut, + class Err > + void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt, + time_type t, StateOut &out, time_type dt, Err &xerr ) + { + // std::cout << "dt: " << dt << std::endl; + // normal step + do_step_impl( system, in, dxdt, t, out, dt ); + + static const value_type val1( 1.0 ); + // additionally, perform the error calculation + stepper_base_type::m_algebra.for_each3( + xerr, out, m_table[0].m_v, + typename operations_type::template scale_sum2< + value_type, value_type >( val1, -val1 ) ); + } + + template < class System, class StateInOut, class DerivIn, class Err > + void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt, + time_type t, time_type dt, Err &xerr ) + { + // normal step + do_step_impl_io( system, inout, dxdt, t, dt ); + + static const value_type val1( 1.0 ); + // additionally, perform the error calculation + stepper_base_type::m_algebra.for_each3( + xerr, inout, m_table[0].m_v, + typename operations_type::template scale_sum2< + value_type, value_type >( val1, -val1 ) ); + } + + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt, + time_type t, StateOut &out, time_type dt ) + { + m_resizer.adjust_size( + in, detail::bind( &stepper_type::template resize_impl< StateIn >, + detail::ref( *this ), detail::_1 ) ); + size_t k = 0; + m_midpoint.set_steps( m_interval_sequence[k] ); + m_midpoint.do_step( system, in, dxdt, t, out, dt ); + for ( k = 1; k <= m_k_max; ++k ) + { + m_midpoint.set_steps( m_interval_sequence[k] ); + m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt ); + extrapolate( k, m_table, m_coeff, out ); + } + } + + template < class System, class StateInOut, class DerivIn > + void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt, + time_type t, time_type dt ) + { + // special care for inout + m_xout_resizer.adjust_size( + inout, + detail::bind( &stepper_type::template resize_m_xout< StateInOut >, + detail::ref( *this ), detail::_1 ) ); + do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt ); + boost::numeric::odeint::copy( m_xout.m_v, inout ); + } + + template < class System, class StateInOut, class DerivIn > + void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt, + time_type t, time_type dt ) + { + do_step_impl_io( system , x , dxdt , t , dt ); + } + + template < class System, class StateIn, class DerivIn, class StateOut > + void do_step_dxdt_impl( System system, const StateIn &in, + const DerivIn &dxdt, time_type t, StateOut &out, + time_type dt ) + { + do_step_impl( system , in , dxdt , t , out , dt ); + } + + + template < class StateIn > void adjust_size( const StateIn &x ) + { + resize_impl( x ); + m_midpoint.adjust_size( x ); + } + + private: + template < class StateIn > bool resize_impl( const StateIn &x ) + { + bool resized( false ); + for ( size_t i = 0; i < m_k_max; ++i ) + resized |= adjust_size_by_resizeability( + m_table[i], x, typename is_resizeable< state_type >::type() ); + return resized; + } + + template < class StateIn > bool resize_m_xout( const StateIn &x ) + { + return adjust_size_by_resizeability( + m_xout, x, typename is_resizeable< state_type >::type() ); + } + + template < class StateInOut > + void extrapolate( size_t k, state_table_type &table, + const value_matrix &coeff, StateInOut &xest ) + /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf + uses the obtained intermediate results to extrapolate to dt->0 + */ + { + static const value_type val1 = static_cast< value_type >( 1.0 ); + + for ( int j = k - 1; j > 0; --j ) + { + stepper_base_type::m_algebra.for_each3( + table[j - 1].m_v, table[j].m_v, table[j - 1].m_v, + typename operations_type::template scale_sum2< + value_type, value_type >( val1 + coeff[k][j], + -coeff[k][j] ) ); + } + stepper_base_type::m_algebra.for_each3( + xest, table[0].m_v, xest, + typename operations_type::template scale_sum2< + value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) ); + } + + private: + midpoint_stepper_type m_midpoint; + + resizer_type m_resizer; + resizer_type m_xout_resizer; + + int_vector m_interval_sequence; // stores the successive interval counts + value_matrix m_coeff; + + wrapped_state_type m_xout; + state_table_type m_table; // sequence of states for extrapolation +}; + +/******** DOXYGEN *******/ + +/** + * \class extrapolation_stepper + * \brief Extrapolation stepper with configurable order, and error estimation. + * + * The extrapolation stepper is a stepper with error estimation and configurable + * order. The order is given as template parameter and needs to be an _odd_ + * number. The stepper is based on several executions of the modified midpoint + * method and a Richardson extrapolation. This is essentially the same technique + * as for bulirsch_stoer, but without the variable order. + * + * \note The Order parameter has to be an even number greater 2. + */ +} +} +} +#endif diff --git a/boost/numeric/odeint/stepper/implicit_euler.hpp b/boost/numeric/odeint/stepper/implicit_euler.hpp index f3019889c..e1c64164b 100644 --- a/boost/numeric/odeint/stepper/implicit_euler.hpp +++ b/boost/numeric/odeint/stepper/implicit_euler.hpp @@ -144,7 +144,7 @@ private: void solve( state_type &x , matrix_type &m ) { int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v ); - if( res != 0 ) exit(0); + if( res != 0 ) std::exit(0); boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x ); } diff --git a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp index 4b3824fc0..df4e6c48b 100644 --- a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp +++ b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp @@ -66,6 +66,7 @@ public: { BOOST_USING_STD_MAX(); using std::abs; + using std::sqrt; const size_t n = x.size(); value_type err = 0.0 , sk = 0.0; diff --git a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp index 2b8e70b16..32bda0bd7 100644 --- a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp +++ b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp @@ -121,11 +121,20 @@ public : // dt * m_dxh = k4 sys( m_x_tmp.m_v , m_dxh.m_v , t + dt ); + //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm ) time_type dt6 = dt / static_cast< value_type >( 6 ); time_type dt3 = dt / static_cast< value_type >( 3 ); stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_dxt.m_v , m_dxm.m_v , m_dxh.m_v , - typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) ); + typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) ); + + // x += dt/6 * m_dxdt + dt/3 * m_dxt ) + // stepper_base_type::m_algebra.for_each4( out , in , dxdt , m_dxt.m_v , + // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt6 , dt3 ) ); + // // x += dt/3 * m_dxm + dt/6 * m_dxh ) + // stepper_base_type::m_algebra.for_each4( out , out , m_dxm.m_v , m_dxh.m_v , + // typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt3 , dt6 ) ); + } template< class StateType > diff --git a/boost/numeric/odeint/stepper/velocity_verlet.hpp b/boost/numeric/odeint/stepper/velocity_verlet.hpp index 24f3c0e3b..3a20fc25c 100644 --- a/boost/numeric/odeint/stepper/velocity_verlet.hpp +++ b/boost/numeric/odeint/stepper/velocity_verlet.hpp @@ -114,7 +114,7 @@ public: algebra_stepper_base_type::m_algebra.for_each4( qout , qin , pin , ain , - typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one_half * dt , one * dt * dt ) ); + typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) ); typename odeint::unwrap_reference< System >::type & sys = system; diff --git a/boost/numeric/odeint/util/bind.hpp b/boost/numeric/odeint/util/bind.hpp index 966fd0397..1201afab8 100644 --- a/boost/numeric/odeint/util/bind.hpp +++ b/boost/numeric/odeint/util/bind.hpp @@ -41,9 +41,12 @@ using namespace ::std::placeholders; #else +// unnamed namespace to avoid multiple declarations (#138) +namespace { using ::boost::bind; boost::arg<1> _1; boost::arg<2> _2; +} // using ::boost::bind; // using ::_1; // using ::_2; diff --git a/boost/numeric/odeint/util/resize.hpp b/boost/numeric/odeint/util/resize.hpp index db2ce208b..75e9eff44 100644 --- a/boost/numeric/odeint/util/resize.hpp +++ b/boost/numeric/odeint/util/resize.hpp @@ -33,8 +33,8 @@ namespace boost { namespace numeric { namespace odeint { - - + + template< class StateOut , class StateIn , class Enabler = void > struct resize_impl_sfinae { diff --git a/boost/numeric/odeint/util/same_instance.hpp b/boost/numeric/odeint/util/same_instance.hpp index dc55db30a..a889ee192 100644 --- a/boost/numeric/odeint/util/same_instance.hpp +++ b/boost/numeric/odeint/util/same_instance.hpp @@ -25,7 +25,7 @@ namespace odeint { template< class T1 , class T2 , class Enabler=void > struct same_instance_impl { - static bool same_instance( const T1 &x1 , const T2 &x2 ) + static bool same_instance( const T1& /* x1 */ , const T2& /* x2 */ ) { return false; } diff --git a/boost/numeric/odeint/util/unwrap_reference.hpp b/boost/numeric/odeint/util/unwrap_reference.hpp index d622a5c00..bc7d423a0 100644 --- a/boost/numeric/odeint/util/unwrap_reference.hpp +++ b/boost/numeric/odeint/util/unwrap_reference.hpp @@ -31,9 +31,9 @@ namespace boost { #if BOOST_NUMERIC_ODEINT_CXX11 -template<typename T> class reference_wrapper; +template<typename T> struct reference_wrapper; -template<typename T> class unwrap_reference; +template<typename T> struct unwrap_reference; #endif namespace numeric { @@ -43,24 +43,21 @@ namespace odeint { #if BOOST_NUMERIC_ODEINT_CXX11 template<typename T> -class unwrap_reference +struct unwrap_reference { -public: typedef typename std::remove_reference<T>::type type; }; template<typename T> -class unwrap_reference< std::reference_wrapper<T> > +struct unwrap_reference< std::reference_wrapper<T> > { -public: typedef typename std::remove_reference<T>::type type; }; template<typename T> -class unwrap_reference< boost::reference_wrapper<T> > +struct unwrap_reference< boost::reference_wrapper<T> > { -public: - typedef typename boost::unwrap_reference<T>::type type; + typedef typename boost::unwrap_reference<T>::type type; }; #else diff --git a/boost/numeric/ublas/detail/matrix_assign.hpp b/boost/numeric/ublas/detail/matrix_assign.hpp index 5dea5920f..be172dd63 100644 --- a/boost/numeric/ublas/detail/matrix_assign.hpp +++ b/boost/numeric/ublas/detail/matrix_assign.hpp @@ -1259,20 +1259,20 @@ namespace detail { void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, row_major_tag) { typedef F<typename M::iterator2::reference, typename E::reference> functor_type; // R unnecessary, make_conformant not required - typedef typename M::size_type size_type; + //typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right typedef typename M::difference_type difference_type; typename M::iterator1 it1 (m.begin1 ()); typename E::iterator1 it1e (e ().begin1 ()); - difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (e ().end1 () - it1e))); + difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (e ().end1 () - it1e))); while (-- size1 >= 0) { #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION typename M::iterator2 it2 (it1.begin ()); typename E::iterator2 it2e (it1e.begin ()); - difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (it1e.end () - it2e))); + difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (it1e.end () - it2e))); #else typename M::iterator2 it2 (begin (it1, iterator1_tag ())); typename E::iterator2 it2e (begin (it1e, iterator1_tag ())); - difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (end (it1e, iterator1_tag ()) - it2e))); + difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (end (it1e, iterator1_tag ()) - it2e))); #endif while (-- size2 >= 0) functor_type::apply (*it2, *it2e), ++ it2, ++ it2e; @@ -1285,20 +1285,20 @@ namespace detail { void matrix_swap (M &m, matrix_expression<E> &e, dense_proxy_tag, column_major_tag) { typedef F<typename M::iterator1::reference, typename E::reference> functor_type; // R unnecessary, make_conformant not required - typedef typename M::size_type size_type; + // typedef typename M::size_type size_type; // gcc is complaining that this is not used, although this is not right typedef typename M::difference_type difference_type; typename M::iterator2 it2 (m.begin2 ()); typename E::iterator2 it2e (e ().begin2 ()); - difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), size_type (e ().end2 () - it2e))); + difference_type size2 (BOOST_UBLAS_SAME (m.size2 (), typename M::size_type (e ().end2 () - it2e))); while (-- size2 >= 0) { #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION typename M::iterator1 it1 (it2.begin ()); typename E::iterator1 it1e (it2e.begin ()); - difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (it2e.end () - it1e))); + difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (it2e.end () - it1e))); #else typename M::iterator1 it1 (begin (it2, iterator2_tag ())); typename E::iterator1 it1e (begin (it2e, iterator2_tag ())); - difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), size_type (end (it2e, iterator2_tag ()) - it1e))); + difference_type size1 (BOOST_UBLAS_SAME (m.size1 (), typename M::size_type (end (it2e, iterator2_tag ()) - it1e))); #endif while (-- size1 >= 0) functor_type::apply (*it1, *it1e), ++ it1, ++ it1e; diff --git a/boost/numeric/ublas/experimental/sparse_view.hpp b/boost/numeric/ublas/experimental/sparse_view.hpp index 52ad509e9..3a1641084 100644 --- a/boost/numeric/ublas/experimental/sparse_view.hpp +++ b/boost/numeric/ublas/experimental/sparse_view.hpp @@ -19,6 +19,7 @@ #include <boost/next_prior.hpp> #include <boost/type_traits/remove_cv.hpp> +#include <boost/numeric/ublas/storage.hpp> namespace boost { namespace numeric { namespace ublas { diff --git a/boost/numeric/ublas/functional.hpp b/boost/numeric/ublas/functional.hpp index e6f92731f..c7c92614c 100644 --- a/boost/numeric/ublas/functional.hpp +++ b/boost/numeric/ublas/functional.hpp @@ -15,6 +15,8 @@ #include <functional> +#include <boost/core/ignore_unused.hpp> + #include <boost/numeric/ublas/traits.hpp> #ifdef BOOST_UBLAS_USE_DUFF_DEVICE #include <boost/numeric/ublas/detail/duff.hpp> @@ -1161,6 +1163,7 @@ namespace boost { namespace numeric { namespace ublas { #elif defined(BOOST_UBLAS_HAVE_BINDINGS) return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j)); #else + boost::ignore_unused(j); return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i)); #endif } diff --git a/boost/numeric/ublas/matrix_expression.hpp b/boost/numeric/ublas/matrix_expression.hpp index 3953fad4c..a36313096 100644 --- a/boost/numeric/ublas/matrix_expression.hpp +++ b/boost/numeric/ublas/matrix_expression.hpp @@ -3946,7 +3946,9 @@ namespace boost { namespace numeric { namespace ublas { // (m / t) [i] [j] = m [i] [j] / t template<class E1, class T2> BOOST_UBLAS_INLINE + typename enable_if< is_convertible<T2, typename E1::value_type>, typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type + >::type operator / (const matrix_expression<E1> &e1, const T2 &e2) { typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type; diff --git a/boost/numeric/ublas/storage.hpp b/boost/numeric/ublas/storage.hpp index 9dfaaee07..8821309dd 100644 --- a/boost/numeric/ublas/storage.hpp +++ b/boost/numeric/ublas/storage.hpp @@ -370,7 +370,7 @@ namespace boost { namespace numeric { namespace ublas { // Random Access Container BOOST_UBLAS_INLINE size_type max_size () const { - return ALLOC ().max_size(); + return N; } BOOST_UBLAS_INLINE @@ -537,6 +537,10 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE array_adaptor (size_type size, pointer data): size_ (size), own_ (false), data_ (data) {} + + template <size_t N> + BOOST_UBLAS_INLINE array_adaptor (T (&data)[N]): + size_ (N), own_ (false), data_ (data) {} BOOST_UBLAS_INLINE array_adaptor (const array_adaptor &a): storage_array<self_type> (), @@ -602,6 +606,16 @@ namespace boost { namespace numeric { namespace ublas { resize_internal (size, data, init, true); } + template <size_t N> + BOOST_UBLAS_INLINE void resize (T (&data)[N]) { + resize_internal (N, data, value_type (), false); + } + + template <size_t N> + BOOST_UBLAS_INLINE void resize (T (&data)[N], value_type init) { + resize_internal (N, data, init, true); + } + BOOST_UBLAS_INLINE size_type size () const { return size_; @@ -734,7 +748,7 @@ namespace boost { namespace numeric { namespace ublas { typedef TT *argument_type; BOOST_UBLAS_INLINE - result_type operator () (argument_type x) {} + result_type operator () (argument_type /* x */) {} }; public: @@ -763,6 +777,10 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE shallow_array_adaptor (size_type size, pointer data): size_ (size), own_ (false), data_ (data, leaker<value_type> ()) {} + BOOST_UBLAS_INLINE + template <size_t N> + shallow_array_adaptor (T (&data)[N]): + size_ (N), own_ (false), data_ (data, leaker<value_type> ()) {} BOOST_UBLAS_INLINE shallow_array_adaptor (const shallow_array_adaptor &a): @@ -784,6 +802,7 @@ namespace boost { namespace numeric { namespace ublas { std::fill (data.get () + (std::min) (size, size_), data.get () + size, init); } size_ = size; + own_ = true; data_ = data; } } @@ -794,7 +813,8 @@ namespace boost { namespace numeric { namespace ublas { std::fill (data + (std::min) (size, size_), data + size, init); } size_ = size; - data_ = data; + own_ = false; + data_.reset(data, leaker<value_type> ()); } public: BOOST_UBLAS_INLINE @@ -813,6 +833,16 @@ namespace boost { namespace numeric { namespace ublas { void resize (size_type size, pointer data, value_type init) { resize_internal (size, data, init, true); } + BOOST_UBLAS_INLINE + template <size_t N> + void resize (T (&data)[N]) { + resize_internal (N, data, value_type (), false); + } + BOOST_UBLAS_INLINE + template <size_t N> + void resize (T (&data)[N], value_type init) { + resize_internal (N, data, init, true); + } BOOST_UBLAS_INLINE size_type size () const { diff --git a/boost/numeric/ublas/traits.hpp b/boost/numeric/ublas/traits.hpp index 86ae9b50f..ecd52cad8 100644 --- a/boost/numeric/ublas/traits.hpp +++ b/boost/numeric/ublas/traits.hpp @@ -20,6 +20,9 @@ #include <boost/numeric/ublas/detail/config.hpp> #include <boost/numeric/ublas/detail/iterator.hpp> #include <boost/numeric/ublas/detail/returntype_deduction.hpp> +#ifdef BOOST_UBLAS_USE_INTERVAL +#include <boost/numeric/interval.hpp> +#endif #include <boost/type_traits.hpp> #include <complex> diff --git a/boost/numeric/ublas/vector.hpp b/boost/numeric/ublas/vector.hpp index c90e42961..c1384c145 100644 --- a/boost/numeric/ublas/vector.hpp +++ b/boost/numeric/ublas/vector.hpp @@ -27,6 +27,9 @@ #ifdef BOOST_UBLAS_CPP_GE_2011 #include <array> #include <initializer_list> +#if defined(BOOST_MSVC) // For std::forward in fixed_vector +#include <utility> +#endif #endif // Iterators based on ideas of Jeremy Siek @@ -903,10 +906,18 @@ namespace boost { namespace numeric { namespace ublas { /// \brief Construct a fixed_vector from a list of values /// This constructor enables initialization by using any of: /// fixed_vector<double, 3> v = { 1, 2, 3 } or fixed_vector<double,3> v( {1, 2, 3} ) or fixed_vector<double,3> v( 1, 2, 3 ) +#if defined(BOOST_MSVC) + // This may or may not work. Maybe use this for all instead only for MSVC + template <typename... U> + fixed_vector(U&&... values) : + vector_container<self_type> (), + data_{{ std::forward<U>(values)... }} {} +#else template <typename... Types> fixed_vector(value_type v0, Types... vrest) : vector_container<self_type> (), data_{ { v0, vrest... } } {} +#endif // ----------------------- // Random Access Container |