diff --git a/Jamroot b/Jamroot index e7b92234..882251e2 100644 --- a/Jamroot +++ b/Jamroot @@ -19,7 +19,6 @@ project # tests, regression tests and examples build-project libs/numeric/odeint/test ; build-project libs/numeric/odeint/examples ; -build-project libs/numeric/odeint/regression_test ; # additional tests with external libraries : @@ -28,14 +27,6 @@ build-project libs/numeric/odeint/regression_test ; # build-project libs/numeric/odeint/test_external/gsl ; -# ideas -# build-project libs/numeric/odeint/ideas/butcher ; -# build-project libs/numeric/odeint/ideas/generic_stepper ; -#build-project libs/numeric/odeint/ideas/rosenbrock4 ; -#build-project libs/numeric/odeint/ideas/units ; -#build-project libs/numeric/odeint/ideas/algebra ; - - # docs: # build-project libs/numeric/odeint/doc ; diff --git a/boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp b/boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp new file mode 100644 index 00000000..3588ed20 --- /dev/null +++ b/boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp @@ -0,0 +1,162 @@ +/* +[begin_description] +Modification of the implicit Euler method, works with the MTL4 matrix library only. +[end_description] + +Copyright 2009-2011 Karsten Ahnert +Copyright 2009-2011 Mario Mulansky +Copyright 2012 Andreas Angelopoulos + +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_IMPLICIT_EULER_MTL4_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED + + +#include + +#include +#include +#include + +#include + +#include +#include + + + + +namespace boost { +namespace numeric { +namespace odeint { + + +template< class ValueType , class Resizer = initially_resizer > +class implicit_euler_mtl4 +{ + +public: + + typedef ValueType value_type; + typedef value_type time_type; + typedef mtl::dense_vector state_type; + + typedef state_wrapper< state_type > wrapped_state_type; + typedef state_type deriv_type; + typedef state_wrapper< deriv_type > wrapped_deriv_type; + typedef mtl::compressed2D< value_type > matrix_type; + typedef state_wrapper< matrix_type > wrapped_matrix_type; + + typedef Resizer resizer_type; + typedef stepper_tag stepper_category; + + typedef implicit_euler_mtl4< ValueType , Resizer > stepper_type; + + + implicit_euler_mtl4( const value_type epsilon = 1E-6 ) + : m_epsilon( epsilon ) , m_resizer() , + m_dxdt() , m_x() , + m_identity() , m_jacobi() + { } + + + template< class System > + void do_step( System system , state_type &x , time_type t , time_type dt ) + { + typedef typename odeint::unwrap_reference< System >::type system_type; + typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type; + typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type; + system_type &sys = system; + deriv_func_type &deriv_func = sys.first; + jacobi_func_type &jacobi_func = sys.second; + + m_resizer.adjust_size( x , detail::bind( + &stepper_type::template resize_impl< state_type > , detail::ref( *this ) , detail::_1 ) ); + + m_identity.m_v = 1; + + t += dt; + m_x.m_v = x; + + deriv_func( x , m_dxdt.m_v , t ); + jacobi_func( x , m_jacobi.m_v , t ); + + + m_dxdt.m_v *= -dt; + + m_jacobi.m_v *= dt; + m_jacobi.m_v -= m_identity.m_v ; + + + + // using ilu_0 preconditioning -incomplete LU factorisation + // itl::pc::diagonal L(m_jacobi.m_v); + itl::pc::ilu_0 L( m_jacobi.m_v ); + + solve( m_jacobi.m_v , m_x.m_v , m_dxdt.m_v , L ); + x+= m_x.m_v; + + + } + + + template< class StateType > + void adjust_size( const StateType &x ) + { + resize_impl( x ); + } + + +private: + + + /* + Applying approximate iteractive linear solvers + default solver is Biconjugate gradient stabilized method + itl::bicgstab(A, x, b, L, iter); + */ + template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner> + void solve(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, + const Preconditioner& L, int max_iteractions =500) + { + // Termination criterion: r < 1e-6 * b or N iterations + itl::basic_iteration< double > iter( b , max_iteractions , 1e-6 ); + itl::bicgstab( A , x , b , L , iter ); + + } + + + template< class StateIn > + bool resize_impl( const StateIn &x ) + { + bool resized = false; + resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable::type() ); + resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable::type() ); + resized |= adjust_size_by_resizeability( m_identity , x , typename is_resizeable::type() ); + resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable::type() ); + return resized; + } + + +private: + + value_type m_epsilon; + resizer_type m_resizer; + wrapped_deriv_type m_dxdt; + wrapped_state_type m_x; + wrapped_matrix_type m_identity; + wrapped_matrix_type m_jacobi; +}; + + +} // odeint +} // numeric +} // boost + + +#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED diff --git a/boost/numeric/odeint/external/mtl4/mtl4_resize.hpp b/boost/numeric/odeint/external/mtl4/mtl4_resize.hpp new file mode 100644 index 00000000..93569747 --- /dev/null +++ b/boost/numeric/odeint/external/mtl4/mtl4_resize.hpp @@ -0,0 +1,133 @@ +/* +[begin_description] +Modification of the implicit Euler method, works with the MTL4 matrix library only. +[end_description] + +Copyright 2009-2011 Karsten Ahnert +Copyright 2009-2011 Mario Mulansky +Copyright 2012 Andreas Angelopoulos + +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_MTL4_RESIZE_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_EXTERNAL_MTL4_RESIZE_HPP_INCLUDED + +#include +#include +#include + +#include +#include +#include + + +namespace boost { +namespace numeric { +namespace odeint { + + +template< class Value , class Parameters > +struct is_resizeable< mtl::dense_vector< Value , Parameters > > +{ + typedef boost::true_type type; + const static bool value = type::value; +}; + +template< class Value , class Parameters > +struct is_resizeable< mtl::dense2D< Value , Parameters > > +{ + typedef boost::true_type type; + const static bool value = type::value; +}; + +template< class Value , class Parameters > +struct is_resizeable< mtl::compressed2D< Value , Parameters > > +{ + typedef boost::true_type type; + const static bool value = type::value; +}; + + + + +template< class Value , class Parameters > +struct same_size_impl< mtl::dense_vector< Value , Parameters > , mtl::dense_vector< Value , Parameters > > +{ + static bool same_size( const mtl::dense_vector< Value , Parameters > &v1 , + const mtl::dense_vector< Value , Parameters > &v2 ) + { + return mtl::size( v1 ) == mtl::size( v2 ); + } +}; + +template< class Value , class Parameters > +struct resize_impl< mtl::dense_vector< Value , Parameters > , mtl::dense_vector< Value , Parameters > > +{ + static void resize( mtl::dense_vector< Value , Parameters > &v1 , + const mtl::dense_vector< Value , Parameters > &v2 ) + { + v1.change_dim( mtl::size( v2 ) ); + } +}; + + + +template< class Value , class MatrixParameters , class VectorParameters > +struct same_size_impl< mtl::dense2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > > +{ + static bool same_size( const mtl::dense2D< Value , MatrixParameters > &m , + const mtl::dense_vector< Value , VectorParameters > &v ) + { + return ( ( mtl::size( v ) == m.num_cols() ) && ( mtl::size( v ) == m.num_rows() ) ); + } +}; + +template< class Value , class MatrixParameters , class VectorParameters > +struct resize_impl< mtl::dense2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > > +{ + static void resize( mtl::dense2D< Value , MatrixParameters > &m , + const mtl::dense_vector< Value , VectorParameters > &v ) + { + m.change_dim( mtl::size( v ) , mtl::size( v ) , false ); + } +}; + + + + +template< class Value , class MatrixParameters , class VectorParameters > +struct same_size_impl< mtl::compressed2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > > +{ + static bool same_size( const mtl::compressed2D< Value , MatrixParameters > &m , + const mtl::dense_vector< Value , VectorParameters > &v ) + { + return ( ( mtl::size( v ) == m.num_cols() ) && ( mtl::size( v ) == m.num_rows() ) ); + } +}; + +template< class Value , class MatrixParameters , class VectorParameters > +struct resize_impl< mtl::compressed2D< Value , MatrixParameters > , mtl::dense_vector< Value , VectorParameters > > +{ + static void resize( mtl::compressed2D< Value , MatrixParameters > &m , + const mtl::dense_vector< Value , VectorParameters > &v ) + { + m.change_dim( mtl::size( v ) , mtl::size( v ) ); + } +}; + + + + + + + + +} // namespace odeint +} // namesapce numeric +} // namespace boost + +#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_MTL4_RESIZE_HPP_INCLUDED diff --git a/boost/numeric/odeint/external/thrust/thrust_operations.hpp b/boost/numeric/odeint/external/thrust/thrust_operations.hpp index be79f17a..710657ae 100644 --- a/boost/numeric/odeint/external/thrust/thrust_operations.hpp +++ b/boost/numeric/odeint/external/thrust/thrust_operations.hpp @@ -46,6 +46,25 @@ struct thrust_operations } }; + template< class Fac1 = double , class Fac2 = Fac1 > + struct scale_sum_swap2 + { + const Fac1 m_alpha1; + const Fac2 m_alpha2; + + scale_sum_swap2( const Fac1 alpha1 , const Fac2 alpha2 ) + : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { } + + template< class Tuple > + __host__ __device__ + void operator()( Tuple t ) const + { + typename thrust::tuple_element<0,Tuple>::type tmp = thrust::get<0>(t); + thrust::get<0>(t) = m_alpha1 * thrust::get<1>(t) + m_alpha2 * thrust::get<2>(t); + thrust::get<1>(t) = tmp; + } + }; + template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 > struct scale_sum3 { diff --git a/boost/numeric/odeint/integrate/detail/integrate_const.hpp b/boost/numeric/odeint/integrate/detail/integrate_const.hpp new file mode 100644 index 00000000..2eb7ae0a --- /dev/null +++ b/boost/numeric/odeint/integrate/detail/integrate_const.hpp @@ -0,0 +1,159 @@ +/* + [auto_generated] + boost/numeric/odeint/integrate/detail/integrate_const.hpp + + [begin_description] + integrate const implementation + [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_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_INTEGRATE_DETAIL_INTEGRATE_CONST_HPP_INCLUDED + +#include +#include +#include +#include + +#include + +namespace boost { +namespace numeric { +namespace odeint { +namespace detail { + +// forward declaration +template< class Stepper , class System , class State , class Time , class Observer > +size_t integrate_adaptive( + Stepper stepper , System system , State &start_state , + Time &start_time , Time end_time , Time &dt , + Observer observer , controlled_stepper_tag +); + + +template< class Stepper , class System , class State , class Time , class Observer > +size_t integrate_const( + Stepper stepper , System system , State &start_state , + Time start_time , Time end_time , Time dt , + Observer observer , stepper_tag +) +{ + typename odeint::unwrap_reference< Observer >::type &obs = observer; + + Time time = start_time; + int step = 0; + + while( less_eq_with_sign( time+dt , end_time , dt ) ) + { + obs( start_state , time ); + stepper.do_step( system , start_state , time , dt ); + // direct computation of the time avoids error propagation happening when using time += dt + // we need clumsy type analysis to get boost units working here + ++step; + time = start_time + static_cast< typename unit_value_type