Subject: [Boost-commit] svn:boost r66656 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper
From: mario.mulansky_at_[hidden]
Date: 2010-11-21 12:29:32
Author: mariomulansky
Date: 2010-11-21 12:29:28 EST (Sun, 21 Nov 2010)
New Revision: 66656
URL: http://svn.boost.org/trac/boost/changeset/66656
Log:
+ fusion stepper
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp | 15 +++++++++++++++
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/test.cpp | 11 ++++++++++-
2 files changed, 25 insertions(+), 1 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp 2010-11-21 12:29:28 EST (Sun, 21 Nov 2010)
@@ -42,4 +42,19 @@
};
+template< typename state_type>
+struct algebra< state_type , 2 >
+{
+
+ typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
+ typedef typename boost::numeric::odeint::standard_operations< double > std_op;
+
+ static void foreach( state_type &x_tmp , const state_type &x , const std::vector< double > &a , const state_type *k_vector , const double dt )
+ {
+ std_algebra::for_each4( x_tmp , x , k_vector[0] , k_vector[1] , std_op::scale_sum3( 1.0 , a[0]*dt , a[1]*dt ) );
+ }
+
+};
+
+
#endif /* ALGEBRA_HPP_ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp 2010-11-21 12:29:28 EST (Sun, 21 Nov 2010)
@@ -0,0 +1,273 @@
+/*
+ * fusion_stepper.hpp
+ *
+ * Created on: Nov 21, 2010
+ * Author: mario
+ */
+
+#ifndef FUSION_STEPPER_HPP_
+#define FUSION_STEPPER_HPP_
+
+#include <boost/mpl/vector.hpp>
+#include <boost/mpl/size.hpp>
+#include <boost/mpl/push_back.hpp>
+#include <boost/mpl/for_each.hpp>
+#include <boost/mpl/zip_view.hpp>
+#include <boost/mpl/range_c.hpp>
+#include <boost/mpl/int.hpp>
+#include <boost/mpl/at.hpp>
+#include <boost/mpl/copy.hpp>
+#include <boost/mpl/insert_range.hpp>
+
+#include <vector>
+
+#include <boost/fusion/container.hpp>
+#include <boost/fusion/algorithm.hpp>
+#include <boost/fusion/include/mpl.hpp>
+
+#include <typeinfo>
+
+
+#include "algebra.hpp"
+
+namespace mpl = boost::mpl;
+namespace fusion = boost::fusion;
+
+using namespace std;
+
+template< typename State , size_t stage >
+class stepper_stage
+{
+
+public:
+
+ typedef double value_type;
+ typedef State state_type;
+ typedef vector< vector< double > > parameter_array_type2D;
+ typedef vector< double > parameter_array_type1D;
+
+
+ void init( const parameter_array_type2D &a , const parameter_array_type1D &c )
+ {
+ m_a_row = a[stage-1];
+ m_c = c[stage-1];
+ }
+
+ template< typename System >
+ void operator() ( System &system , state_type &x_tmp , const state_type &x , state_type *k_vector ,
+ double t , const double dt ) const
+ {
+ system( x_tmp , k_vector[stage-1] , t + m_c * dt );
+ algebra<state_type , stage>::foreach< stage >( x_tmp , x , m_a_row , k_vector , dt);
+ }
+
+private:
+
+ vector< value_type > m_a_row;
+ value_type m_c;
+
+};
+
+template< typename State >
+class stepper_stage< State , 1 >
+{
+
+public:
+
+ typedef double value_type;
+ typedef State state_type;
+ typedef vector< vector< double > > parameter_array_type2D;
+ typedef vector< double > parameter_array_type1D;
+
+
+ void init( const parameter_array_type2D &a , const parameter_array_type1D &c )
+ {
+ m_a_row = a[0];
+ m_c = c[0];
+ }
+
+ template< typename System >
+ void operator() ( System &system , state_type &x_tmp , const state_type &x , state_type *k_vector ,
+ double t , const double dt ) const
+ {
+ system( x , k_vector[0] , t + m_c * dt );
+ algebra<state_type , 1>::foreach( x_tmp , x , m_a_row , k_vector , dt);
+ }
+
+private:
+
+ vector< value_type > m_a_row;
+ value_type m_c;
+
+};
+
+
+template< typename State , size_t stage >
+class stepper_last_stage
+{
+
+public:
+
+ typedef double value_type;
+ typedef State state_type;
+ typedef vector< vector< double > > parameter_array_type2D;
+ typedef vector< double > parameter_array_type1D;
+
+
+ void init( const parameter_array_type1D &b , const parameter_array_type1D &c )
+ {
+ m_b = b;
+ m_c = c[stage-1];
+ }
+
+ template< typename System >
+ void operator() ( System &system , state_type &x_tmp , state_type &x , state_type *k_vector , double t , const double dt )
+ {
+ system( x_tmp , k_vector[stage-1] , t + m_c * dt );
+ algebra<state_type , stage>::foreach( x , x , m_b , k_vector , dt);
+ }
+
+private:
+
+ vector< value_type > m_b;
+ value_type m_c;
+
+};
+
+
+template< typename State >
+class stepper_last_stage< State , 1 >
+{
+
+public:
+
+ typedef double value_type;
+ typedef State state_type;
+ typedef vector< vector< double > > parameter_array_type2D;
+ typedef vector< double > parameter_array_type1D;
+
+
+ void init( const parameter_array_type1D &b , const parameter_array_type1D &c )
+ {
+ m_b = b;
+ m_c = c[0];
+ }
+
+ template< typename System >
+ void operator() ( System &system , state_type &x_tmp , state_type &x , state_type k_vector[1] , double t , const double dt )
+ {
+ system( x , k_vector[0] , t + m_c * dt );
+ algebra<state_type , 1>::foreach( x , x , m_b , k_vector , dt);
+ }
+
+private:
+
+ vector< value_type > m_b;
+ value_type m_c;
+
+};
+
+
+
+
+
+
+
+
+template< class T , class Constant >
+struct stepper_stage_wrapper
+{
+ typedef stepper_stage< T , Constant::value > type;
+};
+
+
+
+
+/* Runge Kutta Stepper - consisting of several stages */
+
+template< typename State , size_t stage_count >
+class runge_kutta_stepper
+{
+
+ typedef State state_type;
+ typedef vector< vector< double > > parameter_array_type2D;
+ typedef vector< double > parameter_array_type1D;
+
+ template< class System >
+ struct calculate_stage
+ {
+ System &system;
+ state_type &x , &x_tmp;
+ state_type *k_vector;
+ const double t;
+ const double dt;
+
+ calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector ,
+ const double _t , const double _dt )
+ : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
+ {}
+
+ template< class Stage >
+ void operator()( Stage const &stage ) const
+ {
+ stage( system , x_tmp , x , k_vector , t , dt );
+ }
+
+ };
+
+ struct init_stage {
+
+ const parameter_array_type2D &m_a;
+ const parameter_array_type1D &m_c;
+
+ init_stage( const parameter_array_type2D &a , const parameter_array_type1D &c )
+ : m_a( a ) , m_c( c )
+ { }
+
+ template< typename Stage >
+ void operator() ( Stage &stage ) const
+ {
+ stage.init( m_a , m_c );
+ }
+ };
+
+public:
+
+ typedef mpl::range_c< size_t , 1 , stage_count > indices;
+ typedef typename fusion::result_of::as_vector <
+ typename mpl::copy // create mpl::vector< stepper_stage< 0 >, stepper_stage< 1 > , .. stepper_stage< stage_count-1 > >
+ <
+ indices ,
+ mpl::inserter
+ <
+ mpl::vector0<> ,
+ mpl::push_back< mpl::_1 , stepper_stage_wrapper< State , mpl::_2 > >
+ >
+ >::type
+ >::type stage_vector;
+ typedef stepper_last_stage< State , stage_count > last_stage_type;
+
+ runge_kutta_stepper( const parameter_array_type2D &a ,
+ const parameter_array_type1D &b ,
+ const parameter_array_type1D &c )
+ {
+ fusion::for_each( m_stages , init_stage( a , c ) );
+ m_last_stage.init( b , c );
+ }
+
+ template< class System >
+ void do_step( System &system , state_type &x , double t , const double dt )
+ {
+ fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
+ m_last_stage( system , m_x_tmp , x , m_k_vector , t , dt );
+ }
+
+private:
+
+ state_type m_k_vector[stage_count];
+ state_type m_x_tmp;
+ stage_vector m_stages;
+ last_stage_type m_last_stage;
+};
+
+#endif /* FUSION_STEPPER_HPP_ */
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/test.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/test.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/test.cpp 2010-11-21 12:29:28 EST (Sun, 21 Nov 2010)
@@ -9,12 +9,13 @@
#include <vector>
#include <tr1/array>
-#include "generic_stepper.hpp"
+#include "fusion_stepper.hpp"
using namespace std;
typedef tr1::array< double , 3 > state_type;
typedef runge_kutta_stepper< state_type , 1 > euler_stepper;
+typedef runge_kutta_stepper< state_type , 2 > midpoint_stepper;
const double sigma = 10.0;
@@ -37,9 +38,17 @@
vector< double > c( 1 , 0.0 );
euler_stepper euler( a , b , c );
+ vector< vector< double > > a2( 1 , vector<double>( 1 , 0.5 ) );
+ vector< double > b2( 2 , 0.0 ); b2[1] = 0.5;
+ vector< double > c2( 1 , 0.0 ); c2[1] = 1.0;
+ midpoint_stepper midpoint( a2 , b2 , c2 );
+
state_type x = {{ 1.0 , 1.0 , 2.0 }};
+ state_type x2 = {{ 1.0 , 1.0 , 2.0 }};
double t = 0.0;
euler.do_step( lorenz , x , t , dt );
+ midpoint.do_step( lorenz , x2 , t , dt );
cout << x[0] << " , " << x[1] << " , " << x[2] << endl;
+ cout << x2[0] << " , " << x2[1] << " , " << x2[2] << endl;
}