[Boost-commit] svn:boost r72385 - in sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper: . detail

Subject: [Boost-commit] svn:boost r72385 - in sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper: . detail
From: mario.mulansky_at_[hidden]
Date: 2011-06-04 07:57:03

Author: mariomulansky
Date: 2011-06-04 07:57:03 EDT (Sat, 04 Jun 2011)
New Revision: 72385
URL: http://svn.boost.org/trac/boost/changeset/72385

Log:
refactored generic steppers
now all tests pass nicely
Added:
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/generic_rk_algorithm.hpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp |    94 +++++++++++++---                        
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp       |   223 +++------------------------------------ 
   2 files changed, 97 insertions(+), 220 deletions(-)

Added: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/generic_rk_algorithm.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/detail/generic_rk_algorithm.hpp	2011-06-04 07:57:03 EDT (Sat, 04 Jun 2011)
@@ -0,0 +1,258 @@
+/*
+ * generic_rk_algorithm.hpp
+ *
+ *  Created on: Jun 3, 2011
+ *      Author: mario
+ */
+
+#ifndef GENERIC_RK_ALGORITHM_HPP_
+#define GENERIC_RK_ALGORITHM_HPP_
+
+#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/iterator.hpp>
+#include <boost/fusion/mpl.hpp>
+#include <boost/fusion/sequence.hpp>
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+#include <boost/numeric/odeint/stepper/detail/generic_rk_call_algebra.hpp>
+#include <boost/numeric/odeint/stepper/detail/generic_rk_operations.hpp>
+
+namespace mpl = boost::mpl;
+namespace fusion = boost::fusion;
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+namespace detail {
+
+template< class T , class Constant >
+struct array_wrapper
+{
+    typedef const typename boost::array< T , Constant::value > type;
+};
+
+template< class T , size_t i >
+struct stage
+{
+    T c;
+    boost::array< T , i > a;
+};
+
+
+template< class T , class Constant >
+struct stage_wrapper
+{
+    typedef stage< T , Constant::value > type;
+};
+
+
+template<
+    size_t StageCount,
+    class Value = double ,
+    class Algebra = range_algebra,
+    class Operations = default_operations
+    >
+class generic_rk_algorithm {
+
+public:
+    typedef mpl::range_c< size_t , 1 , StageCount > 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< Value , mpl::_2 > >
+            >
+        >::type
+    >::type coef_a_type;
+
+    typedef boost::array< double , StageCount > coef_b_type;
+    typedef boost::array< double , StageCount > 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_wrapper< Value , mpl::_2 > >
+                >
+            >::type ,
+            stage< Value , StageCount >
+        >::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< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , fusion::at< Index >( m_a ) );
+                fusion::at< Index >( m_base ).c  = m_c[ Index::value ];
+                fusion::at< Index >( m_base ).a = fusion::at< Index >( m_a );
+            }
+        };
+
+
+        struct do_insertion_from_stage
+        {
+            stage_vector_base &m_base;
+            const stage_vector_base &m_source;
+
+            do_insertion_from_stage( stage_vector_base &base, const stage_vector_base &source )
+                    : m_base(base) , m_source( source )
+            { }
+
+            template<class Index>
+            void operator()(Index) const {
+                //fusion::at< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , fusion::at< Index >( m_a ) );
+                fusion::at<Index>(m_base).c = fusion::at<Index>(m_source).c;
+                fusion::at<Index>(m_base).a = fusion::at<Index>(m_source).a;
+            }
+        };
+
+        struct print_butcher
+        {
+            const stage_vector_base &m_base;
+            std::ostream &m_os;
+
+            print_butcher( const stage_vector_base &base , std::ostream &os )
+                : m_base( base ) , m_os( os )
+            { }
+
+            template<class Index>
+            void operator()(Index) const {
+                m_os << fusion::at<Index>(m_base).c << " | ";
+                for( size_t i=0 ; i<Index::value ; ++i )
+                    m_os << fusion::at<Index>(m_base).a[i] << " ";
+                m_os << std::endl;
+            }
+        };
+
+
+        stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
+        {
+            typedef mpl::range_c< size_t , 0 , StageCount-1 > indices;
+            mpl::for_each< indices >( do_insertion( *this , a , c ) );
+            fusion::at_c< StageCount - 1 >( *this ).c = c[ StageCount - 1 ];
+            fusion::at_c< StageCount - 1 >( *this ).a = b;
+        }
+
+        stage_vector( const stage_vector &s )
+        {
+            typedef mpl::range_c< size_t , 0 , StageCount > indices;
+            mpl::for_each< indices >( do_insertion_from_stage( *this , s ) );
+        }
+
+        void print( std::ostream &os ) const
+        {
+            typedef mpl::range_c< size_t , 0 , StageCount > indices;
+            mpl::for_each< indices >( print_butcher( *this , os ) );
+        }
+    };
+
+
+
+    template< class System , class StateIn , class StateTemp , class DerivIn , class Deriv ,
+               class StateOut , class Time >
+    struct calculate_stage
+    {
+        System &system;
+        const StateIn &x;
+        StateTemp &x_tmp;
+        StateOut &x_out;
+        const DerivIn &dxdt;
+        Deriv *F;
+        const Time t;
+        const Time dt;
+
+        calculate_stage( System &_system , const StateIn &_x , const DerivIn &_dxdt , StateOut &_out ,
+            StateTemp &_x_tmp , Deriv *_F , const Time &_t , const Time &_dt )
+        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , x_out( _out) , dxdt( _dxdt ) , F( _F ) , t( _t ) , dt( _dt )
+        {}
+
+
+        template< typename T , size_t stage_number >
+        void inline operator()( stage< T , stage_number > const &stage ) const
+        //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
+        {
+            if( stage_number > 1 )
+            {
+                #ifdef BOOST_MSVC
+                #pragma warning( disable : 4307 34 )
+                #endif
+                system( x_tmp , F[stage_number-2] , t + stage.c * dt );
+                #ifdef BOOST_MSVC
+                #pragma warning( default : 4307 34 )
+                #endif
+            }
+            //std::cout << stage_number-2 << ", t': " << t + stage.c * dt << std::endl;
+
+            if( stage_number < StageCount )
+                detail::template generic_rk_call_algebra< stage_number , Algebra >()( x_tmp , x , dxdt , F ,
+                            detail::generic_rk_scale_sum< stage_number , Operations , Time >( stage.a , dt) );
+//                  algebra_type::template for_eachn<stage_number>( x_tmp , x , dxdt , F ,
+//                          typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
+            else
+                detail::template generic_rk_call_algebra< stage_number , Algebra >()( x_out , x , dxdt , F ,
+                            detail::generic_rk_scale_sum< stage_number , Operations , Time >( stage.a , dt) );
+//                algebra_type::template for_eachn<stage_number>( x_out , x , dxdt , F ,
+//                            typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
+        }
+
+    };
+
+    generic_rk_algorithm( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
+        : m_stages( a , b , c )
+    { }
+
+    template< class System , class StateIn , class DerivIn , class Time , class StateOut , class StateTemp , class Deriv >
+    void inline do_step( System system , const StateIn &in , const DerivIn &dxdt ,
+                    const Time &t , StateOut &out , const Time &dt ,
+                    StateTemp &x_tmp , Deriv F[StageCount-1] )
+    {
+        fusion::for_each( m_stages , calculate_stage<
+                System , StateIn , StateTemp , DerivIn , Deriv , StateOut , Time >
+            ( system , in , dxdt , out , x_tmp , F , t , dt ) );
+    }
+
+private:
+    stage_vector m_stages;
+};
+
+
+}
+}
+}
+}
+
+#endif /* GENERIC_RK_ALGORITHM_HPP_ */

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp	(original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp	2011-06-04 07:57:03 EDT (Sat, 04 Jun 2011)
@@ -8,12 +8,11 @@
 #ifndef EXPLICIT_ERROR_GENERIC_RK_HPP_
 #define EXPLICIT_ERROR_GENERIC_RK_HPP_
 
-#include <boost/numeric/odeint/stepper/explicit_generic_rk.hpp>
 #include <boost/numeric/odeint/stepper/base/explicit_stepper_and_error_stepper_base.hpp>
 
 #include <boost/numeric/odeint/algebra/default_operations.hpp>
 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
-#include <boost/numeric/odeint/algebra/default_operations.hpp>
+#include <boost/numeric/odeint/stepper/detail/generic_rk_algorithm.hpp>
 #include <boost/numeric/odeint/stepper/detail/generic_rk_call_algebra.hpp>
 #include <boost/numeric/odeint/stepper/detail/generic_rk_operations.hpp>
 
@@ -63,43 +62,81 @@
     typedef typename stepper_base_type::adjust_size_policy adjust_size_policy;
     typedef typename stepper_base_type::stepper_type stepper_type;
 
-    typedef explicit_generic_rk< StageCount , Order , State , Value , Deriv ,
-                                  Time , Algebra , Operations , AdjustSizePolicy > generic_rk_stepper_type;
+    typedef detail::generic_rk_algorithm< StageCount , Value , Algebra , Operations > rk_algorithm_type;
 
-    typedef typename generic_rk_stepper_type::coef_a_type coef_a_type;
-    typedef typename generic_rk_stepper_type::coef_b_type coef_b_type;
-    typedef typename generic_rk_stepper_type::coef_c_type coef_c_type;
+    typedef typename rk_algorithm_type::coef_a_type coef_a_type;
+    typedef typename rk_algorithm_type::coef_b_type coef_b_type;
+    typedef typename rk_algorithm_type::coef_c_type coef_c_type;
 
     static const size_t stage_count = StageCount;
 
+private:
+
+    void initialize( void )
+    {
+        boost::numeric::odeint::construct( m_x_tmp );
+        m_state_adjuster.register_state( 0 , m_x_tmp );
+        for( size_t i = 0 ; i < StageCount-1 ; ++i )
+        {
+            boost::numeric::odeint::construct( m_F[i] );
+            m_deriv_adjuster.register_state( i , m_F[i] );
+        }
+    }
+
+    void copy( const explicit_error_generic_rk &rk )
+    {
+        boost::numeric::odeint::copy( rk.m_x_tmp , m_x_tmp );
+        for( size_t i = 0 ; i < StageCount-1 ; ++i )
+        {
+            boost::numeric::odeint::copy( rk.m_F[i] , m_F[i] );
+        }
+    }
+
+public:
+
     // we use an explicit_generic_rk to do the normal rk step
     // and add a separate calculation of the error estimate afterwards
     explicit_error_generic_rk( const coef_a_type &a ,
                                   const coef_b_type &b ,
                                   const coef_b_type &b2 ,
                                   const coef_c_type &c )
-        : m_rk_stepper( a , b , c ) , m_b2( b2 )
-    { }
+        : m_rk_algorithm( a , b , c ) , m_b2( b2 ) , m_x_tmp()
+    {
+        initialize();
+    }
 
     explicit_error_generic_rk( const explicit_error_generic_rk &rk )
-        : m_rk_stepper( rk.m_rk_stepper ) , m_b2( rk.m_b2 )
-    { }
+        : stepper_base_type( rk )  , m_rk_algorithm( rk.m_rk_algorithm ) , m_b2( rk.m_b2 ) , m_x_tmp()
+    {
+        initialize();
+        copy( rk );
+    }
 
     explicit_error_generic_rk& operator=( const explicit_error_generic_rk &rk )
     {
-        //stepper_base_type::operator=( rk );
-        m_rk_stepper = rk.m_rk_stepper;
+        stepper_base_type::operator=( rk );
+        copy( rk );
         return *this;
     }
 
+    ~explicit_error_generic_rk( void )
+   {
+       boost::numeric::odeint::destruct( m_x_tmp );
+       for( size_t i = 0 ; i < StageCount-1 ; ++i )
+       {
+           boost::numeric::odeint::destruct( m_F[i] );
+       }
+   }
+
     template< class System , class StateIn , class DerivIn , class StateOut , class Err >
     void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
             const time_type &t , StateOut &out , const time_type &dt , Err &xerr )
     {
-        m_rk_stepper.do_step_impl( system , in , dxdt , t , out , dt );
+        // normal step
+        do_step_impl( system , in , dxdt , t , out , dt );
 
-        // additionally perform the error calculation
-        detail::template generic_rk_call_algebra< StageCount , algebra_type >()( xerr , dxdt , m_rk_stepper.m_F ,
+        // additionally, perform the error calculation
+        detail::template generic_rk_call_algebra< StageCount , algebra_type >()( xerr , dxdt , m_F ,
                     detail::generic_rk_scale_sum_err< StageCount , operations_type , time_type >( m_b2 , dt) );
     }
 
@@ -107,14 +144,35 @@
     void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt ,
             const time_type &t , StateOut &out , const time_type &dt )
     {
-        m_rk_stepper.do_step_impl( system , in , dxdt , t , out , dt );
+        typedef typename boost::unwrap_reference< System >::type unwrapped_system_type;
+        unwrapped_system_type &sys = system;
+
+        m_deriv_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+        m_state_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
+
+        // actual calculation done in generic_rk.hpp
+        m_rk_algorithm.do_step( sys , in , dxdt , t , out , dt , m_x_tmp , m_F );
+    }
+
+    template< class StateType >
+    void adjust_size( const StateType &x )
+    {
+        m_deriv_adjuster.adjust_size( x );
+        m_state_adjuster.adjust_size( x );
+        stepper_base_type::adjust_size( x );
     }
 
 private:
 
