|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r66949 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas: butcher generic_stepper
From: mario.mulansky_at_[hidden]
Date: 2010-12-01 10:52:22
Author: mariomulansky
Date: 2010-12-01 10:52:17 EST (Wed, 01 Dec 2010)
New Revision: 66949
URL: http://svn.boost.org/trac/boost/changeset/66949
Log:
rk fusion working
performance odeint, mpl , fusion
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp
- copied, changed from r66944, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance.cpp (contents, props changed)
Removed:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile | 2
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp | 113 ++++++++++++++++++++++-----------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile | 6 ++
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp | 65 +++++++++++++++++-----
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/rk_test.cpp | 47 ++++++++++++++--
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp | 66 ++++++++++++++++++++++
6 files changed, 225 insertions(+), 74 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile 2010-12-01 10:52:17 EST (Wed, 01 Dec 2010)
@@ -13,7 +13,7 @@
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
<include>$(BOOST_ROOT)
- <include>$(HOME)/boost/chrono
+ <include>$(HOME)/boost_sandbox/chrono
;
exe butcher_performance
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp 2010-12-01 10:52:17 EST (Wed, 01 Dec 2010)
@@ -28,8 +28,8 @@
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 0 > >
{
- typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
- typedef typename boost::numeric::odeint::standard_operations< double > std_op;
+ typedef boost::numeric::odeint::standard_algebra std_algebra;
+ typedef boost::numeric::odeint::standard_operations< double > std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
@@ -58,8 +58,8 @@
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 1 > >
{
- typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
- typedef typename boost::numeric::odeint::standard_operations< double > std_op;
+ typedef boost::numeric::odeint::standard_algebra std_algebra;
+ typedef boost::numeric::odeint::standard_operations< double > std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
@@ -69,10 +69,13 @@
const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
- iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
- const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin();
- while( first1 != last1 )
- *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++;
+// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
+// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin();
+// while( first1 != last1 )
+// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++;
+
+ std_algebra::for_each4( x_tmp , x , k_vector[0] , k_vector[1] ,
+ std_op::scale_sum3( 1.0 , a1 , a2 ) );
// for( size_t i=0 ; i<x.size() ; ++i )
// {
@@ -86,48 +89,51 @@
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 2 > >
{
- typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
- typedef typename boost::numeric::odeint::standard_operations< double > std_op;
+ typedef boost::numeric::odeint::standard_algebra std_algebra;
+ typedef boost::numeric::odeint::standard_operations< double > std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
-// const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
-// const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
-// const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
+ const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
+ const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
+ const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
//
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin() , first5 = k_vector[2].begin();
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++;
- for( size_t i=0 ; i<x.size() ; ++i )
- {
- x_tmp[i] = x[i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i];
- }
+ std_algebra::for_each5( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] ,
+ std_op::scale_sum4( 1.0 , a1 , a2 , a3 ) );
+
+// for( size_t i=0 ; i<x.size() ; ++i )
+// {
+// x_tmp[i] = x[i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i];
+// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 3 > >
{
- typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
- typedef typename boost::numeric::odeint::standard_operations< double > std_op;
+ typedef boost::numeric::odeint::standard_algebra std_algebra;
+ typedef boost::numeric::odeint::standard_operations< double > std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
-// const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
-// const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
-// const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
-// const static double a4 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value();
+ const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
+ const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
+ const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
+ const static double a4 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value();
//
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin() , first5 = k_vector[2].begin();
@@ -135,33 +141,36 @@
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++ + a4 * *first6++;
- for( size_t i=0 ; i<x.size() ; ++i )
- {
- x_tmp[i] = x[i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i];
- }
+ std_algebra::for_each6( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] ,
+ std_op::scale_sum5( 1.0 , a1 , a2 , a3 , a4 ) );
+
+// for( size_t i=0 ; i<x.size() ; ++i )
+// {
+// x_tmp[i] = x[i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i];
+// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 4 > >
{
- typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
- typedef typename boost::numeric::odeint::standard_operations< double > std_op;
+ typedef boost::numeric::odeint::standard_algebra std_algebra;
+ typedef boost::numeric::odeint::standard_operations< double > std_op;
typedef typename state_type::iterator iterator;
typedef typename state_type::const_iterator const_iterator;
static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{
-// const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
-// const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
-// const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
-// const static double a4 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value();
-// const static double a5 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 4 > >::type >::get_value();
+ const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
+ const static double a2 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();
+ const static double a3 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value();
+ const static double a4 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value();
+ const static double a5 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 4 > >::type >::get_value();
//
// iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
// const_iterator first2 = x.begin() , first3 = k_vector[0].begin() , first4 = k_vector[1].begin() , first5 = k_vector[2].begin();
@@ -169,15 +178,19 @@
// while( first1 != last1 )
// *first1++ = *first2++ + a1 * *first3++ + a2 * *first4++ + a3 * *first5++ + a4 * *first6++ + a5 * *first7++;
- for( size_t i=0 ; i<x.size() ; ++i )
- {
- x_tmp[i] = x[i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i] +
- dt * convert_value< typename mpl::at< coef_type , mpl::int_< 4 > >::type >::get_value() * k_vector[4][i];
- }
+ std_algebra::for_each7( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] , k_vector[3] , k_vector[4] ,
+ std_op::scale_sum6( 1.0 , a1 , a2 , a3 , a4 , a5 ) );
+
+
+// for( size_t i=0 ; i<x.size() ; ++i )
+// {
+// x_tmp[i] = x[i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value() * k_vector[0][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value() * k_vector[1][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 2 > >::type >::get_value() * k_vector[2][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 3 > >::type >::get_value() * k_vector[3][i] +
+// dt * convert_value< typename mpl::at< coef_type , mpl::int_< 4 > >::type >::get_value() * k_vector[4][i];
+// }
}
};
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-01 10:52:17 EST (Wed, 01 Dec 2010)
@@ -12,6 +12,8 @@
: requirements
<define>BOOST_ALL_NO_LIB=1
<include>../../../../..
+ <include>../butcher
+ <include>$(HOME)/boost_sandbox/chrono
<include>$(BOOST_ROOT)
;
@@ -21,4 +23,8 @@
exe rk_test
: rk_test.cpp
+ ;
+
+exe performance
+ : performance.cpp
;
\ No newline at end of file
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/algebra.hpp 2010-12-01 10:52:17 EST (Wed, 01 Dec 2010)
+++ (empty file)
@@ -1,65 +0,0 @@
-/*
- * algebra.hpp
- *
- * Created on: Nov 13, 2010
- * Author: mario
- */
-
-#ifndef ALGEBRA_HPP_
-#define ALGEBRA_HPP_
-
-#include <vector>
-
-#include <boost/array.hpp>
-
-#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 boost::array< double , n > &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 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 ) );
- }
-
-};
-
-
-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 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] , std_op::scale_sum3( 1.0 , a[0]*dt , a[1]*dt ) );
- }
-
-};
-
-
-#endif /* ALGEBRA_HPP_ */
Copied: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp (from r66944, /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/fusion_algebra.hpp 2010-12-01 10:52:17 EST (Wed, 01 Dec 2010)
@@ -1,12 +1,12 @@
/*
- * algebra.hpp
+ * fusion_algebra.hpp
*
* Created on: Nov 13, 2010
* Author: mario
*/
-#ifndef ALGEBRA_HPP_
-#define ALGEBRA_HPP_
+#ifndef FUSION_ALGEBRA_HPP_
+#define FUSION_ALGEBRA_HPP_
#include <vector>
@@ -15,13 +15,12 @@
#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
+struct fusion_algebra
{
- typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
- typedef typename boost::numeric::odeint::standard_operations< double > std_op;
+ 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 ,
@@ -31,11 +30,11 @@
};
template< typename state_type>
-struct algebra< state_type , 1 >
+struct fusion_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;
+ 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 ,
const state_type *k_vector , const double dt )
@@ -47,19 +46,55 @@
template< typename state_type>
-struct algebra< state_type , 2 >
+struct fusion_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;
+ 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 ,
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 ) );
+ 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_ */
+template< typename state_type>
+struct fusion_algebra< state_type , 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 ,
+ const state_type *k_vector , const double dt )
+ {
+ std_algebra::for_each5( x_tmp , x , k_vector[0] , k_vector[1] , k_vector[2] ,
+ std_op::scale_sum4( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt) );
+ }
+
+};
+
+
+
+template< typename state_type>
+struct fusion_algebra< state_type , 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 ,
+ 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] ,
+ std_op::scale_sum5( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt , a[3]*dt ) );
+ }
+
+};
+
+
+#endif /* FUSION_ALGEBRA_HPP_ */
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance.cpp 2010-12-01 10:52:17 EST (Wed, 01 Dec 2010)
@@ -0,0 +1,130 @@
+/*
+ * 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>
+
+#include "predefined_steppers.hpp"
+#include "runge_kutta_stepper.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;
+typedef mpl_rk4_stepper< state_type > rk4_mpl_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 )
+{
+ rk4_odeint_type rk4_odeint;
+ rk4_mpl_type rk4_mpl;
+
+ 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 acc1 , acc2 , acc3;
+ timer_type timer;
+
+ srand48( 12312354 );
+
+ while( true )
+ {
+ state_type x1 = {{ 10.0 * drand48() , 10.0 * drand48() , 10.0 * drand48() }};
+ state_type x2 = x1;
+ state_type x3 = x1;
+ double t1 = 0.0 ;
+ double t2 = 0.0;
+ double t3 = 0.0;
+
+ timer.restart();
+ for( size_t i=0 ; i<num_of_steps ; ++i,t1+=dt )
+ rk4_odeint.do_step( lorenz , x1 , t1 , dt );
+ acc1( timer.elapsed() );
+
+ timer.restart();
+ for( size_t i=0 ; i<num_of_steps ; ++i,t2+=dt )
+ rk4_mpl.do_step( lorenz , x2 , t2 , dt );
+ acc2( timer.elapsed() );
+
+ timer.restart();
+ for( size_t i=0 ; i<num_of_steps ; ++i,t3+=dt )
+ rk4_fusion.do_step( lorenz , x3 , t3 , dt );
+ acc3( timer.elapsed() );
+
+ clog.precision( 3 );
+ clog.width( 5 );
+ clog << acc1 << " " << acc2 << " " << acc3 << " " << x1[0] << " " << x2[0] << " " << x3[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/rk_test.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/rk_test.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/rk_test.cpp 2010-12-01 10:52:17 EST (Wed, 01 Dec 2010)
@@ -1,21 +1,56 @@
+#include <iostream>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
#include "runge_kutta_stepper.hpp"
+using namespace std;
+
+namespace odeint = boost::numeric::odeint;
+
+typedef boost::array< double , 3 > state_type;
+
+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];
+}
+
int main( int argc , char **argv )
{
- typedef runge_kutta_stepper< 4 > rk_type;
+ typedef runge_kutta_stepper< state_type , 4 > rk_type;
typedef rk_type::coef_a_type coef_a_type;
typedef rk_type::coef_b_type coef_b_type;
typedef rk_type::coef_c_type coef_c_type;
- const boost::array< double , 1 > a1 = {{ 5.1 }};
- const boost::array< double , 2 > a2 = {{ 5.2 , 5.3 }};
- const boost::array< double , 3 > a3 = {{ 5.4 , 5.5 , 5.6 }};
+ 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 = {{ 6.1 , 6.2 , 6.3 , 6.4 }};
- const coef_c_type c = {{ 1.1 , 1.2 , 1.3 , 1.4 }};
+ 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 }};
rk_type rk( a , b , c );
rk.print_vals();
+ state_type x = {{ 1.0 , 1.0 , 2.0 }};
+ double t = 0.0;
+ rk.do_step( lorenz , x , t , 0.1 );
+
+ cout.precision(16);
+
+ cout << x[0] << " , " << x[1] << " , " << x[2] << endl;
+
+
+ odeint::explicit_rk4< state_type > rk4;
+ state_type x2 = {{ 1.0 , 1.0 , 2.0 }};
+ rk4.do_step( lorenz , x2 , t , 0.1 );
+ cout << x2[0] << " , " << x2[1] << " , " << x2[2] << endl;
}
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-01 10:52:17 EST (Wed, 01 Dec 2010)
@@ -25,6 +25,8 @@
#include <typeinfo>
+#include "fusion_algebra.hpp"
+
namespace mpl = boost::mpl;
namespace fusion = boost::fusion;
@@ -64,7 +66,7 @@
template< class T , class Constant , class StageCategory >
struct stage_fusion_wrapper
{
- typedef typename fusion::vector< size_t , T , boost::array< T , Constant::value > , StateCategory > type;
+ typedef typename fusion::vector< size_t , T , boost::array< T , Constant::value > , StageCategory > type;
};
struct print_stage
@@ -125,11 +127,14 @@
};
-template< size_t stage_count >
+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
@@ -198,6 +203,54 @@
+ 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()( 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<state_type , stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+ }
+
+
+ template< typename T , size_t stage_number >
+ 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<state_type , stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+ }
+
+
+ };
+
+
public:
@@ -211,6 +264,13 @@
{ }
+ 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 ) );
+ }
+
+
void print_vals()
{
cout << "Typeof state_vector_base : " << endl;
@@ -237,4 +297,6 @@
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];
};
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk