$include_dir="/home/hyper-archives/boost-commit/include"; include("$include_dir/msg-header.inc") ?>
Subject: [Boost-commit] svn:boost r66984 - in sandbox/odeint/branches/karsten/libs/numeric/odeint: ideas/generic_stepper test
From: mario.mulansky_at_[hidden]
Date: 2010-12-03 09:50:12
Author: mariomulansky
Date: 2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
New Revision: 66984
URL: http://svn.boost.org/trac/boost/changeset/66984
Log:
more performance tests
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_fusion.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_mpl.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_odeint.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper_fast.hpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile                 |    12 ++++++++++++                            
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp      |    34 +++++++++++++++++++---------------      
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp |    10 ++--------                              
   sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp               |    12 ++++--------                            
   4 files changed, 37 insertions(+), 31 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -27,4 +27,16 @@
     
 exe performance
     : performance.cpp
+    ;
+    
+exe performance_fusion
+    : performance_fusion.cpp runge_kutta_stepper_fast.hpp
+    ;
+    
+exe performance_mpl
+    : performance_mpl.cpp
+    ;
+    
+exe performance_odeint
+    : performance_odeint.cpp
     ;
\ No newline at end of file
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -15,28 +15,29 @@
 #include <boost/numeric/odeint/algebra/standard_algebra.hpp>
 #include <boost/numeric/odeint/algebra/standard_operations.hpp>
 
-template< typename state_type , size_t n >
+template< size_t n >
 struct fusion_algebra
 {
 
         typedef boost::numeric::odeint::standard_algebra std_algebra;
         typedef boost::numeric::odeint::standard_operations< double > std_op;
 
-
-	static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , n > &a ,
+	template< class state_type >
+	inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , n > &a ,
                             const state_type *k_vector , const double dt )
         {}
 
 };
 
-template< typename state_type>
-struct fusion_algebra< state_type , 1 >
+template<>
+struct fusion_algebra< 1 >
 {
 
         typedef boost::numeric::odeint::standard_algebra std_algebra;
         typedef boost::numeric::odeint::standard_operations< double > std_op;
 
-	static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 1 > &a ,
+	template< class state_type >
+	inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 1 > &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 ) );
@@ -45,14 +46,15 @@
 };
 
 