-    generic_rk_stepper_type m_rk_stepper;
+    rk_algorithm_type m_rk_algorithm;
     const coef_b_type m_b2;
 
+    size_adjuster< deriv_type , StageCount-1 > m_deriv_adjuster;
+    size_adjuster< state_type , 1 > m_state_adjuster;
+
+    state_type m_x_tmp;
+    deriv_type m_F[StageCount-1];
+
 };
 
 }

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp	(original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp	2011-06-04 07:57:03 EDT (Sat, 04 Jun 2011)
@@ -8,29 +8,14 @@
 #ifndef EXPLICIT_GENERIC_RK_HPP_
 #define EXPLICIT_GENERIC_RK_HPP_
 
-#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/iterator.hpp>
-#include <boost/fusion/mpl.hpp>
-#include <boost/fusion/sequence.hpp>
-
 #include <boost/array.hpp>
 
 #include <boost/ref.hpp>
 
 #include <boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp>
-#include <boost/numeric/odeint/stepper/detail/macros.hpp>
 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
 #include <boost/numeric/odeint/algebra/default_operations.hpp>
-#include <boost/numeric/odeint/stepper/detail/generic_rk_call_algebra.hpp>
-#include <boost/numeric/odeint/stepper/detail/generic_rk_operations.hpp>
+#include <boost/numeric/odeint/stepper/detail/generic_rk_algorithm.hpp>
 
 //#include "fusion_foreach_performance.hpp"
 
@@ -39,8 +24,12 @@
 namespace mpl = boost::mpl;
 namespace fusion = boost::fusion;
 
+namespace boost {
+namespace numeric {
+namespace odeint {
+
 //forward declarations
-namespace boost { namespace numeric { namespace odeint {
+
 template<
     size_t StageCount,
     size_t Order,
@@ -53,15 +42,8 @@
     class AdjustSizePolicy = adjust_size_initially_tag
     >
 class explicit_generic_rk;
-struct stage_vector;
-} } }
-
-#include <boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp>
-
-namespace boost {
-namespace numeric {
-namespace odeint {
 
+struct stage_vector;
 
 template< class T , class Constant >
 struct array_wrapper
@@ -124,8 +106,9 @@
 public:
 
         typedef explicit_stepper_base<
-explicit_generic_rk< StageCount , Order , State , Value , Deriv ,Time , Algebra , Operations , AdjustSizePolicy > ,
-Order , State , Value , Deriv , Time , Algebra , Operations , AdjustSizePolicy > stepper_base_type;
+	        explicit_generic_rk< StageCount , Order , State , Value , Deriv ,Time , Algebra , Operations , AdjustSizePolicy > ,
+	        Order , State , Value , Deriv , Time , Algebra ,
+	        Operations , AdjustSizePolicy > stepper_base_type;
 
         typedef typename stepper_base_type::state_type state_type;
         typedef typename stepper_base_type::value_type value_type;
@@ -136,170 +119,13 @@
         typedef typename stepper_base_type::adjust_size_policy adjust_size_policy;
         typedef typename stepper_base_type::stepper_type stepper_type;
 
+	typedef detail::generic_rk_algorithm< StageCount , Value , Algebra , Operations > rk_algorithm_type;
 
-	typedef mpl::range_c< size_t , 1 , StageCount > 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< Value , mpl::_2 > >
-            >
-        >::type
-    >::type coef_a_type;
-
-    typedef boost::array< double , StageCount > coef_b_type;
-    typedef boost::array< double , StageCount > 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_wrapper< Value , mpl::_2 > >
-                >
-            >::type ,
-            stage< Value , StageCount >
-        >::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< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , fusion::at< Index >( m_a ) );
-                fusion::at< Index >( m_base ).c  = m_c[ Index::value ];
-                fusion::at< Index >( m_base ).a = fusion::at< Index >( m_a );
-            }
-        };
-
-
-        struct do_insertion_from_stage
-        {
-            stage_vector_base &m_base;
-            const stage_vector_base &m_source;
-
-            do_insertion_from_stage( stage_vector_base &base, const stage_vector_base &source )
-                    : m_base(base) , m_source( source )
-            { }
-
-            template<class Index>
-            void operator()(Index) const {
-                //fusion::at< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , fusion::at< Index >( m_a ) );
-                fusion::at<Index>(m_base).c = fusion::at<Index>(m_source).c;
-                fusion::at<Index>(m_base).a = fusion::at<Index>(m_source).a;
-            }
-        };
-
-        struct print_butcher
-        {
-            const stage_vector_base &m_base;
-            std::ostream &m_os;
-
-            print_butcher( const stage_vector_base &base , std::ostream &os )
-                : m_base( base ) , m_os( os )
-            { }
-
-            template<class Index>
-            void operator()(Index) const {
-                m_os << fusion::at<Index>(m_base).c << " | ";
-                for( size_t i=0 ; i<Index::value ; ++i )
-                    m_os << fusion::at<Index>(m_base).a[i] << " ";
-                m_os << std::endl;
-            }
-        };
-
-
-        stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
-        {
-            typedef mpl::range_c< size_t , 0 , StageCount-1 > indices;
-            mpl::for_each< indices >( do_insertion( *this , a , c ) );
-            fusion::at_c< StageCount - 1 >( *this ).c = c[ StageCount - 1 ];
-            fusion::at_c< StageCount - 1 >( *this ).a = b;
-        }
-
-        stage_vector( const stage_vector &s )
-        {
-            typedef mpl::range_c< size_t , 0 , StageCount > indices;
-            mpl::for_each< indices >( do_insertion_from_stage( *this , s ) );
-        }
-
-        void print( std::ostream &os ) const
-        {
-            typedef mpl::range_c< size_t , 0 , StageCount > indices;
-            mpl::for_each< indices >( print_butcher( *this , os ) );
-        }
-    };
-
+	typedef typename rk_algorithm_type::coef_a_type coef_a_type;
+    typedef typename rk_algorithm_type::coef_b_type coef_b_type;
+    typedef typename rk_algorithm_type::coef_c_type coef_c_type;
 
-
-    template< class System , class StateIn , class DerivIn , class StateOut >
-    struct calculate_stage
-    {
-        System &system;
-        const StateIn &x;
-		state_type &x_tmp;
-		StateOut &x_out;
-		const DerivIn &dxdt;
-        state_type *F;
-        const double t;
-        const double dt;
-
-        calculate_stage( System &_system , const StateIn &_x , const DerivIn &_dxdt , StateOut &_out ,
-			state_type &_x_tmp , state_type *_F , const time_type &_t , const time_type &_dt )
-        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , x_out( _out) , dxdt( _dxdt ) , F( _F ) , t( _t ) , dt( _dt )
-        {}
-
-
-        template< typename T , size_t stage_number >
-        void inline operator()( stage< T , stage_number > const &stage ) const
-        //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
-        {
-            if( stage_number > 1 )
-            {
-				#ifdef BOOST_MSVC
-				#pragma warning( disable : 4307 34 )
-				#endif
-                system( x_tmp , F[stage_number-2] , t + stage.c * dt );
-				#ifdef BOOST_MSVC
-				#pragma warning( default : 4307 34 )
-				#endif
-            }
-			//std::cout << stage_number-2 << ", t': " << t + stage.c * dt << std::endl;
-
-			if( stage_number < StageCount )
-			    detail::template generic_rk_call_algebra< stage_number , algebra_type >()( x_tmp , x , dxdt , F ,
-                            detail::generic_rk_scale_sum< stage_number , operations_type , time_type >( stage.a , dt) );
-//				    algebra_type::template for_eachn<stage_number>( x_tmp , x , dxdt , F ,
-//				            typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
-			else
-			    detail::template generic_rk_call_algebra< stage_number , algebra_type >()( x_out , x , dxdt , F ,
-                            detail::generic_rk_scale_sum< stage_number , operations_type , time_type >( stage.a , dt) );
-//                algebra_type::template for_eachn<stage_number>( x_out , x , dxdt , F ,
-//                            typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
-        }
-
-    };
+    static const size_t stage_count = StageCount;
 
 
 private:
@@ -327,17 +153,15 @@
 
 public:
 
-    static const size_t stage_count = StageCount;
-
     explicit_generic_rk( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
-        : m_stages( a , b , c ) , m_x_tmp()
+        : m_rk_algorithm( a , b , c ) , m_x_tmp()
 
     {
         initialize();
     }
 
     explicit_generic_rk( const explicit_generic_rk &rk )
-        : stepper_base_type( rk ) , m_stages( rk.m_stages ) , m_x_tmp()
+        : stepper_base_type( rk ) , m_rk_algorithm( rk.m_rk_algorithm) , m_x_tmp()
     {
         initialize();
         copy( rk );
@@ -370,8 +194,8 @@
         m_deriv_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
         m_state_adjuster.adjust_size_by_policy( in , adjust_size_policy() );
 
-        fusion::for_each( m_stages , calculate_stage< unwrapped_system_type , StateIn , DerivIn , StateOut >
-			( sys , in , dxdt , out , m_x_tmp , m_F , t , dt ) );
+        // actual calculation done in generic_rk.hpp
+        m_rk_algorithm.do_step( sys , in , dxdt , t , out , dt , m_x_tmp , m_F );
     }
 
     template< class StateType >
@@ -384,19 +208,14 @@
 
     friend std::ostream& operator << <>( std::ostream &os , const explicit_generic_rk &rk );
 
-    template< size_t , size_t , size_t , size_t ,
-        typename , typename , typename , typename, typename , typename, typename >
-    friend class explicit_error_generic_rk;
-
 private:
 
-    const stage_vector m_stages;
-    state_type m_x_tmp;
+    rk_algorithm_type m_rk_algorithm;
 
     size_adjuster< deriv_type , StageCount-1 > m_deriv_adjuster;
     size_adjuster< state_type , 1 > m_state_adjuster;
 
-protected:
+    state_type m_x_tmp;
     deriv_type m_F[StageCount-1];
 
 };