[Boost-commit] svn:boost r66656 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper

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;
 }