|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r57436 - in sandbox/odeint: boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/examples libs/numeric/odeint/stuff/iterator_algebra
From: karsten.ahnert_at_[hidden]
Date: 2009-11-06 11:56:01
Author: karsten
Date: 2009-11-06 11:56:00 EST (Fri, 06 Nov 2009)
New Revision: 57436
URL: http://svn.boost.org/trac/boost/changeset/57436
Log:
small changes
Added:
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp (contents, props changed)
Text files modified:
sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp | 152 ++++++++++++++++++++++++---------------
sandbox/odeint/boost/numeric/odeint/euler.hpp | 11 +-
sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp | 19 +---
sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp | 2
sandbox/odeint/libs/numeric/odeint/stuff/iterator_algebra/increment_decrement.cc | 90 ++++++++++++++++++++---
5 files changed, 182 insertions(+), 92 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp 2009-11-06 11:56:00 EST (Fri, 06 Nov 2009)
@@ -22,82 +22,114 @@
namespace detail {
namespace it_algebra { // iterator algebra
- /* computes x_n += alpha * y_n for all n
- where x and y are given by the iterators x_start, x_end, y_start
- make sure that x and y have the same range
- */
- template< class Iterator1 ,
- class Iterator2 ,
- class T >
- void scale_and_add( Iterator1 x_start ,
- Iterator1 x_end ,
- Iterator2 y_start ,
- T alpha )
+
+ // computes y += alpha * x1
+ template <
+ class InOutIterator ,
+ class InputIterator ,
+ class T
+ >
+ void increment(
+ InOutIterator first1 ,
+ InOutIterator last1 ,
+ InputIterator first2 ,
+ T alpha
+ )
{
- while( x_start != x_end )
- (*x_start++) += alpha * (*y_start++);
+ while( first1 < last1 )
+ (*first1++) += alpha * (*first2++);
}
- /* computes out_n = x_n - y_n for all n
- where x,y, and out are given as iterators
- make sure that x, y and out have the same range
- */
- template< class InputIterator1 ,
- class InputIterator2 ,
- class OutputIterator >
- void substract_and_assign( InputIterator1 x_start ,
- InputIterator1 x_end ,
- InputIterator2 y_start ,
- OutputIterator out_start )
+ // computes y = x1 - x2
+ template <
+ class OutputIterator ,
+ class InputIterator1 ,
+ class InputIterator2
+ >
+ void assign_diff(
+ OutputIterator first1 ,
+ OutputIterator last1 ,
+ InputIterator1 first2 ,
+ InputIterator2 first3 )
{
- while( x_start != x_end )
- ( *out_start++ ) = ( *x_start++ ) - ( *y_start++ );
+ while( first1 < last1 )
+ (*first1++) = (*first2++) - (*first3++);
}
-
- /* computes out_n = x_n + alpha * y_n for all n
- where x,y and out are given as iterators.
- make sure, that x,y and out have the same range
- */
- template< class InputIterator1 ,
- class InputIterator2 ,
- class OutputIterator ,
- class T >
- void scale_and_add_and_assign( InputIterator1 x_start ,
- InputIterator1 x_end ,
- InputIterator2 y_start ,
- OutputIterator out_start ,
- T alpha )
+ // computes y = x1 + alpha * x2
+ template <
+ class OutputIterator ,
+ class InputIterator1 ,
+ class InputIterator2 ,
+ class T
+ >
+ void assign_sum(
+ OutputIterator first1 ,
+ OutputIterator last1 ,
+ InputIterator1 first2 ,
+ InputIterator2 first3 ,
+ T alpha )
{
- while( x_start != x_end )
- ( *out_start++) = ( *x_start++ ) + alpha * ( *y_start++ );
+ while( first1 < last1 )
+ (*first1++) = (*first2++) + alpha * (*first3++);
}
+ // computes y = alpha1 * ( x1 + x2 + alpha2*x3 )
+ template <
+ class OutputIterator ,
+ class InputIterator1 ,
+ class InputIterator2 ,
+ class InputIterator3 ,
+ class T
+ >
+ void increment_sum_sum(
+ OutputIterator first1 ,
+ OutputIterator last1 ,
+ InputIterator1 first2 ,
+ InputIterator2 first3 ,
+ InputIterator3 first4 ,
+ T alpha1 ,
+ T alpha2
+ )
+ {
+ while( first1 < last1 )
+ (*first1++) += alpha1 *
+ ( (*first2++) + (*first3++) + alpha2*(*first4++) );
+ }
- /* computes out_n = alpha2* ( x_n + y_n + alpha1*z_n )for all n
- where x,y,z and out are given as iterators.
- make sure, that x,y,z and out have the same range
- */
- template< class InputIterator1 ,
- class InputIterator2 ,
- class InputIterator3 ,
- class OutputIterator ,
- class T >
- void scale_and_add_and_add_and_assign( InputIterator1 x_start ,
- InputIterator1 x_end ,
- InputIterator2 y_start ,
- InputIterator3 z_start ,
- OutputIterator out_start ,
- T alpha1 ,
- T alpha2 )
+ // computes y = x1 + alpha * x2 ; x2 += x3
+ template<
+ class OutputIterator ,
+ class InputIterator1 ,
+ class InOutIterator ,
+ class InputIterator2 ,
+ class T
+ >
+ void assign_sum_increment(
+ OutputIterator first1 ,
+ OutputIterator last1 ,
+ InputIterator1 first2 ,
+ InOutIterator first3 ,
+ InputIterator2 first4 ,
+ T alpha
+ )
{
- while( x_start != x_end )
- (*out_start++) += alpha2 * ( (*x_start++) + (*y_start++) + alpha1*(*z_start++) );
+ while( first1 < last1 )
+ {
+ (*first1++) = (*first2++) + alpha * (*first3);
+ (*first3++) += (*first4++);
+ }
}
+
+
+
+
+
+
}
}
}
Modified: sandbox/odeint/boost/numeric/odeint/euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/euler.hpp 2009-11-06 11:56:00 EST (Fri, 06 Nov 2009)
@@ -83,7 +83,7 @@
time_type t ,
time_type dt )
{
- detail::it_algebra::scale_and_add( x.begin() , x.end() , dxdt.begin() , dt );
+ detail::it_algebra::increment( x.begin() , x.end() , dxdt.begin() , dt );
}
@@ -118,10 +118,11 @@
next_step( system , m_xtemp , dxdt , t , dt2 );
next_step( system , m_xtemp , t+dt2 , dt2 );
- detail::it_algebra::substract_and_assign(x.begin() ,
- x.end() ,
- m_xtemp.begin() ,
- xerr.begin() );
+ detail::it_algebra::assign_diff(
+ xerr.begin() ,
+ xerr.end() ,
+ x.begin() ,
+ m_xtemp.begin() );
}
Modified: sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp 2009-11-06 11:56:00 EST (Fri, 06 Nov 2009)
@@ -90,29 +90,20 @@
m_resizer.adjust_size( x , m_xt );
time_type dh = time_type( 0.5 ) * dt;
- time_type d6 = dt / time_type( 6.0 );
time_type th = t + dh;
- iterator iter1 , iter2 ,iter3 , iter4;
- iterator x_end = x.end() , xt_end = m_xt.end();
-
system( x , m_dxdt , t );
- scale_and_add_and_assign( x.begin() , x.end() , m_dxdt.begin() , m_xt.begin() , dh );
+ assign_sum( m_xt.begin() , m_xt.end() , x.begin() , m_dxdt.begin() , dh );
system( m_xt , m_dxt , th );
- scale_and_add_and_assign( x.begin() , x.end() , m_dxt.begin() , m_xt.begin() , dh );
+ assign_sum( m_xt.begin() , m_xt.end() , x.begin() , m_dxt.begin() , dh );
system( m_xt , m_dxm , th );
- iter1 = m_xt.begin() ; iter2 = x.begin() ; iter3 = m_dxm.begin() ; iter4 = m_dxt.begin();
- while( iter1 != xt_end )
- {
- (*iter1++) = (*iter2++) + dt * (*iter3);
- (*iter3++) += (*iter4++);
- }
+ assign_sum_increment( m_xt.begin() , m_xt.end() , x.begin() , m_dxm.begin() , m_dxt.begin() , dt );
system( m_xt , m_dxt , value_type( t + dt ) );
- scale_and_add_and_add_and_assign( m_dxdt.begin() , m_dxdt.end() , m_dxt.begin() , m_dxm.begin() , x.begin() , val2 , d6 );
-
+ increment_sum_sum( x.begin() , x.end() , m_dxdt.begin() , m_dxt.begin() , m_dxm.begin() , dt / time_type( 6.0 ) , val2
+ );
}
Added: sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp 2009-11-06 11:56:00 EST (Fri, 06 Nov 2009)
@@ -0,0 +1,136 @@
+/* Boost odeint/stepper_half_stepr.hpp header file
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ This file includes a stepper which calculates the
+ error during one step from performing two steps with
+ the halt stepsize. It works with arbitray steppers
+
+ 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)
+*/
+
+#ifndef BOOST_NUMERIC_ODEINT_EULER_HPP
+#define BOOST_NUMERIC_ODEINT_EULER_HPP
+
+#include <boost/concept_check.hpp>
+
+#include <boost/numeric/odeint/detail/iterator_algebra.hpp>
+#include <boost/numeric/odeint/concepts/state_concept.hpp>
+#include <boost/numeric/odeint/resizer.hpp>
+
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+ template<
+ class Stepper
+ >
+ class ode_step_half_step
+ {
+ // provide basic typedefs
+ public:
+
+ typedef Stepper stepper_type;
+ typedef typename Stepper::container_type container_type;
+ typedef typename Stepper::resizer_type resizer_type;
+ typedef typename Stepper::time_type time_type;
+ typedef typename container_type::value_type value_type;
+ typedef typename container_type::iterator iterator;
+
+
+
+
+
+
+
+ // private members
+ private:
+
+ container_type m_dxdt;
+ container_type m_xtemp;
+ resizer_type m_resizer;
+ stepper_type m_stepper;
+
+
+ const unsigned int order() { return m_stepper.order(); }
+
+
+
+ template< class DynamicalSystem >
+ void next_step( DynamicalSystem system ,
+ container_type &x ,
+ const container_type &dxdt ,
+ time_type t ,
+ time_type dt )
+ {
+ m_stepper.next_step( system , x , dxdt , t , dt );
+ }
+
+
+
+ template< class DynamicalSystem >
+ void next_step( DynamicalSystem system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt )
+ {
+ m_stepper.next_step( system , x , t , dt );
+ }
+
+ /*
+
+ template< class DynamicalSystem >
+ void next_step( DynamicalSystem system ,
+ container_type &x ,
+ const container_type &dxdt ,
+ time_type t ,
+ time_type dt ,
+ container_type &xerr )
+ {
+ m_resizer.adjust_size( x , xerr );
+
+ m_xtemp = x;
+ time_type dt2 = 0.5 * dt;
+
+ next_step( system , x , dxdt , t , dt );
+ next_step( system , m_xtemp , dxdt , t , dt2 );
+ next_step( system , m_xtemp , t+dt2 , dt2 );
+
+ detail::it_algebra::assign_diff(
+ xerr.begin() ,
+ xerr.end() ,
+ x.begin() ,
+ m_xtemp.begin() );
+ }
+
+
+
+ template< class DynamicalSystem >
+ void next_step( DynamicalSystem system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt ,
+ container_type &xerr )
+ {
+ m_resizer.check_size_and_resize( x , m_dxdt );
+ system( x , m_dxdt , t );
+ next_step( system , x , m_dxdt , t , dt , xerr );
+ }
+ */
+ }
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_EULER_HPP
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp 2009-11-06 11:56:00 EST (Fri, 06 Nov 2009)
@@ -37,7 +37,7 @@
const double eps_abs = 1E-3;
const double eps_rel = 1E-3;
-const size_t time_points = 101;
+const size_t time_points = 201;
typedef array<double, 3> state_type;
Modified: sandbox/odeint/libs/numeric/odeint/stuff/iterator_algebra/increment_decrement.cc
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/stuff/iterator_algebra/increment_decrement.cc (original)
+++ sandbox/odeint/libs/numeric/odeint/stuff/iterator_algebra/increment_decrement.cc 2009-11-06 11:56:00 EST (Fri, 06 Nov 2009)
@@ -12,43 +12,61 @@
using namespace std;
using namespace boost::lambda;
-template< class Iterator1 , class Iterator2 >
-void increment( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 )
+/*
+ Computes y += x
+ */
+template< class InOutIterator , class InputIterator >
+void increment( InOutIterator first1 , InOutIterator last1 , InputIterator first2 )
{
while( first1 < last1 )
(*first1++) += (*first2++);
}
-template< class Iterator1 , class Iterator2 , class T >
-void increment( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , T val )
+/*
+ Computes y += a*x
+ */
+template< class InOutIterator , class InputIterator , class T >
+void increment( InOutIterator first1 , InOutIterator last1 , InputIterator first2 , T val )
{
while( first1 < last1 )
(*first1++) += val * (*first2++);
}
-template< class Iterator1 , class Iterator2 , class Operation >
-void increment_op( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Operation op )
+/*
+ Computes y += op(x)
+ */
+template< class InOutIterator , class InputIterator , class Operation >
+void increment_op( InOutIterator first1 , InOutIterator last1 , InputIterator first2 , Operation op )
{
while( first1 < last1 )
(*first1++) += op(*first2++);
}
-template< class Iterator1 , class Iterator2 >
-void decrement( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 )
+/*
+ Computes y -= x
+ */
+template< class InOutIterator , class InputIterator >
+void decrement( InOutIterator first1 , InOutIterator last1 , InputIterator first2 )
{
while( first1 < last1 )
(*first1++) -= (*first2++);
}
-template< class Iterator1 , class Iterator2 , class T >
-void decrement( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , T val )
+/*
+ Computes y -= a*x
+ */
+template< class InOutIterator , class InputIterator , class T >
+void decrement( InOutIterator first1 , InOutIterator last1 , InputIterator first2 , T val )
{
while( first1 < last1 )
(*first1++) -= val * (*first2++);
}
-template< class Iterator1 , class Iterator2 , class Operation >
-void decrement_op( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Operation op )
+/*
+ Computes y -= f(x)
+ */
+template< class InOutIterator , class InputIterator , class Operation >
+void decrement_op( InOutIterator first1 , InOutIterator last1 , InputIterator first2 , Operation op )
{
while( first1 < last1 )
(*first1++) -= op(*first2++);
@@ -57,6 +75,54 @@
+
+/*
+ Coputes y += x1 + x2
+*/
+template< class InOutIterator , class InputIterator1 , class InputIterator2 >
+void increment_add( InOutIterator first1 , InOutIterator last1 , InputIterator1 first2 , InputIterator2 first3 )
+{
+ while( first1 < last1 )
+ (*first1++) += (*first2++) + (*first3++);
+}
+
+/*
+ Computes y += a * ( x1 + x2 )
+*/
+template< class InOutIterator , class InputIterator1 , class InputIterator2 , class T >
+void increment_add( InOutIterator first1 , InOutIterator last1 , InputIterator1 first2 , InputIterator2 first3 , T val )
+{
+ while( first1 < last1 )
+ (*first1++) += val * ( (*first2++) + (*first3++) );
+}
+
+/*
+ Computes y += a * x1 + b * x2
+*/
+template< class InOutIterator , class InputIterator1 , class InputIterator2 , class T >
+void increment_add( InOutIterator first1 , InOutIterator last1 , InputIterator1 first2 , InputIterator2 first3 , T val1 , T val2 )
+{
+ while( first1 < last1 )
+ (*first1++) += val1 * (*first2++) + val2 * (*first3++);
+}
+
+/*
+ Computes y += op( x1 , x2 )
+*/
+template< class InOutIterator , class InputIterator1 , class InputIterator2 , class Operation >
+void increment_add_op( InOutIterator first1 , InOutIterator last1 , InputIterator1 first2 , InputIterator2 first3 , Operation op )
+{
+ while( first1 < last1 )
+ (*first1++) += op( *first2++ , *first3++ );
+}
+
+
+
+
+
+
+
+
const size_t n = 100000;
const size_t cycles = 10000;
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