-template< typename state_type>
-struct fusion_algebra< state_type , 2 >
+template<>
+struct fusion_algebra< 2 >
 {
 
         typedef boost::numeric::odeint::standard_algebra std_algebra;
         typedef boost::numeric::odeint::standard_operations< double > std_op;
 
-	static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 2 > &a ,
+	template< class state_type >
+	inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 2 > &a ,
                             const state_type *k_vector , const double dt )
         {
                 std_algebra::for_each4( x_tmp , x ,  k_vector[0] , k_vector[1] ,
@@ -62,14 +64,15 @@
 };
 
 
-template< typename state_type>
-struct fusion_algebra< state_type , 3 >
+template<>
+struct fusion_algebra< 3 >
 {
 
     typedef boost::numeric::odeint::standard_algebra std_algebra;
     typedef boost::numeric::odeint::standard_operations< double > std_op;
 
-    static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 3 > &a ,
+    template< class state_type >
+    inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 3 > &a ,
                             const state_type *k_vector , const double dt )
     {
         std_algebra::for_each5( x_tmp , x ,  k_vector[0] , k_vector[1] , k_vector[2] ,
@@ -80,14 +83,15 @@
 
 
 
-template< typename state_type>
-struct fusion_algebra< state_type , 4 >
+template<>
+struct fusion_algebra< 4 >
 {
 
     typedef boost::numeric::odeint::standard_algebra std_algebra;
     typedef boost::numeric::odeint::standard_operations< double > std_op;
 
-    static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 4 > &a ,
+    template< class state_type >
+    inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 4 > &a ,
                             const state_type *k_vector , const double dt )
     {
         std_algebra::for_each6( x_tmp , x ,  k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] ,
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_fusion.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_fusion.cpp	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -0,0 +1,109 @@
+/*
+ * performance.cpp
+ *
+ *  Created on: Dec 1, 2010
+ *      Author: mario
+ */
+
+/*
+ * butcher_test.cpp
+ *
+ *  Created on: Nov 5, 2010
+ *      Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <tr1/array>
+
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+#include "runge_kutta_stepper_fast.hpp"
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::accumulators;
+
+typedef accumulator_set<
+    double , stats< tag::mean , tag::variance >
+    > accumulator_type;
+
+ostream& operator<<( ostream& out , accumulator_type &acc )
+{
+    out << boost::accumulators::mean( acc ) << tab;
+//    out << boost::accumulators::variance( acc ) << tab;
+    return out;
+}
+
+typedef boost::timer timer_type;
+
+
+typedef std::tr1::array< double , 3 > state_type;
+typedef runge_kutta_stepper< state_type , 4 > rk4_fusion_type;
+
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+    const double sigma = 10.0;
+    const double R = 28.0;
+    const double b = 8.0 / 3.0;
+    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];
+}
+
+
+
+
+int main( int argc , char **argv )
+{
+    typedef rk4_fusion_type::coef_a_type coef_a_type;
+    typedef rk4_fusion_type::coef_b_type coef_b_type;
+    typedef rk4_fusion_type::coef_c_type coef_c_type;
+
+    const boost::array< double , 1 > a1 = {{ 0.5 }};
+    const boost::array< double , 2 > a2 = {{ 0.0 , 0.5 }};
+    const boost::array< double , 3 > a3 = {{ 0.0 , 0.0 , 1.0 }};
+
+    const coef_a_type a = fusion::make_vector( a1 , a2 , a3 );
+    const coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
+    const coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
+
+    rk4_fusion_type rk4_fusion( a , b , c );
+
+    const size_t num_of_steps = 20000000;
+    const size_t dt = 0.01;
+
+    accumulator_type acc;
+    timer_type timer;
+
+    srand48( 12312354 );
+
+    while( true )
+    {
+        state_type x = {{ 10.0 * drand48() , 10.0 * drand48() , 10.0 * drand48() }};
+        double t = 0.0;
+
+        timer.restart();
+        for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
+            rk4_fusion.do_step( lorenz , x , t , dt );
+        acc( timer.elapsed() );
+
+        clog.precision( 3 );
+        clog.width( 5 );
+        clog << acc << " " << x[0] << endl;
+    }
+
+
+
+    return 0;
+}
+
+/*
+ * Compile with :
+ * g++ -O3 -I$BOOST_ROOT -I$HOME/boost/chrono -I$ODEINT_ROOT butcher_performance.cpp
+ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_mpl.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_mpl.cpp	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -0,0 +1,97 @@
+/*
+ * performance.cpp
+ *
+ *  Created on: Dec 1, 2010
+ *      Author: mario
+ */
+
+/*
+ * butcher_test.cpp
+ *
+ *  Created on: Nov 5, 2010
+ *      Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <tr1/array>
+
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+#include "predefined_steppers.hpp"
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::accumulators;
+
+typedef accumulator_set<
+    double , stats< tag::mean , tag::variance >
+    > accumulator_type;
+
+ostream& operator<<( ostream& out , accumulator_type &acc )
+{
+    out << boost::accumulators::mean( acc ) << tab;
+//    out << boost::accumulators::variance( acc ) << tab;
+    return out;
+}
+
+typedef boost::timer timer_type;
+
+
+typedef std::tr1::array< double , 3 > state_type;
+typedef mpl_rk4_stepper< state_type > rk4_mpl_type;
+
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+    const double sigma = 10.0;
+    const double R = 28.0;
+    const double b = 8.0 / 3.0;
+    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];
+}
+
+
+
+
+int main( int argc , char **argv )
+{
+    rk4_mpl_type rk4_mpl;
+
+    const size_t num_of_steps = 20000000;
+    const size_t dt = 0.01;
+
+    accumulator_type acc;
+    timer_type timer;
+
+    srand48( 12312354 );
+
+    while( true )
+    {
+        state_type x = {{ 10.0 * drand48() , 10.0 * drand48() , 10.0 * drand48() }};
+        double t = 0.0;
+
+        timer.restart();
+        for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
+            rk4_mpl.do_step( lorenz , x , t , dt );
+        acc( timer.elapsed() );
+
+        clog.precision( 3 );
+        clog.width( 5 );
+        clog << acc << " " << x[0] << endl;
+    }
+
+
+
+    return 0;
+}
+
+/*
+ * Compile with :
+ * g++ -O3 -I$BOOST_ROOT -I$HOME/boost/chrono -I$ODEINT_ROOT butcher_performance.cpp
+ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_odeint.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_odeint.cpp	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -0,0 +1,96 @@
+/*
+ * performance.cpp
+ *
+ *  Created on: Dec 1, 2010
+ *      Author: mario
+ */
+
+/*
+ * butcher_test.cpp
+ *
+ *  Created on: Nov 5, 2010
+ *      Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <tr1/array>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::accumulators;
+
+typedef accumulator_set<
+    double , stats< tag::mean , tag::variance >
+    > accumulator_type;
+
+ostream& operator<<( ostream& out , accumulator_type &acc )
+{
+    out << boost::accumulators::mean( acc ) << tab;
+//    out << boost::accumulators::variance( acc ) << tab;
+    return out;
+}
+
+typedef boost::timer timer_type;
+
+
+typedef std::tr1::array< double , 3 > state_type;
+typedef boost::numeric::odeint::explicit_rk4< state_type > rk4_odeint_type;
+
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+    const double sigma = 10.0;
+    const double R = 28.0;
+    const double b = 8.0 / 3.0;
+    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];
+}
+
+
+
+
+int main( int argc , char **argv )
+{
+    rk4_odeint_type rk4_odeint;
+
+    const size_t num_of_steps = 20000000;
+    const size_t dt = 0.01;
+
+    accumulator_type acc;
+    timer_type timer;
+
+    srand48( 12312354 );
+
+    while( true )
+    {
+        state_type x = {{ 10.0 * drand48() , 10.0 * drand48() , 10.0 * drand48() }};
+        double t = 0.0;
+
+        timer.restart();
+        for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
+            rk4_odeint.do_step( lorenz , x , t , dt );
+        acc( timer.elapsed() );
+
+        clog.precision( 3 );
+        clog.width( 5 );
+        clog << acc << " " << x[0] << endl;
+    }
+
+
+
+    return 0;
+}
+
+/*
+ * Compile with :
+ * g++ -O3 -I$BOOST_ROOT -I$HOME/boost/chrono -I$ODEINT_ROOT butcher_performance.cpp
+ */
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -1,13 +1,8 @@
 #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 <boost/mpl/size_t.hpp>
 
 #include <boost/fusion/container.hpp>
@@ -16,7 +11,6 @@
 #include <boost/fusion/adapted.hpp>
 
 
-#include <vector>
 #include <algorithm>
 #include <iostream>
 #include <string>
@@ -229,7 +223,7 @@
             else
                 system( x_tmp , k_vector[stage_number-1] , t + c * dt );
 
-            fusion_algebra<state_type , stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+            fusion_algebra<stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
         }
 
 
@@ -244,7 +238,7 @@
             else
                 system( x_tmp , k_vector[stage_number-1] , t + c * dt );
 
-            fusion_algebra<state_type , stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+            fusion_algebra<stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
         }
 
 
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper_fast.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper_fast.hpp	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -0,0 +1,201 @@
+#include <boost/mpl/vector.hpp>
+#include <boost/mpl/push_back.hpp>
+#include <boost/mpl/for_each.hpp>
+#include <boost/mpl/range_c.hpp>
+#include <boost/mpl/copy.hpp>
+#include <boost/mpl/size_t.hpp>
+
+#include <boost/fusion/container.hpp>
+#include <boost/fusion/algorithm.hpp>
+#include <boost/fusion/sequence.hpp>
+#include <boost/fusion/adapted.hpp>
+
+
+#include <algorithm>
+#include <iostream>
+#include <string>
+
+#include <boost/array.hpp>
+#include <typeinfo>
+
+
+#include "fusion_algebra.hpp"
+
+namespace mpl = boost::mpl;
+namespace fusion = boost::fusion;
+
+using namespace std;
+
+
+struct first_stage {};
+struct intermediate_stage {};
+struct last_stage {};
+
+
+template< class T , class Constant >
+struct array_wrapper
+{
+    typedef typename boost::array< T , Constant::value > type;
+};
+
+template< class T , class Constant , class StageCategory >
+struct stage_fusion_wrapper
+{
+    typedef typename fusion::vector< size_t , T , boost::array< T , Constant::value > , StageCategory > type;
+};
+
+template< class StateType , size_t stage_count >
+class runge_kutta_stepper
+{
+
+public:
+
+    typedef StateType state_type;
+
+    typedef mpl::range_c< size_t , 1 , stage_count > stage_indices;
+
+    typedef typename fusion::result_of::as_vector
+    <
+        typename mpl::copy
+        <
+            stage_indices ,
+            mpl::inserter
+            <
+                mpl::vector0< > ,
+                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
+            >
+        >::type
+    >::type coef_a_type;
+
+    typedef boost::array< double , stage_count > coef_b_type;
+    typedef boost::array< double , stage_count > coef_c_type;
+
+    typedef typename fusion::result_of::as_vector
+    <
+        typename mpl::push_back
+        <
+			typename mpl::copy
+            <
+				stage_indices,
+                mpl::inserter
+                <
+                    mpl::vector0<> ,
+                    mpl::push_back< mpl::_1 , stage_fusion_wrapper< double , mpl::_2 , intermediate_stage > >
+                >
+            >::type ,
+            typename stage_fusion_wrapper< double , mpl::size_t< stage_count > , last_stage >::type
+        >::type
+    >::type stage_vector_base;
+
+
+    struct stage_vector : public stage_vector_base
+    {
+    	struct do_insertion
+    	{
+    		stage_vector_base &m_base;
+    		const coef_a_type &m_a;
+    		const coef_c_type &m_c;
+
+    		do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
+    		: m_base( base ) , m_a( a ) , m_c( c ) { }
+
+    		template< class Index >
+    		void operator()( Index ) const
+    		{
+    			fusion::at_c< 0 >( fusion::at< Index >( m_base ) ) = Index::value;
+    			fusion::at_c< 1 >( fusion::at< Index >( m_base ) ) = m_c[ Index::value ];
+    			fusion::at_c< 2 >( fusion::at< Index >( m_base ) ) = fusion::at< Index >( m_a );
+    		}
+    	};
+
+    	stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
+    	{
+    		typedef mpl::range_c< size_t , 0 , stage_count - 1 > indices;
+    		mpl::for_each< indices >( do_insertion( *this , a , c ) );
+			fusion::at_c< 0 >( fusion::at_c< stage_count - 1 >( *this ) ) = stage_count - 1 ;
+			fusion::at_c< 1 >( fusion::at_c< stage_count - 1 >( *this ) ) = c[ stage_count - 1 ];
+			fusion::at_c< 2 >( fusion::at_c< stage_count - 1 >( *this ) ) = b;
+    	}
+    };
+
+
+
+    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< typename T , size_t stage_number >
+        inline void operator()( fusion::vector< size_t , T , boost::array< T , stage_number > , intermediate_stage > const &stage ) const
+        //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
+        {
+            double c = fusion::at_c< 1 >( stage );
+
+            if( stage_number == 1 )
+                system( x , k_vector[stage_number-1] , t + c * dt );
+            else
+                system( x_tmp , k_vector[stage_number-1] , t + c * dt );
+
+            fusion_algebra<stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+        }
+
+
+        template< typename T , size_t stage_number >
+        inline void operator()( fusion::vector< size_t , T , boost::array< T , stage_number > , last_stage > const &stage ) const
+        //void operator()( typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , last_stage >::type const &stage ) const
+        {
+            double c = fusion::at_c< 1 >( stage );
+
+            if( stage_number == 1 )
+                system( x , k_vector[stage_number-1] , t + c * dt );
+            else
+                system( x_tmp , k_vector[stage_number-1] , t + c * dt );
+
+            fusion_algebra<stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+        }
+
+
+    };
+
+
+
+
+public:
+
+    runge_kutta_stepper( const coef_a_type &a ,
+                            const coef_b_type &b ,
+                            const coef_c_type &c )
+    : m_a( a ) , m_b( b ) , m_c( c ) ,
+      m_stages( a , b , c )
+
+    { }
+
+
+    template< class System >
+    inline 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 ) );
+    }
+
+
+
+
+private:
+
+    const coef_a_type m_a;
+    const coef_b_type m_b;
+    const coef_c_type m_c;
+    const stage_vector m_stages;
+    state_type m_x_tmp;
+    state_type m_k_vector[stage_count];
+};
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/check_stepper_concepts.cpp	2010-12-03 09:50:08 EST (Fri, 03 Dec 2010)
@@ -36,6 +36,8 @@
 #include <boost/numeric/odeint.hpp>
 #include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
 
