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

Subject: [Boost-commit] svn:boost r66564 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper
From: mario.mulansky_at_[hidden]
Date: 2010-11-13 18:38:12

Author: mariomulansky
Date: 2010-11-13 18:38:06 EST (Sat, 13 Nov 2010)
New Revision: 66564
URL: http://svn.boost.org/trac/boost/changeset/66564

Log:
it conpiles, but quite messy
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/test.cpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/generic_stepper.hpp |   187 +++++++++++++++++++++++++++++---------- 
   1 files changed, 139 insertions(+), 48 deletions(-)

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile	2010-11-13 18:38:06 EST (Sat, 13 Nov 2010)
@@ -0,0 +1,20 @@
+# (C) Copyright 2010 : Karsten Ahnert, 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)
+
+import os ;
+import modules ;
+import path ; 
+
+path-constant HOME : [ os.environ HOME ] ;
+
+project
+    : requirements
+      <define>BOOST_ALL_NO_LIB=1
+      <include>../../../../..
+      <include>$(BOOST_ROOT)
+    ;
+
+exe test 
+	: test.cpp
+	;

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp	2010-11-13 18:38:06 EST (Sat, 13 Nov 2010)
@@ -0,0 +1,45 @@
+/*
+ * algebra.hpp
+ *
+ *  Created on: Nov 13, 2010
+ *      Author: mario
+ */
+
+#ifndef ALGEBRA_HPP_
+#define ALGEBRA_HPP_
+
+#include <vector>
+
+#include <boost/numeric/odeint/algebra/standard_algebra.hpp>
+#include <boost/numeric/odeint/algebra/standard_operations.hpp>
+
+
+template< typename state_type , size_t n >
+struct algebra
+{
+
+	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 )
+	{}
+
+};
+
+template< typename state_type>
+struct algebra< state_type , 1 >
+{
+
+	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_each3( x_tmp , x ,  k_vector[0] , std_op::scale_sum2( 1.0 , a[0]*dt ) );
+	}
+
+};
+
+
+#endif /* ALGEBRA_HPP_ */

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/generic_stepper.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/generic_stepper.hpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/generic_stepper.hpp	2010-11-13 18:38:06 EST (Sat, 13 Nov 2010)
@@ -9,6 +9,37 @@
 #define GENERIC_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 "algebra.hpp"
+
+namespace mpl = boost::mpl;
+
+using namespace std;
+
+
+/*
+ * c_1 |
+ * c_2 | a_{2,1}
+ * c_3 | a_{3,1} a_{3,2}
+ * ...
+ * c_s | a_{s,1} a_{s,2} ... a_{s,s-1}
+ * -----------------------------------
+ *     | b_1     b_2         b_{s-1}    b_s
+ */
+
+
 /* Runge Kutta Stages */
 
 
@@ -26,14 +57,15 @@
 
     void init( parameter_array_type2D &a , parameter_array_type1D &c )
     {
-        m_a_row( a[stage-1] );
-        m_c( c[stage] );
+        m_a_row = a[stage-1];
+        m_c = c[stage-1];
     }
 
-    operator() ( state_type &x_tmp , const state_type &x , const state_type *k_vector , double t , const double dt )
+    template< typename System >
+    void operator() ( System &system , state_type &x_tmp , const state_type &x , const state_type *k_vector , double t , const double dt )
     {
-        system( x_tmp , k_vector[stage] , t + m_c * dt );
-        algebra::foreach< stage >( x_tmp , x , m_a_row , k_vector , dt);
+        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:
@@ -44,8 +76,8 @@
 };
 
 
-template<>
-class stepper_stage< 0 >
+template< typename State >
+class stepper_stage< State , 1 >
 {
 
 public:
@@ -58,17 +90,20 @@
 
     void init( parameter_array_type2D &a , parameter_array_type1D &c )
     {
-        m_c( c[0] );
+    	m_a_row = a[0];
+        m_c = c[0];
     }
 
-    operator() ( state_type &x_tmp , const state_type &x , const state_type *k_vector , double t , const double dt )
+    template< typename System >
+    void operator() ( System &system , state_type &x_tmp , const state_type &x , const state_type *k_vector , double t , const double dt )
     {
-        system( x , k_vector[0] , t + c * dt );
-        algebra::foreach< stage >( x_tmp , x , m_a_row , k_vector , dt);
+        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;
 
 };
@@ -86,82 +121,138 @@
     typedef vector< double > parameter_array_type1D;
 
 
-    void init( parameter_array_type2D &a , parameter_array_type1D &c )
+    void init( const parameter_array_type1D &b , const parameter_array_type1D &c )
     {
-        m_a_row( a[stage-1] );
-        m_c( c[stage] );
+        m_b = b;
+        m_c = c[stage-1];
     }
 
-    operator() ( state_type &x_tmp , const state_type &x , const state_type *k_vector , double t , const double dt )
+    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] , t + m_c * dt );
-        algebra::foreach< stage >( x , x , m_a_row , k_vector , 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_a_row;
+    vector< value_type > m_b;
     value_type m_c;
 
 };
 
 
