Subject: [Boost-commit] svn:boost r66481 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher
From: karsten.ahnert_at_[hidden]
Date: 2010-11-10 05:34:05
Author: karsten
Date: 2010-11-10 05:34:01 EST (Wed, 10 Nov 2010)
New Revision: 66481
URL: http://svn.boost.org/trac/boost/changeset/66481
Log:
Adding bjam support for the butcher idea
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/algebra.hpp | 119 +++++++++++++++++++++++++++++++++------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp | 2
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_performance.cpp | 6 +-
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp | 4 -
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp | 8 +-
5 files changed, 110 insertions(+), 29 deletions(-)
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/Jamfile 2010-11-10 05:34:01 EST (Wed, 10 Nov 2010)
@@ -0,0 +1,29 @@
+# (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)
+ <include>$(HOME)/boost/chrono
+ ;
+
+exe butcher_performance
+ : butcher_performance.cpp
+ ;
+
+exe butcher_lorenz
+ : butcher_lorenz.cpp
+ ;
+
+exe butcher_vdp
+ : butcher_vdp.cpp
+ ;
\ No newline at end of file
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-11-10 05:34:01 EST (Wed, 10 Nov 2010)
@@ -11,6 +11,9 @@
#include <boost/mpl/int.hpp>
#include <boost/mpl/at.hpp>
+#include <boost/numeric/odeint/algebra/standard_algebra.hpp>
+#include <boost/numeric/odeint/algebra/standard_operations.hpp>
+
#include "convert_value.hpp"
namespace mpl = boost::mpl;
@@ -18,42 +21,88 @@
template< class state_type , class coef_tye , class index >
struct algebra
{
- static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
{}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 0 > >
{
- static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
+ typedef typename 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 )
{
- 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];
- }
+ const static double a1 = dt * convert_value< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
+
+ iterator first1 = x_tmp.begin() , last1 = x_tmp.end();
+ const_iterator first2 = x.begin() , first3 = k_vector[0].begin();
+ while( first1 != last1 )
+ *first1++ = *first2++ + a1 * *first3++ ;
+
+
+// std_algebra::for_each3( x_tmp , x , k_vector[0] ,
+// std_op::scale_sum2( 1.0 , a1 ) );
+
+
+// for( size_t i=0 ; i<x.size() ; ++i )
+// {
+// x_tmp[i] = x[i] + a1 * k_vector[0][i];
+// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 1 > >
{
- static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
+ typedef typename 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 )
{
- 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];
- }
+ 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++;
+
+// 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];
+// }
}
};
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 2 > >
{
- static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
+ typedef typename 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();
+//
+// 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] +
@@ -61,15 +110,31 @@
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 > >
{
- static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
+ typedef typename 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();
+//
+// 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();
+// const_iterator first6 = k_vector[3].begin();
+// 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] +
@@ -84,8 +149,26 @@
template< class state_type , class coef_type >
struct algebra< state_type , coef_type , mpl::int_< 4 > >
{
- static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , double dt )
+ typedef typename boost::numeric::odeint::standard_algebra< state_type > std_algebra;
+ typedef typename 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();
+//
+// 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();
+// const_iterator first6 = k_vector[3].begin() , first7 = k_vector[4].begin();
+// 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] +
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_lorenz.cpp 2010-11-10 05:34:01 EST (Wed, 10 Nov 2010)
@@ -20,7 +20,7 @@
typedef std::tr1::array< double , 3 > state_type;
typedef mpl_rk4_stepper< state_type > stepper_type;
-typedef boost::numeric::odeint::stepper_rk4< state_type > stepper_std_type;
+typedef boost::numeric::odeint::explicit_rk4< state_type > stepper_std_type;
const double sigma = 10.0;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_performance.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_performance.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_performance.cpp 2010-11-10 05:34:01 EST (Wed, 10 Nov 2010)
@@ -37,8 +37,8 @@
typedef std::tr1::array< double , 3 > state_type;
-typedef mpl_euler_stepper< state_type > stepper_type;
-typedef boost::numeric::odeint::stepper_euler< state_type > stepper_std_type;
+typedef mpl_rk4_stepper< state_type > stepper_type;
+typedef boost::numeric::odeint::explicit_rk4< state_type > stepper_std_type;
void lorenz( const state_type &x , state_type &dxdt , double t )
@@ -59,7 +59,7 @@
stepper_type stepper;
stepper_std_type stepper_std;
- const size_t num_of_steps = 10000000;
+ const size_t num_of_steps = 20000000;
const size_t dt = 0.01;
accumulator_type acc1 , acc2;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/butcher_vdp.cpp 2010-11-10 05:34:01 EST (Wed, 10 Nov 2010)
@@ -20,7 +20,7 @@
typedef std::tr1::array< double , 2 > state_type;
typedef mpl_rk4_stepper< state_type > stepper_type;
-typedef boost::numeric::odeint::stepper_rk4< state_type > stepper_std_type;
+typedef boost::numeric::odeint::explicit_rk4< state_type > stepper_std_type;
void vdp( const state_type &x , state_type &dxdt , double t )
{
@@ -29,8 +29,6 @@
dxdt[1] = -x[0] + mu * ( 1.0 - x[0]*x[0] ) * x[1];
}
-typedef mpl_rk4_stepper< state_type > stepper_type;
-typedef boost::numeric::odeint::stepper_rk4< state_type > stepper_std_type;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/butcher/explicit_runge_kutta.hpp 2010-11-10 05:34:01 EST (Wed, 10 Nov 2010)
@@ -72,10 +72,10 @@
System &system;
state_type &x , &x_tmp;
state_type *k_vector;
- double t;
- double dt;
+ const double t;
+ const double dt;
- calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , double _t , 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 )
{}
@@ -116,7 +116,7 @@
}
template< class System >
- void do_step( System &system , state_type &x , double t , double dt )
+ void do_step( System &system , state_type &x , double t , const double dt )
{
mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
}