Subject: [Boost-commit] svn:boost r71532 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta
From: mario.mulansky_at_[hidden]
Date: 2011-04-27 10:51:08
Author: mariomulansky
Date: 2011-04-27 10:51:04 EDT (Wed, 27 Apr 2011)
New Revision: 71532
URL: http://svn.boost.org/trac/boost/changeset/71532
Log:
+new
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp | 20 ++++++++++++--------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp | 1 -
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp | 7 ++++---
3 files changed, 16 insertions(+), 12 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp 2011-04-27 10:51:04 EDT (Wed, 27 Apr 2011)
@@ -10,18 +10,17 @@
#include <boost/array.hpp>
-#include <boost/numeric/odeint/algebra/range_algebra.hpp>
-#include <boost/numeric/odeint/algebra/default_operations.hpp>
template< size_t n >
struct fusion_algebra
{
- 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 T , size_t dim >
+ inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
+ const boost::array< double , n > &a ,
+ const boost::array< T , dim > *k_vector , const double dt )
{
- for( size_t i=0 ; i<x.size() ; ++i )
+ for( size_t i=0 ; i<dim ; ++i )
{
x_tmp[i] = x[i];
for( size_t j = 0 ; j<n ; ++j )
@@ -31,7 +30,13 @@
};
-/*
+
+
+
+/** hand-wise implementation for performance improvement for n = 1..4 **/
+
+/* !!!!!!! Actually, this is factor 3 slower with intel compiler, so we don'y use it !!!!!
+ *
template<>
struct fusion_algebra< 1 >
@@ -83,7 +88,6 @@
};
-
template<>
struct fusion_algebra< 4 >
{
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp 2011-04-27 10:51:04 EDT (Wed, 27 Apr 2011)
@@ -53,7 +53,6 @@
};
-
template< class StateType , size_t stage_count >
class explicit_rk
{
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp 2011-04-27 10:51:04 EDT (Wed, 27 Apr 2011)
@@ -0,0 +1,237 @@
+/*
+ * fusion_runge_kutta.hpp
+ *
+ * Created on: Apr 26, 2011
+ * Author: mario
+ */
+
+#ifndef FUSION_EXPLICIT_RK_HPP_
+#define FUSION_EXPLICIT_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/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 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 T , size_t n , class StageCategory >
+struct stage
+{
+ static const size_t stage_number = n - 1;
+ T m_c;
+ boost::array< T , n > m_a;
+ typedef StageCategory category;
+};
+
+template< class T , size_t n>
+struct stage< T , n , last_stage >
+{
+ static const size_t stage_number = n - 1;
+ T m_c;
+ boost::array< T , n > m_b;
+ typedef last_stage category;
+};
+
+
+
+template< class T , class Constant , class StageCategory >
+struct stage_wrapper
+{
+ typedef stage< T , Constant::value , StageCategory > type;
+};
+
+
+template< class StateType , size_t stage_count >
+class explicit_rk
+{
+
+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_wrapper< double , mpl::_2 , intermediate_stage > >
+ >
+ >::type ,
+ stage< double , stage_count , last_stage >
+ >::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 ):: = Index::value;
+ fusion::at< Index >( m_base ).m_c = m_c[ Index::value ];
+ fusion::at< Index >( m_base ).m_a = 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< stage_count - 1 >( *this ).m_c = c[ stage_count - 1 ];
+ fusion::at_c< stage_count - 1 >( *this ).m_b = 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 >
+ void operator()( stage< 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 = stage.m_c;
+
+ 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 , stage.m_a , k_vector , dt);
+ }
+
+
+ template< typename T , size_t stage_number >
+ void operator()( stage< 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 = stage.m_c;
+
+ 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 , stage.m_b , k_vector , dt);
+ }
+
+
+ };
+
+
+
+
+public:
+
+ explicit_rk( 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 >
+ 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];
+};
+
+
+#endif /* FUSION_EXPLICIT_RK_HPP_ */
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp 2011-04-27 10:51:04 EDT (Wed, 27 Apr 2011)
@@ -8,13 +8,14 @@
#include <iostream>
#include <fstream>
-#include <tr1/array>
+#include <boost/array.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
-#include "fusion_explicit_rk.hpp"
+//#include "fusion_explicit_rk.hpp"
+#include "fusion_explicit_rk_new.hpp"
#define tab "\t"
@@ -35,7 +36,7 @@
typedef boost::timer timer_type;
-typedef std::tr1::array< double , 3 > state_type;
+typedef boost::array< double , 3 > state_type;
typedef explicit_rk< state_type , 4 > rk4_fusion_type;