+#include <boost/type_traits/add_reference.hpp>
+
 #include "vector_space_1d.hpp"
 
 
@@ -55,15 +57,13 @@
 
 template< class State > struct algebra_dispatcher { typedef standard_algebra type; };
 template<> struct algebra_dispatcher< vector_space_type > { typedef vector_space_algebra type; };
+template<> struct algebra_dispatcher< double > { typedef vector_space_algebra type; };
 
 
 void constant_system_vector( const vector_type &x , vector_type &dxdt , double t ) { dxdt[0] = 1.0; }
 void constant_system_vector_space( const vector_space_type &x , vector_space_type &dxdt , double t ) { dxdt.m_x = 1.0; }
 void constant_system_array( const array_type &x , array_type &dxdt , double t ) { dxdt[0] = 1.0; }
 
-
-
-
 const double eps = 1.0e-14;
 
 template< class Stepper , class System >
@@ -85,7 +85,7 @@
     typedef typename stepper_type::order_type order_type;
     typedef typename stepper_type::time_type time_type;
 
-    stepper.do_step( system , x , 0.0 , 0.1 , xerr );
+    stepper.do_step( system , typename boost::add_reference< container_type>::type( x ), 0.0 , 0.1 ,  typename boost::add_reference< container_type>::type( xerr ) );
 }
 
 template< class Stepper , class System >
@@ -146,8 +146,6 @@
 };
 
 
-
-
 template< class State > class stepper_methods : public mpl::vector<
         explicit_euler< State , double , typename algebra_dispatcher< State >::type > ,
         explicit_rk4< State , double , typename algebra_dispatcher< State >::type > ,
@@ -243,7 +241,6 @@
 };
 
 
-
 template< class State > class error_stepper_methods : public mpl::vector<
         explicit_error_rk54_ck< State , double , typename algebra_dispatcher< State >::type > ,
         explicit_error_dopri5< State , double , typename algebra_dispatcher< State >::type >
@@ -340,7 +337,6 @@
         }
 };
 
-
 template< class State > class controlled_stepper_methods : public mpl::vector<
         controlled_error_stepper< explicit_error_rk54_ck< State , double , typename algebra_dispatcher< State >::type > > ,
         controlled_error_stepper< explicit_error_dopri5< State , double , typename algebra_dispatcher< State >::type > >