-/* Runge Kutta Stepper - consisting of several stages */
-
-template< size_t stage_count >
-class runge_kutta_stepper
+template< typename State >
+class stepper_last_stage< State , 1 >
 {
-    struct initialize
-    {
-        typedef double value_type;
-        typedef vector< vector< double > > parameter_array_type2D;
-        typedef vector< double > parameter_array_type1D;
 
-        parameter_array_type2D &a;
-        parameter_array_type1D &c;
+public:
 
-        initialize( parameter_array_type2D &_a , parameter_array_type1D &_c ) : a( _a ) , c( _c ) { }
+    typedef double value_type;
+    typedef State state_type;
+    typedef vector< vector< double > > parameter_array_type2D;
+    typedef vector< double > parameter_array_type1D;
 
 
-        template< typename Stage >
-        operator() ( Stage &s ) { s.init( a , c ); }
-    };
+    void init( const parameter_array_type1D &b , const parameter_array_type1D &c )
+    {
+        m_b = b;
+        m_c = c[0];
+    }
 
     template< typename System >
+    void do_step( 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;
+
+};
+
+
+/* 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
     {
-        calculate_stage() {}
+    	System &system;
+    	state_type &x , &x_tmp;
+    	state_type *k_vector;
+    	const double t;
+    	const double dt;
+    	const parameter_array_type2D &m_a;
+    	const parameter_array_type1D &m_b;
+    	const parameter_array_type1D &m_c;
+
+    	calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector ,
+    						const double _t , const double _dt ,
+    						const parameter_array_type2D &a ,
+    						const parameter_array_type1D &b ,
+    						const parameter_array_type1D &c )
+    	: system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt ),
+    	  m_a( a ) , m_b( b ) , m_c( c )
+    	{}
+
+    	template< class Stage >
+    	void operator()( Stage &s )
+    	{
+    		s.init( m_a , m_b , m_c );
+    		s( x_tmp , x , k_vector , t , dt );
+    	}
+
     };
 
 public:
 
-    typedef mpl::range_c< size_t , 0 , stage_count-1 > indices;
-    typedef mpl::copy // create mpl::vector< stepper_stage< 0 >, stepper_stage< 1 > , .. stepper_stage< stage_count-1 > >
+    typedef mpl::range_c< size_t , 1 , stage_count > indices;
+    typedef typename mpl::copy // create mpl::vector< stepper_stage< 0 >, stepper_stage< 1 > , .. stepper_stage< stage_count-1 > >
     <
       indices ,
       mpl::inserter
       <
         mpl::vector0<> ,
-        mpl::insert_range
-        <
-          mpl::_1 ,
-          mpl::end< mpl::_1 > ,
-          stage< mpl::_2 >
-        >
+        stepper_stage< State , mpl::_2::value >
       >
-    >::type stages;
-    typedef stepper_last_stage< stage_count > last_stage;
+    >::type stage_vector;
+    typedef stepper_last_stage< State , stage_count > last_stage_type;
 
-    runge_kutta_stepper( parameter_array_type2D &a , parameter_array_type1D &c ) : last_stage( a , c )
+    runge_kutta_stepper( const parameter_array_type2D &a ,
+							const parameter_array_type1D &b ,
+							const parameter_array_type1D &c )
+    	: m_a( a ) , m_b( b ) , m_c( c )
     {
-        mpl::for_each< stages >( initialize( a , c ) );
+    	last_stage.init( b , c );
     }
 
-    void do_step()
+    template< class System >
+    void do_step( System &system , state_type &x , double t , const double dt )
     {
-        mpl::for_each< stages >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
+        mpl::for_each< stage_vector >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ,
+        													 	   m_a , m_b , m_c ) );
+        last_stage.do_step( system , m_x_tmp , x , m_k_vector , t , dt );
     }
 
+private:
+
+	state_type m_k_vector[stage_count];
+	state_type m_x_tmp;
+	last_stage_type last_stage;
+	const parameter_array_type2D &m_a;
+	const parameter_array_type1D &m_b;
+	const parameter_array_type1D &m_c;
 };
 
 

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/test.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/test.cpp	2010-11-13 18:38:06 EST (Sat, 13 Nov 2010)
@@ -0,0 +1,45 @@
+/*
+ * test.cpp
+ *
+ *  Created on: Nov 13, 2010
+ *      Author: mario
+ */
+
+#include <iostream>
+#include <vector>
+#include <tr1/array>
+
+#include "generic_stepper.hpp"
+
+using namespace std;
+
+typedef tr1::array< double , 3 > state_type;
+typedef runge_kutta_stepper< state_type , 1 > euler_stepper;
+
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+    dxdt[0] = sigma * ( x[1] - x[0] );
+    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+    dxdt[2] = x[0]*x[1] - b * x[2];
+}
+
+const double dt = 0.001;
+
+int main( void )
+{
+	vector< vector <double> > a( 0, vector<double>( 0 ) );
+	vector< double > b( 1 , 1.0 );
+	vector< double > c( 1 , 0.0 );
+	euler_stepper euler( a , b , c );
+
+	state_type x = {{ 1.0 , 1.0 , 2.0 }};
+	double t = 0.0;
+	euler.do_step( lorenz , x , t , dt );
+
+	cout << x[0] << " , " << x[1] << " , " << x[2] << endl;
+}