Squashed 'boost/' content from commit b4feb19f2

git-subtree-dir: boost
git-subtree-split: b4feb19f287ee92d87a9624b5d36b7cf46aeadeb

libs/numeric/odeint/examples/2d_lattice/Jamfile.v2

@@ -0,0 +1,13 @@
+# Copyright 2011 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# 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)
+
+project
+    : requirements 
+      <include>../../../../..
+      <define>BOOST_ALL_NO_LIB=1
+    ;
+    
+exe spreading : spreading.cpp ;

libs/numeric/odeint/examples/2d_lattice/lattice2d.hpp

@@ -0,0 +1,165 @@
+/*
+ Copyright 2011 Mario Mulansky
+ Copyright 2012-2013 Karsten Ahnert
+
+ 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)
+ */
+
+
+/* strongly nonlinear hamiltonian lattice in 2d */
+
+#ifndef LATTICE2D_HPP
+#define LATTICE2D_HPP
+
+#include <vector>
+
+#include <boost/math/special_functions/pow.hpp>
+
+using boost::math::pow;
+
+template< int Kappa , int Lambda >
+struct lattice2d {
+
+    const double m_beta;
+    std::vector< std::vector< double > > m_omega;
+
+    lattice2d( const double beta )
+        : m_beta( beta )
+    { }
+
+    template< class StateIn , class StateOut >
+    void operator()( const StateIn &q , StateOut &dpdt )
+    {
+        // q and dpdt are 2d
+        const int N = q.size();
+
+        int i;
+        for( i = 0 ; i < N ; ++i )
+        {
+            const int i_l = (i-1+N) % N;
+            const int i_r = (i+1) % N;
+            for( int j = 0 ; j < N ; ++j )
+            {
+            const int j_l = (j-1+N) % N;
+            const int j_r = (j+1) % N;
+            dpdt[i][j] = - m_omega[i][j] * pow<Kappa-1>( q[i][j] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i][j_l] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i][j_r] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i_l][j] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i_r][j] );
+            }
+        }
+    }
+
+    template< class StateIn >
+    double energy( const StateIn &q , const StateIn &p )
+    {
+        // q and dpdt are 2d
+        const int N = q.size();
+        double energy = 0.0;
+        int i;
+        for( i = 0 ; i < N ; ++i )
+        {
+            const int i_l = (i-1+N) % N;
+            const int i_r = (i+1) % N;
+            for( int j = 0 ; j < N ; ++j )
+            {
+            const int j_l = (j-1+N) % N;
+            const int j_r = (j+1) % N;
+            energy += p[i][j]*p[i][j] / 2.0
+                        + m_omega[i][j] * pow<Kappa>( q[i][j] ) / Kappa
+                + m_beta * pow<Lambda>( q[i][j] - q[i][j_l] ) / Lambda / 2
+                + m_beta * pow<Lambda>( q[i][j] - q[i][j_r] ) / Lambda / 2
+                + m_beta * pow<Lambda>( q[i][j] - q[i_l][j] ) / Lambda / 2
+                + m_beta * pow<Lambda>( q[i][j] - q[i_r][j] ) / Lambda / 2;
+            }
+        }
+        return energy;
+    }
+
+
+    template< class StateIn , class StateOut >
+    double local_energy( const StateIn &q , const StateIn &p , StateOut &energy )
+    {
+        // q and dpdt are 2d
+        const int N = q.size();
+        double e = 0.0;
+        int i;
+        for( i = 0 ; i < N ; ++i )
+        {
+            const int i_l = (i-1+N) % N;
+            const int i_r = (i+1) % N;
+            for( int j = 0 ; j < N ; ++j )
+            {
+                const int j_l = (j-1+N) % N;
+                const int j_r = (j+1) % N;
+                energy[i][j] = p[i][j]*p[i][j] / 2.0
+                    + m_omega[i][j] * pow<Kappa>( q[i][j] ) / Kappa
+                    + m_beta * pow<Lambda>( q[i][j] - q[i][j_l] ) / Lambda / 2
+                    + m_beta * pow<Lambda>( q[i][j] - q[i][j_r] ) / Lambda / 2
+                    + m_beta * pow<Lambda>( q[i][j] - q[i_l][j] ) / Lambda / 2
+                    + m_beta * pow<Lambda>( q[i][j] - q[i_r][j] ) / Lambda / 2;
+                e += energy[i][j];
+            }
+        }
+        //rescale
+        e = 1.0/e;
+        for( i = 0 ; i < N ; ++i )
+            for( int j = 0 ; j < N ; ++j )
+                energy[i][j] *= e;
+        return 1.0/e;
+    }
+
+    void load_pot( const char* filename , const double W , const double gap , 
+                   const size_t dim )
+    {
+        std::ifstream in( filename , std::ios::in | std::ios::binary );
+        if( !in.is_open() ) {
+            std::cerr << "pot file not found: " << filename << std::endl;
+            exit(0);
+        } else {
+            std::cout << "using pot file: " << filename << std::endl;
+        }
+
+        m_omega.resize( dim );
+        for( int i = 0 ; i < dim ; ++i )
+        {
+            m_omega[i].resize( dim );
+            for( size_t j = 0 ; j < dim ; ++j )
+            {
+                if( !in.good() )
+                {
+                    std::cerr << "I/O Error: " << filename << std::endl;
+                    exit(0);
+                }
+                double d;
+                in.read( (char*) &d , sizeof(d) );
+                if( (d < 0) || (d > 1.0) )
+                {
+                    std::cerr << "ERROR: " << d << std::endl;
+                    exit(0);
+                }
+                m_omega[i][j] = W*d + gap;
+            }
+        }
+
+    }
+
+    void generate_pot( const double W , const double gap , const size_t dim )
+    {
+        m_omega.resize( dim );
+        for( size_t i = 0 ; i < dim ; ++i )
+        {
+            m_omega[i].resize( dim );
+            for( size_t j = 0 ; j < dim ; ++j )
+            {
+                m_omega[i][j] = W*static_cast<double>(rand())/RAND_MAX + gap;
+            }
+        }
+    }
+
+};
+
+#endif

libs/numeric/odeint/examples/2d_lattice/nested_range_algebra.hpp

@@ -0,0 +1,46 @@
+/*
+ Copyright 2011 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ 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)
+ */
+
+
+/* nested range algebra */
+
+#ifndef NESTED_RANGE_ALGEBRA
+#define NESTED_RANGE_ALGEBRA
+
+namespace detail {
+
+    template< class Iterator1 , class Iterator2 , class Iterator3 , class Operation , class Algebra >
+    void for_each3( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3, Operation op , Algebra &algebra )
+{
+    for( ; first1 != last1 ; )
+        algebra.for_each3( *first1++ , *first2++ , *first3++ , op );
+}
+}
+
+
+template< class InnerAlgebra >
+struct nested_range_algebra
+{
+    
+    nested_range_algebra() 
+        : m_inner_algebra()
+    { }
+
+    template< class S1 , class S2 , class S3 , class Op >
+    void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op )
+    {
+        detail::for_each3( boost::begin( s1 ) , boost::end( s1 ) , boost::begin( s2 ) , boost::begin( s3 ) , op , m_inner_algebra );
+    }
+
+
+private:
+    InnerAlgebra m_inner_algebra;
+};
+
+#endif

libs/numeric/odeint/examples/2d_lattice/spreading.cpp

@@ -0,0 +1,122 @@
+/*
+ Copyright 2011 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ 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)
+ */
+
+
+/*
+ * Example of a 2D simulation of nonlinearly coupled oscillators.
+ * Program just prints final energy which should be close to the initial energy (1.0).
+ * No parallelization is employed here.
+ * Run time on a 2.3GHz Intel Core-i5: about 10 seconds for 100 steps.
+ * Compile simply via bjam or directly:
+ * g++ -O3 -I${BOOST_ROOT} -I../../../../.. spreading.cpp
+ */
+
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+#include <sys/time.h>
+
+#include <boost/ref.hpp>
+#include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp>
+
+// we use a vector< vector< double > > as state type,
+// for that some functionality has to be added for odeint to work
+#include "nested_range_algebra.hpp"
+#include "vector_vector_resize.hpp"
+
+// defines the rhs of our dynamical equation
+#include "lattice2d.hpp"
+/* dynamical equations (Hamiltonian structure):
+dqdt_{i,j} = p_{i,j}
+dpdt_{i,j} = - omega_{i,j}*q_{i,j} - \beta*[ (q_{i,j} - q_{i,j-1})^3
+                                            +(q_{i,j} - q_{i,j+1})^3
+                                            +(q_{i,j} - q_{i-1,j})^3
+                                            +(q_{i,j} - q_{i+1,j})^3 ]
+*/
+
+
+using namespace std;
+
+static const int MAX_N = 1024;//2048;
+
+static const size_t KAPPA = 2;
+static const size_t LAMBDA = 4;
+static const double W = 1.0;
+static const double gap = 0.0;
+static const size_t steps = 100;
+static const double dt = 0.1;
+
+double initial_e = 1.0;
+double beta = 1.0;
+int realization_index = 0;
+
+//the state type
+typedef vector< vector< double > > state_type;
+
+//the stepper, choose a symplectic one to account for hamiltonian structure
+//use nested_range_algebra for calculations on vector< vector< ... > >
+typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan< 
+    state_type , state_type , double , state_type , state_type , double , 
+    nested_range_algebra< boost::numeric::odeint::range_algebra > ,
+    boost::numeric::odeint::default_operations > stepper_type;
+
+double time_diff_in_ms( timeval &t1 , timeval &t2 )
+{ return (t2.tv_sec - t1.tv_sec)*1000.0 + (t2.tv_usec - t1.tv_usec)/1000.0 + 0.5; }
+
+
+int main( int argc, const char* argv[] ) {
+
+    srand( time(NULL) );
+
+    lattice2d< KAPPA , LAMBDA > lattice( beta );
+
+
+    lattice.generate_pot( W , gap , MAX_N );
+
+    state_type q( MAX_N , vector< double >( MAX_N , 0.0 ) );
+
+    state_type p( q );
+
+    state_type energy( q );
+
+    p[MAX_N/2][MAX_N/2] = sqrt( 0.5*initial_e );
+    p[MAX_N/2+1][MAX_N/2] = sqrt( 0.5*initial_e );
+    p[MAX_N/2][MAX_N/2+1] = sqrt( 0.5*initial_e );
+    p[MAX_N/2+1][MAX_N/2+1] = sqrt( 0.5*initial_e );
+
+    cout.precision(10);
+
+    lattice.local_energy( q , p , energy );
+    double e=0.0;
+    for( size_t i=0 ; i<energy.size() ; ++i )
+        for( size_t j=0 ; j<energy[i].size() ; ++j )
+        {
+            e += energy[i][j];
+        }
+
+    cout << "initial energy: " << lattice.energy( q , p ) << endl;
+
+    timeval elapsed_time_start , elapsed_time_end;
+    gettimeofday(&elapsed_time_start , NULL);
+
+    stepper_type stepper;
+
+    for( size_t step=0 ; step<=steps ; ++step )
+    {
+        stepper.do_step( boost::ref( lattice ) , 
+                         make_pair( boost::ref( q ) , boost::ref( p ) ) , 
+                         0.0 , 0.1 );
+    }
+
+    gettimeofday(&elapsed_time_end , NULL);
+    double elapsed_time = time_diff_in_ms( elapsed_time_start , elapsed_time_end );
+    cout << steps << " steps in " << elapsed_time/1000 << " s (energy: " << lattice.energy( q , p ) << ")" << endl;
+}

libs/numeric/odeint/examples/2d_lattice/vector_vector_resize.hpp

@@ -0,0 +1,105 @@
+/*
+ Copyright 2011 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ 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)
+ */
+
+
+/* reserved vector */
+
+#ifndef VECTOR_VECTOR_RESIZE_HPP
+#define VECTOR_VECTOR_RESIZE_HPP
+
+#include <vector>
+
+#include <boost/range.hpp>
+
+namespace boost { namespace numeric { namespace odeint {
+
+template<>
+struct is_resizeable< std::vector< std::vector< double > > >
+{
+    typedef boost::true_type type;
+    const static bool value = type::value;
+};
+
+template<>
+struct same_size_impl< std::vector< std::vector< double > > , std::vector< std::vector< double > > >
+{
+    typedef std::vector< std::vector< double > > state_type;
+
+    static bool same_size( const state_type &x1 ,
+                           const state_type &x2 )
+    {
+        bool same = ( boost::size( x1 ) == boost::size( x2 ) );
+        if( !same )
+            return false;
+        typename state_type::const_iterator begin1 = boost::begin( x1 );
+        typename state_type::const_iterator begin2 = boost::begin( x2 );
+        while( begin1 != boost::end( x1 ) )
+            same &= ( boost::size( *begin1++ ) == boost::size( *begin2++ ) );
+        return same;
+    }
+};
+
+template<>
+struct resize_impl< std::vector< std::vector< double > > , std::vector< std::vector< double > > >
+{
+    typedef std::vector< std::vector< double > > state_type;
+
+    static void resize( state_type &x1 , const state_type &x2 )
+    {
+        x1.resize( boost::size( x2 ) );
+        typename state_type::iterator begin1 = boost::begin( x1 );
+        typename state_type::const_iterator begin2 = boost::begin( x2 );
+        while( begin1 != boost::end( x1 ) )
+            (*begin1++).resize( boost::size( *begin2++ ) );
+    }
+};
+
+template<>
+struct state_wrapper< std::vector< std::vector< double > > >
+{
+    typedef std::vector< std::vector< double > > state_type;
+    typedef state_wrapper< state_type > state_wrapper_type;
+    typedef boost::true_type is_resizeable;
+
+    state_type m_v;
+
+    template< class State >
+    bool same_size( const State &x )
+    {
+        bool same = ( boost::size( m_v ) == boost::size( x ) );
+        if( !same )
+            return false;
+        typename state_type::iterator begin1 = boost::begin( m_v );
+        typename State::const_iterator begin2 = boost::begin( x );
+        while( begin1 != boost::end( m_v ) )
+            same &= ( boost::size( *begin1++ ) == boost::size( *begin2++ ) );
+        return same;
+    }
+
+    template< class State >
+    bool resize( const State &x )
+    {
+        if( !same_size( x ) )
+        {
+            m_v.resize( boost::size( x ) );
+            typename state_type::iterator begin1 = boost::begin( m_v );
+            typename State::const_iterator begin2 = boost::begin( x );
+            while( begin1 != boost::end( m_v ) )
+                (*begin1++).resize( boost::size( *begin2++ ) );
+
+            return true;
+        } else
+            return false;
+    }
+ 
+};
+
+} } }
+
+#endif

libs/numeric/odeint/examples/Jamfile.v2

@@ -0,0 +1,51 @@
+# Copyright 2009-2013 Karsten Ahnert
+# Copyright 2010-2013 Mario Mulansky
+# Copyright 2013 Pascal Germroth
+# 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)
+
+
+project
+    : requirements 
+      <define>BOOST_ALL_NO_LIB=1
+    : 
+    ;
+    
+
+exe harmonic_oscillator : harmonic_oscillator.cpp ;
+exe solar_system : solar_system.cpp ;
+exe chaotic_system : chaotic_system.cpp ;
+exe stiff_system : stiff_system.cpp ;
+exe fpu : fpu.cpp ;
+exe phase_oscillator_ensemble : phase_oscillator_ensemble.cpp ;
+exe harmonic_oscillator_units : harmonic_oscillator_units.cpp : <toolset>clang:<build>no ;
+exe stuart_landau : stuart_landau.cpp ;
+exe two_dimensional_phase_lattice : two_dimensional_phase_lattice.cpp ;
+exe bulirsch_stoer : bulirsch_stoer.cpp ;
+exe elliptic_functions : elliptic_functions.cpp ;
+exe resizing_lattice : resizing_lattice.cpp ;
+exe list_lattice : list_lattice.cpp ;
+exe stepper_details : stepper_details.cpp ;
+exe my_vector : my_vector.cpp ;
+exe lorenz : lorenz.cpp ;
+exe lorenz_point : lorenz_point.cpp ;
+exe van_der_pol_stiff : van_der_pol_stiff.cpp ;
+exe simple1d : simple1d.cpp ;
+exe stochastic_euler : stochastic_euler.cpp ;
+exe generation_functions : generation_functions.cpp ;
+exe heun : heun.cpp ;
+exe bind_member_functions : bind_member_functions.cpp ;
+exe bind_member_functions_cpp11 : bind_member_functions_cpp11.cpp : <cxxflags>-std=c++0x ;
+exe molecular_dynamics : molecular_dynamics.cpp : <cxxflags>-std=c++0x ;
+exe molecular_dynamics_cells : molecular_dynamics_cells.cpp : <cxxflags>-std=c++0x ;
+exe abm_precision : abm_precision.cpp ;
+exe integrate_times : integrate_times.cpp ;
+exe find_crossing : find_crossing.cpp : <cxxflags>-std=c++0x ;
+
+build-project multiprecision ;
+# build-project mtl ;
+# build-project ublas ;
+# build-project gmpxx ;
+# build-project openmp ;
+# build-project mpi ;

libs/numeric/odeint/examples/abm_precision.cpp

@@ -0,0 +1,84 @@
+/*
+ * abm_precision.cpp
+ *
+ * example to check the order of the multi-step methods
+ *
+ * Copyright 2009-2013 Karsten Ahnert
+ * Copyright 2009-2013 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)
+ */
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/array.hpp>
+#include <boost/numeric/odeint.hpp>
+
+using namespace boost::numeric::odeint;
+
+const int Steps = 4;
+
+typedef double value_type;
+
+typedef boost::array< double , 2 > state_type;
+
+typedef runge_kutta_fehlberg78<state_type> initializing_stepper_type;
+typedef adams_bashforth_moulton< Steps , state_type > stepper_type;
+//typedef adams_bashforth< Steps , state_type > stepper_type;
+
+// harmonic oscillator, analytic solution x[0] = sin( t )
+struct osc
+{
+    void operator()( const state_type &x , state_type &dxdt , const double t ) const
+    {
+        dxdt[0] = x[1];
+        dxdt[1] = -x[0];
+    }
+};
+
+int main()
+{
+    stepper_type stepper;
+    initializing_stepper_type init_stepper;
+    const int o = stepper.order()+1; //order of the error is order of approximation + 1
+
+    const state_type x0 = {{ 0.0 , 1.0 }};
+    state_type x1 = x0;
+    double t = 0.0;
+    double dt = 0.25;
+    // initialization, does a number of steps already to fill internal buffer, t is increased
+    // we use the rk78 as initializing stepper
+    stepper.initialize( boost::ref(init_stepper) , osc() , x1 , t , dt );
+    // do a number of steps to fill the buffer with results from adams bashforth
+    for( size_t n=0 ; n < stepper.steps ; ++n )
+    {
+        stepper.do_step( osc() , x1 , t , dt );
+        t += dt;
+    }
+    double A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
+    double phi = std::asin(x1[0]/A) - t;
+    // now we do the actual step
+    stepper.do_step( osc() , x1 , t , dt );
+    // only examine the error of the adams-bashforth-moulton step, not the initialization
+    const double f = 2.0 * std::abs( A*sin(t+dt+phi) - x1[0] ) / std::pow( dt , o ); // upper bound
+
+    std::cout << "# " << o << " , " << f << std::endl;
+
+    /* as long as we have errors above machine precision */
+    while( f*std::pow( dt , o ) > 1E-16 )
+    {
+        x1 = x0;
+        t = 0.0;
+        stepper.initialize( boost::ref(init_stepper) , osc() , x1 , t , dt );
+        A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
+        phi = std::asin(x1[0]/A) - t;
+        // now we do the actual step
+        stepper.do_step( osc() , x1 , t , dt );
+        // only examine the error of the adams-bashforth-moulton step, not the initialization
+        std::cout << dt << '\t' <<  std::abs( A*sin(t+dt+phi) - x1[0] ) << std::endl;
+        dt *= 0.5;
+    }
+}

libs/numeric/odeint/examples/adaptive_iterator.cpp

@@ -0,0 +1,361 @@
+/*
+ * adaptive_iterator.cpp
+ *
+ * Copyright 2012-2013 Karsten Ahnert
+ * Copyright 2012 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)
+ */
+
+
+
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <algorithm>
+#include <cassert>
+
+#include <boost/array.hpp>
+
+#include <boost/range/algorithm.hpp>
+#include <boost/range/adaptor/filtered.hpp>
+#include <boost/range/numeric.hpp>
+
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
+#include <boost/numeric/odeint/stepper/generation.hpp>
+
+#include <boost/numeric/odeint/iterator/adaptive_iterator.hpp>
+#include <boost/numeric/odeint/iterator/adaptive_time_iterator.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+struct lorenz
+{
+    template< class State , class Deriv >
+    void operator()( const State &x , Deriv &dxdt , double t ) const
+    {
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+
+#include <typeinfo>
+
+int main( int argc , char **argv )
+{
+    typedef boost::array< double , 3 > state_type;
+
+    /*
+     * Controlled steppers with time iterator
+     */
+
+    // std::for_each
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        std::for_each( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                       make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
+                       []( const std::pair< const state_type&, double > &x ) {
+                           std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
+    }
+
+    // std::copy_if
+    {
+        std::vector< pair< state_type , double > > res;
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        std::copy_if( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                      make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
+                      std::back_inserter( res ) ,
+                      []( const pair< const state_type& , double > &x ) {
+                          return ( x.first[0] > 0.0 ) ? true : false; } );
+        for( size_t i=0 ; i<res.size() ; ++i )
+            cout << res[i].first[0] << tab << res[i].first[1] << tab << res[i].first[2] << "\n";
+    }
+
+    // std::accumulate
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        double res = std::accumulate( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                      make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
+                                      0.0 ,
+                                      []( double sum , const pair< const state_type& , double > &x ) {
+                                          return sum + x.first[0]; } );
+        cout << res << endl;
+    }
+
+
+    // std::transform
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        vector< double > weights;
+        std::transform( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                        make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
+                        back_inserter( weights ) ,
+                        []( const pair< const state_type& , double > &x ) {
+                            return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
+        for( size_t i=0 ; i<weights.size() ; ++i )
+            cout << weights[i] << "\n";
+    }
+
+
+
+
+
+
+
+
+
+
+
+
+
+    /*
+     * Boost.Range versions of controlled stepper with time iterator
+     */
+
+
+    // boost::range::for_each
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::for_each( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                []( const std::pair< const state_type& , double > &x ) {
+                                    std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
+    }
+
+
+    // boost::range::copy with filtered adaptor (simulating std::copy_if)
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        std::vector< std::pair< state_type , double > > res;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::copy( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
+                            boost::adaptors::filtered( [] ( const pair< const state_type& , double > &x ) { return ( x.first[0] > 0.0 ); } ) ,
+                            std::back_inserter( res ) );
+        for( size_t i=0 ; i<res.size() ; ++i )
+            cout << res[i].first[0] << tab << res[i].first[1] << tab << res[i].first[2] << "\n";
+    }
+
+    // boost::range::accumulate
+    {
+        //[adaptive_time_iterator_accumulate_range
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        double res = boost::accumulate( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
+                                        []( double sum , const pair< const state_type& , double > &x ) {
+                                            return sum + x.first[0]; } );
+        cout << res << endl;
+        //]
+    }
+
+
+    //  boost::range::transform
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        vector< double > weights;
+        boost::transform( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
+                          []( const pair< const state_type& , double > &x ) {
+                              return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
+        for( size_t i=0 ; i<weights.size() ; ++i )
+            cout << weights[i] << "\n";
+    }
+
+
+    // boost::range::find with time iterator
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        auto iter = boost::find_if( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                    []( const std::pair< const state_type & , double > &x ) {
+                                        return ( x.first[0] < 0.0 ); } );
+        cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";
+    }
+
+
+
+
+
+
+
+
+
+
+
+
+    // /*
+    //  * Boost.Range versions for dense output steppers
+    //  */
+
+    // // boost::range::for_each
+    // {
+    //     runge_kutta_dopri5< state_type > stepper;
+    //     state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    //     boost::range::for_each( make_adaptive_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+    //                             []( const state_type &x ) {
+    //                                 std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
+    // }
+
+    
+    // // boost::range::for_each with time iterator
+    // {
+    //     runge_kutta_dopri5< state_type > stepper;
+    //     state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    //     boost::range::for_each( make_adaptive_time_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+    //                             []( const std::pair< state_type& , double > &x ) {
+    //                                 std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
+
+    // }
+
+
+
+
+
+    /*
+     * Pure iterators for controlled stepper without time iterator
+     */
+
+    // std::for_each
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        std::for_each( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                       make_adaptive_iterator_end( stepper , lorenz() , x ) ,
+                       []( const state_type& x ) {
+                           std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
+    }
+
+    // std::copy_if
+    {
+        std::vector< state_type > res;
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        std::copy_if( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                      make_adaptive_iterator_end( stepper , lorenz() , x ) ,
+                      std::back_inserter( res ) ,
+                      []( const state_type& x ) {
+                          return ( x[0] > 0.0 ) ? true : false; } );
+        for( size_t i=0 ; i<res.size() ; ++i )
+            cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
+    }
+
+    // std::accumulate
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        double res = std::accumulate( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                      make_adaptive_iterator_end( stepper , lorenz() , x ) ,
+                                      0.0 ,
+                                      []( double sum , const state_type& x ) {
+                                          return sum + x[0]; } );
+        cout << res << endl;
+    }
+
+
+    // std::transform
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        vector< double > weights;
+        std::transform( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                        make_adaptive_iterator_end( stepper , lorenz() , x ) ,
+                        back_inserter( weights ) ,
+                        []( const state_type& x ) {
+                            return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
+        for( size_t i=0 ; i<weights.size() ; ++i )
+            cout << weights[i] << "\n";
+    }
+
+
+
+
+
+
+
+
+
+
+    /*
+     * Boost.Range versions of controlled stepper WITHOUT time iterator
+     */
+
+
+    // boost::range::for_each
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::for_each( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                []( const state_type &x ) {
+                                    std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
+    }
+
+
+    // boost::range::copy with filtered adaptor (simulating std::copy_if)
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        std::vector< state_type > res;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::copy( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
+                            boost::adaptors::filtered( [] ( const state_type& x ) { return ( x[0] > 0.0 ); } ) ,
+                            std::back_inserter( res ) );
+        for( size_t i=0 ; i<res.size() ; ++i )
+            cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
+    }
+
+    // boost::range::accumulate
+    {
+        //[adaptive_iterator_accumulate_range
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        double res = boost::accumulate( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
+                                        []( double sum , const state_type& x ) {
+                                            return sum + x[0]; } );
+        cout << res << endl;
+        //]
+    }
+
+
+    //  boost::range::transform
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        vector< double > weights;
+        boost::transform( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
+                          []( const state_type& x ) {
+                              return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
+        for( size_t i=0 ; i<weights.size() ; ++i )
+            cout << weights[i] << "\n";
+    }
+
+
+    // boost::range::find 
+    {
+        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        auto iter = boost::find_if( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                    []( const state_type &x ) {
+                                        return ( x[0] < 0.0 ); } );
+        cout << (*iter)[0] << "\t" << (*iter)[1] << "\t" << (*iter)[2] << "\n";
+    }
+
+
+
+
+
+    return 0;
+}

libs/numeric/odeint/examples/bind_member_functions.cpp

@@ -0,0 +1,126 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/examples/bind_member_functions.hpp
+
+ [begin_description]
+ tba.
+ [end_description]
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 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)
+ */
+
+#include <iostream>
+
+#include <boost/numeric/odeint.hpp>
+
+namespace odeint = boost::numeric::odeint;
+
+typedef boost::array< double , 3 > state_type;
+
+//[ ode_wrapper
+template< class Obj , class Mem >
+class ode_wrapper
+{
+    Obj m_obj;
+    Mem m_mem;
+
+public:
+
+    ode_wrapper( Obj obj , Mem mem ) : m_obj( obj ) , m_mem( mem ) { }
+
+    template< class State , class Deriv , class Time >
+    void operator()( const State &x , Deriv &dxdt , Time t )
+    {
+        (m_obj.*m_mem)( x , dxdt , t );
+    }
+};
+
+template< class Obj , class Mem >
+ode_wrapper< Obj , Mem > make_ode_wrapper( Obj obj , Mem mem )
+{
+    return ode_wrapper< Obj , Mem >( obj , mem );
+}
+//]
+
+
+template< class Obj , class Mem >
+class observer_wrapper
+{
+    Obj m_obj;
+    Mem m_mem;
+
+public:
+
+    observer_wrapper( Obj obj , Mem mem ) : m_obj( obj ) , m_mem( mem ) { }
+
+    template< class State , class Time >
+    void operator()( const State &x , Time t )
+    {
+        (m_obj.*m_mem)( x , t );
+    }
+};
+
+template< class Obj , class Mem >
+observer_wrapper< Obj , Mem > make_observer_wrapper( Obj obj , Mem mem )
+{
+    return observer_wrapper< Obj , Mem >( obj , mem );
+}
+
+
+
+//[ bind_member_function
+struct lorenz
+{
+    void ode( const state_type &x , state_type &dxdt , double t ) const
+    {
+        dxdt[0] = 10.0 * ( x[1] - x[0] );
+        dxdt[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -8.0 / 3.0 * x[2] + x[0] * x[1];
+    }
+};
+
+int main( int argc , char *argv[] )
+{
+    using namespace boost::numeric::odeint;
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    integrate_const( runge_kutta4< state_type >() , make_ode_wrapper( lorenz() , &lorenz::ode ) ,
+                     x , 0.0 , 10.0 , 0.01 );
+    return 0;
+}
+//]
+
+
+/*
+struct lorenz
+{
+    void ode( const state_type &x , state_type &dxdt , double t ) const
+    {
+        dxdt[0] = 10.0 * ( x[1] - x[0] );
+        dxdt[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -8.0 / 3.0 * x[2] + x[0] * x[1];
+    }
+
+    void obs( const state_type &x , double t ) const
+    {
+        std::cout << t << " " << x[0] << " " << x[1] << " " << x[2] << "\n";
+    }
+};
+
+int main( int argc , char *argv[] )
+{
+    using namespace boost::numeric::odeint;
+
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    integrate_const( runge_kutta4< state_type >() ,
+                     make_ode_wrapper( lorenz() , &lorenz::ode ) ,
+                     x , 0.0 , 10.0 , 0.01 ,
+                     make_observer_wrapper( lorenz() , &lorenz::obs ) );
+
+    return 0;
+}
+*/

libs/numeric/odeint/examples/bind_member_functions_cpp11.cpp

@@ -0,0 +1,56 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/examples/bind_member_functions.hpp
+
+ [begin_description]
+ tba.
+ [end_description]
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 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)
+ */
+
+#include <iostream>
+#include <array>
+#include <type_traits>
+
+#include <boost/numeric/odeint.hpp>
+
+namespace odeint = boost::numeric::odeint;
+
+
+
+typedef std::array< double , 3 > state_type;
+
+struct lorenz
+{
+    void ode( const state_type &x , state_type &dxdt , double t ) const
+    {
+        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] = -b * x[2] + x[0] * x[1];
+    }
+};
+
+int main( int argc , char *argv[] )
+{
+    using namespace boost::numeric::odeint;
+    //[ bind_member_function_cpp11
+    namespace pl = std::placeholders;
+
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    integrate_const( runge_kutta4< state_type >() ,
+                     std::bind( &lorenz::ode , lorenz() , pl::_1 , pl::_2 , pl::_3 ) ,
+                     x , 0.0 , 10.0 , 0.01  );
+    //]
+    return 0;
+}
+

libs/numeric/odeint/examples/bulirsch_stoer.cpp

@@ -0,0 +1,100 @@
+/*
+ * bulirsch_stoer.cpp
+ *
+ * Copyright 2011-2013 Mario Mulansky
+ * Copyright 2011-2012 Karsten Ahnert
+ *
+ * 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)
+ */
+
+#include <iostream>
+#include <fstream>
+#define _USE_MATH_DEFINES
+#include <cmath>
+
+#include <boost/array.hpp>
+#include <boost/ref.hpp>
+
+#include <boost/numeric/odeint/config.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
+#include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::array< double , 1 > state_type;
+
+/*
+ * x' = ( - x*sin t  + 2 tan x ) y
+ * with x( pi/6 ) = 2/sqrt(3) the analytic solution is 1/cos t
+ */
+
+void rhs( const state_type &x , state_type &dxdt , const double t )
+{
+    dxdt[0] = ( - x[0] * sin( t ) + 2.0 * tan( t ) ) * x[0];
+}
+
+void rhs2( const state_type &x , state_type &dxdt , const double t )
+{
+    dxdt[0] = sin(t);
+}
+
+
+ofstream out;
+
+void write_out( const state_type &x , const double t )
+{
+    out << t << '\t' << x[0] << endl;
+}
+
+int main()
+{
+    bulirsch_stoer_dense_out< state_type > stepper( 1E-8 , 0.0 , 0.0 , 0.0 );
+    bulirsch_stoer< state_type > stepper2( 1E-8 , 0.0 , 0.0 , 0.0 );
+
+    state_type x = {{ 2.0 / sqrt(3.0) }};
+
+    double t = M_PI/6.0;
+    //double t = 0.0;
+    double dt = 0.01;
+    double t_end = M_PI/2.0 - 0.1;
+    //double t_end = 100.0;
+
+    out.open( "bs.dat" );
+    out.precision(16);
+    integrate_const( stepper , rhs , x , t , t_end , dt , write_out );
+    out.close();
+
+    x[0] = 2.0 / sqrt(3.0);
+
+    out.open( "bs2.dat" );
+    out.precision(16);
+    integrate_adaptive( stepper , rhs , x , t , t_end , dt , write_out );
+    out.close();
+
+    x[0] = 2.0 / sqrt(3.0);
+
+    out.open( "bs3.dat" );
+    out.precision(16);
+    integrate_adaptive( stepper2 , rhs , x , t , t_end , dt , write_out );
+    out.close();
+
+
+    typedef runge_kutta_dopri5< state_type > dopri5_type;
+    typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
+    typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
+
+    dense_output_dopri5_type dopri5 = make_dense_output( 1E-9 , 1E-9 , dopri5_type() );
+
+    x[0] = 2.0 / sqrt(3.0);
+
+    out.open( "bs4.dat" );
+    out.precision(16);
+    integrate_adaptive( dopri5 , rhs , x , t , t_end , dt , write_out );
+    out.close();
+
+}

libs/numeric/odeint/examples/chaotic_system.cpp

@@ -0,0 +1,119 @@
+/*
+ * chaotic_system.cpp
+ *
+ * This example demonstrates how one can use odeint to determine the Lyapunov
+ * exponents of a chaotic system namely the well known Lorenz system. Furthermore,
+ * it shows how odeint interacts with boost.range.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2011-2013 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)
+ */
+
+
+#include <iostream>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+#include "gram_schmidt.hpp"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+//[ system_function_without_perturbations
+struct lorenz
+{
+    template< class State , class Deriv >
+    void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
+    {
+        typename boost::range_iterator< const State >::type x = boost::begin( x_ );
+        typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
+
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+//]
+
+
+
+//[ system_function_with_perturbations
+const size_t n = 3;
+const size_t num_of_lyap = 3;
+const size_t N = n + n*num_of_lyap;
+
+typedef boost::array< double , N > state_type;
+typedef boost::array< double , num_of_lyap > lyap_type;
+
+void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
+{
+    lorenz()( x , dxdt , t );
+
+    for( size_t l=0 ; l<num_of_lyap ; ++l )
+    {
+        const double *pert = x.begin() + 3 + l * 3;
+        double *dpert = dxdt.begin() + 3 + l * 3;
+        dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
+        dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
+        dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
+    }
+}
+//]
+
+
+
+
+
+int main( int argc , char **argv )
+{
+    state_type x;
+    lyap_type lyap;
+
+    fill( x.begin() , x.end() , 0.0 );
+    x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
+
+    const double dt = 0.01;
+
+    //[ integrate_transients_with_range
+    // explicitly choose range_algebra to override default choice of array_algebra
+    runge_kutta4< state_type , double , state_type , double , range_algebra > rk4;
+
+    // perform 10000 transient steps
+    integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );
+    //]
+
+    //[ lyapunov_full_code
+    fill( x.begin()+n , x.end() , 0.0 );
+    for( size_t i=0 ; i<num_of_lyap ; ++i ) x[n+n*i+i] = 1.0;
+    fill( lyap.begin() , lyap.end() , 0.0 );
+
+    double t = 0.0;
+    size_t count = 0;
+    while( true )
+    {
+
+        t = integrate_n_steps( rk4 , lorenz_with_lyap , x , t , dt , 100 );
+        gram_schmidt< num_of_lyap >( x , lyap , n );
+        ++count;
+
+        if( !(count % 100000) )
+        {
+            cout << t;
+            for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
+            cout << endl;
+        }
+    }
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/const_step_iterator.cpp

@@ -0,0 +1,296 @@
+/*
+ * const_step_iterator.cpp
+ *
+ * Copyright 2012-2013 Karsten Ahnert
+ * Copyright 2013 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)
+ *
+ * several examples for using iterators
+ */
+
+
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <algorithm>
+#include <array>
+#include <cassert>
+
+#include <boost/range/algorithm.hpp>
+#include <boost/range/adaptor/filtered.hpp>
+#include <boost/range/numeric.hpp>
+
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/generation.hpp>
+#include <boost/numeric/odeint/iterator/const_step_iterator.hpp>
+#include <boost/numeric/odeint/iterator/const_step_time_iterator.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+struct lorenz
+{
+    template< class State , class Deriv >
+    void operator()( const State &x , Deriv &dxdt , double t ) const
+    {
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+
+
+
+int main( int argc , char **argv )
+{
+    typedef std::array< double , 3 > state_type;
+
+    // std::for_each
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        std::for_each( make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                       make_const_step_time_iterator_end( stepper , lorenz() , x ) ,
+                       []( const std::pair< const state_type&, double > &x ) {
+                           std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
+    }
+
+    // std::copy_if
+    {
+        std::vector< state_type > res;
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        std::copy_if( make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                      make_const_step_iterator_end( stepper , lorenz() , x ) ,
+                      std::back_inserter( res ) ,
+                      []( const state_type& x ) {
+                          return ( x[0] > 0.0 ) ? true : false; } );
+        for( size_t i=0 ; i<res.size() ; ++i )
+            cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
+    }
+
+    // std::accumulate
+    {
+        //[ const_step_iterator_accumulate
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        double res = std::accumulate( make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                      make_const_step_iterator_end( stepper , lorenz() , x ) ,
+                                      0.0 ,
+                                      []( double sum , const state_type &x ) {
+                                          return sum + x[0]; } );
+        cout << res << endl;
+        //]
+    }
+
+
+    // std::transform
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        vector< double > weights;
+        std::transform( make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                        make_const_step_iterator_end( stepper , lorenz() , x ) ,
+                        back_inserter( weights ) ,
+                        []( const state_type &x ) {
+                            return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
+        for( size_t i=0 ; i<weights.size() ; ++i )
+            cout << weights[i] << "\n";
+    }
+
+
+
+    // std::transform with time_iterator
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        vector< double > weights;
+        std::transform( make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                        make_const_step_time_iterator_end( stepper , lorenz() , x ) ,
+                        back_inserter( weights ) ,
+                        []( const std::pair< const state_type &, double > &x ) {
+                            return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
+        for( size_t i=0 ; i<weights.size() ; ++i )
+            cout << weights[i] << "\n";
+    }
+
+
+
+
+
+
+
+
+
+
+    // /*
+    //  * Boost.Range versions
+    //  */
+
+
+    // boost::range::for_each
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::for_each( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                []( const state_type &x ) {
+                                    std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
+    }
+
+    // boost::range::for_each with time iterator
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::for_each( make_const_step_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                []( const std::pair< const state_type& , double > &x ) {
+                                    std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
+    }
+
+
+    // boost::range::copy with filtered adaptor (simulating std::copy_if)
+    {
+        std::vector< state_type > res;
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::copy( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
+                            boost::adaptors::filtered( [] ( const state_type &x ) { return ( x[0] > 0.0 ); } ) ,
+                            std::back_inserter( res ) );
+        for( size_t i=0 ; i<res.size() ; ++i )
+            cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
+    }
+
+    // boost::range::accumulate
+    {
+        //[const_step_iterator_accumulate_range
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        double res = boost::accumulate( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
+                                        []( double sum , const state_type &x ) {
+                                            return sum + x[0]; } );
+        cout << res << endl;
+        //]
+    }
+
+    // boost::range::accumulate with time iterator
+    {
+        //[const_step_time_iterator_accumulate_range
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        double res = boost::accumulate( make_const_step_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
+                                        []( double sum , const std::pair< const state_type &, double > &x ) {
+                                            return sum + x.first[0]; } );
+        cout << res << endl;
+        //]
+    }
+
+
+    //  boost::range::transform
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        vector< double > weights;
+        boost::transform( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
+                          []( const state_type &x ) {
+                              return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
+        for( size_t i=0 ; i<weights.size() ; ++i )
+            cout << weights[i] << "\n";
+    }
+
+
+    // boost::range::find with time iterator
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        auto iter = boost::find_if( make_const_step_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                    []( const std::pair< const state_type & , double > &x ) {
+                                        return ( x.first[0] < 0.0 ); } );
+        cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";
+                                    
+    }
+
+
+
+
+
+
+
+
+
+
+
+
+    /*
+     * Boost.Range versions for dense output steppers
+     */
+
+    // boost::range::for_each
+    {
+        runge_kutta_dopri5< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::for_each( make_const_step_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                []( const state_type &x ) {
+                                    std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
+    }
+
+    
+    // boost::range::for_each with time iterator
+    {
+        runge_kutta_dopri5< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        boost::range::for_each( make_const_step_time_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
+                                []( const std::pair< const state_type& , double > &x ) {
+                                    std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
+
+    }
+
+
+
+
+
+    /*
+     * Pure iterators
+     */
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        auto first = make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 );
+        auto last  = make_const_step_iterator_end( stepper , lorenz() , x );
+        while( first != last )
+        {
+            assert( last != first );
+            cout << (*first)[0] << tab << (*first)[1] << tab << (*first)[2] << "\n";
+            ++first;
+        }
+    }
+
+    {
+        runge_kutta4< state_type > stepper;
+        state_type x = {{ 10.0 , 10.0 , 10.0 }};
+        auto first = make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 );
+        auto last  = make_const_step_time_iterator_end( stepper , lorenz() , x );
+        while( first != last )
+        {
+            assert( last != first );
+            cout << first->second << tab << first->first[0] << tab << first->first[1] << tab << first->first[2] << "\n";
+            ++first;
+        }
+    }
+
+
+
+
+
+
+
+    return 0;
+}

libs/numeric/odeint/examples/elliptic.py

@@ -0,0 +1,31 @@
+"""
+ Copyright 2011 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ Stochastic euler stepper example and Ornstein-Uhlenbeck process
+
+ 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)
+"""
+
+
+from pylab import *
+from scipy import special
+
+data1 = loadtxt("elliptic1.dat")
+data2 = loadtxt("elliptic2.dat")
+data3 = loadtxt("elliptic3.dat")
+
+sn1,cn1,dn1,phi1 = special.ellipj( data1[:,0] , 0.51 )
+sn2,cn2,dn2,phi2 = special.ellipj( data2[:,0] , 0.51 )
+sn3,cn3,dn3,phi3 = special.ellipj( data3[:,0] , 0.51 )
+
+semilogy( data1[:,0] , abs(data1[:,1]-sn1) )
+semilogy( data2[:,0] , abs(data2[:,1]-sn2) , 'ro' )
+semilogy( data3[:,0] , abs(data3[:,1]-sn3) , '--' )
+
+show()
+
+
+

libs/numeric/odeint/examples/elliptic_functions.cpp

@@ -0,0 +1,89 @@
+/*
+ * elliptic_functions.cpp
+ *
+ * Copyright 2011-2013 Mario Mulansky
+ * Copyright 2011-2012 Karsten Ahnert
+ *
+ * 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)
+ */
+
+
+
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/config.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
+#include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::array< double , 3 > state_type;
+
+/*
+ * x1' = x2*x3
+ * x2' = -x1*x3
+ * x3' = -m*x1*x2
+ */
+
+void rhs( const state_type &x , state_type &dxdt , const double t )
+{
+    static const double m = 0.51;
+
+    dxdt[0] = x[1]*x[2];
+    dxdt[1] = -x[0]*x[2];
+    dxdt[2] = -m*x[0]*x[1];
+}
+
+ofstream out;
+
+void write_out( const state_type &x , const double t )
+{
+    out << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl;
+}
+
+int main()
+{
+    bulirsch_stoer_dense_out< state_type > stepper( 1E-9 , 1E-9 , 1.0 , 0.0 );
+
+    state_type x1 = {{ 0.0 , 1.0 , 1.0 }};
+
+    double t = 0.0;
+    double dt = 0.01;
+
+    out.open( "elliptic1.dat" );
+    out.precision(16);
+    integrate_const( stepper , rhs , x1 , t , 100.0 , dt , write_out );
+    out.close();
+
+    state_type x2 = {{ 0.0 , 1.0 , 1.0 }};
+
+    out.open( "elliptic2.dat" );
+    out.precision(16);
+    integrate_adaptive( stepper , rhs , x2 , t , 100.0 , dt , write_out );
+    out.close();
+
+    typedef runge_kutta_dopri5< state_type > dopri5_type;
+    typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
+    typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
+
+
+    dense_output_dopri5_type dopri5 = make_dense_output( 1E-9 , 1E-9 , dopri5_type() );
+    //dense_output_dopri5_type dopri5( controlled_dopri5_type( default_error_checker< double >( 1E-9 , 0.0 , 0.0 , 0.0 )  ) );
+
+    state_type x3 = {{ 0.0 , 1.0 , 1.0 }};
+
+    out.open( "elliptic3.dat" );
+    out.precision(16);
+    integrate_adaptive( dopri5 , rhs , x3 , t , 100.0 , dt , write_out );
+    out.close();
+}

libs/numeric/odeint/examples/find_crossing.cpp

@@ -0,0 +1,134 @@
+/*
+ * find_crossing.cpp
+ *
+ * Finds the energy threshold crossing for a damped oscillator.
+ * The algorithm uses a dense out stepper with find_if to first find an
+ * interval containing the threshold crossing and the utilizes the dense out
+ * functionality with a bisection to further refine the interval until some
+ * desired precision is reached.
+ *
+ * Copyright 2015 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)
+ */
+
+
+
+#include <iostream>
+#include <utility>
+#include <algorithm>
+#include <array>
+
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/generation.hpp>
+#include <boost/numeric/odeint/iterator/adaptive_iterator.hpp>
+
+namespace odeint = boost::numeric::odeint;
+
+typedef std::array<double, 2> state_type;
+
+const double gam = 1.0;  // damping strength
+
+void damped_osc(const state_type &x, state_type &dxdt, const double /*t*/)
+{
+    dxdt[0] = x[1];
+    dxdt[1] = -x[0] - gam * x[1];
+}
+
+
+struct energy_condition {
+
+    // defines the threshold crossing in terms of a boolean functor
+
+    double m_min_energy;
+
+    energy_condition(const double min_energy)
+        : m_min_energy(min_energy) { }
+
+    double energy(const state_type &x) {
+        return 0.5 * x[1] * x[1] + 0.5 * x[0] * x[0];
+    }
+
+    bool operator()(const state_type &x) {
+        // becomes true if the energy becomes smaller than the threshold
+        return energy(x) <= m_min_energy;
+    }
+};
+
+
+template<class System, class Condition>
+std::pair<double, state_type>
+find_condition(state_type &x0, System sys, Condition cond,
+               const double t_start, const double t_end, const double dt,
+               const double precision = 1E-6) {
+
+    // integrates an ODE until some threshold is crossed
+    // returns time and state at the point of the threshold crossing
+    // if no threshold crossing is found, some time > t_end is returned
+
+    auto stepper = odeint::make_dense_output(1.0e-6, 1.0e-6,
+                                             odeint::runge_kutta_dopri5<state_type>());
+
+    auto ode_range = odeint::make_adaptive_range(std::ref(stepper), sys, x0,
+                                                 t_start, t_end, dt);
+
+    // find the step where the condition changes
+    auto found_iter = std::find_if(ode_range.first, ode_range.second, cond);
+
+    if(found_iter == ode_range.second)
+    {
+        // no threshold crossing -> return time after t_end and ic
+        return std::make_pair(t_end + dt, x0);
+    }
+
+    // the dense out stepper now covers the interval where the condition changes
+    // improve the solution by bisection
+    double t0 = stepper.previous_time();
+    double t1 = stepper.current_time();
+    double t_m;
+    state_type x_m;
+    // use odeint's resizing functionality to allocate memory for x_m
+    odeint::adjust_size_by_resizeability(x_m, x0,
+                                         typename odeint::is_resizeable<state_type>::type());
+    while(std::abs(t1 - t0) > precision) {
+        t_m = 0.5 * (t0 + t1);  // get the mid point time
+        stepper.calc_state(t_m, x_m); // obtain the corresponding state
+        if (cond(x_m))
+            t1 = t_m;  // condition changer lies before midpoint
+        else
+            t0 = t_m;  // condition changer lies after midpoint
+    }
+    // we found the interval of size eps, take it's midpoint as final guess
+    t_m = 0.5 * (t0 + t1);
+    stepper.calc_state(t_m, x_m);
+    return std::make_pair(t_m, x_m);
+}
+
+
+int main(int argc, char **argv)
+{
+    state_type x0 = {{10.0, 0.0}};
+    const double t_start = 0.0;
+    const double t_end = 10.0;
+    const double dt = 0.1;
+    const double threshold = 0.1;
+
+    energy_condition cond(threshold);
+    state_type x_cond;
+    double t_cond;
+    std::tie(t_cond, x_cond) = find_condition(x0, damped_osc, cond,
+                                              t_start, t_end, dt, 1E-6);
+    if(t_cond > t_end)
+    {
+        // time after t_end -> no threshold crossing within [t_start, t_end]
+        std::cout << "No threshold crossing found." << std::endl;
+    } else
+    {
+        std::cout.precision(16);
+        std::cout << "Time of energy threshold crossing: " << t_cond << std::endl;
+        std::cout << "State: [" << x_cond[0] << " , " << x_cond[1] << "]" << std::endl;
+        std::cout << "Energy: " << cond.energy(x_cond) << std::endl;
+    }
+}

libs/numeric/odeint/examples/fpu.cpp

@@ -0,0 +1,169 @@
+/*
+ * fpu.cpp
+ *
+ * This example demonstrates how one can use odeint to solve the Fermi-Pasta-Ulam system.
+
+ *  Created on: July 13, 2011
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2011 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)
+ */
+
+#include <iostream>
+#include <numeric>
+#include <cmath>
+#include <vector>
+
+#include <boost/numeric/odeint.hpp>
+
+#ifndef M_PI //not there on windows
+#define M_PI 3.1415927 //...
+#endif
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+//[ fpu_system_function
+typedef vector< double > container_type;
+
+struct fpu
+{
+    const double m_beta;
+
+    fpu( const double beta = 1.0 ) : m_beta( beta ) { }
+
+    // system function defining the ODE
+    void operator()( const container_type &q , container_type &dpdt ) const
+    {
+        size_t n = q.size();
+        double tmp = q[0] - 0.0;
+        double tmp2 = tmp + m_beta * tmp * tmp * tmp;
+        dpdt[0] = -tmp2;
+        for( size_t i=0 ; i<n-1 ; ++i )
+        {
+            tmp = q[i+1] - q[i];
+            tmp2 = tmp + m_beta * tmp * tmp * tmp;
+            dpdt[i] += tmp2;
+            dpdt[i+1] = -tmp2;
+        }
+        tmp = - q[n-1];
+        tmp2 = tmp + m_beta * tmp * tmp * tmp;
+        dpdt[n-1] += tmp2;
+    }
+
+    // calculates the energy of the system
+    double energy( const container_type &q , const container_type &p ) const
+    {
+        // ...
+        //<-
+        double energy = 0.0;
+        size_t n = q.size();
+
+        double tmp = q[0];
+        energy += 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
+        for( size_t i=0 ; i<n-1 ; ++i )
+        {
+            tmp = q[i+1] - q[i];
+            energy += 0.5 * ( p[i] * p[i] + tmp * tmp ) + 0.25 * m_beta * tmp * tmp * tmp * tmp;
+        }
+        energy += 0.5 * p[n-1] * p[n-1];
+        tmp = q[n-1];
+        energy += 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
+
+        return energy;
+        //->
+    }
+
+    // calculates the local energy of the system
+    void local_energy( const container_type &q , const container_type &p , container_type &e ) const
+    {
+        // ...
+        //<-
+        size_t n = q.size();
+        double tmp = q[0];
+        double tmp2 = 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
+        e[0] = tmp2;
+        for( size_t i=0 ; i<n-1 ; ++i )
+        {
+            tmp = q[i+1] - q[i];
+            tmp2 = 0.25 * tmp * tmp + 0.125 * m_beta * tmp * tmp * tmp * tmp;
+            e[i] += 0.5 * p[i] * p[i] + tmp2 ;
+            e[i+1] = tmp2;
+        }
+        tmp = q[n-1];
+        tmp2 = 0.5 * tmp * tmp + 0.25 * m_beta * tmp * tmp * tmp * tmp;
+        e[n-1] += 0.5 * p[n-1] * p[n-1] + tmp2;
+        //->
+    }
+};
+//]
+
+
+
+//[ fpu_observer
+struct streaming_observer
+{
+    std::ostream& m_out;
+    const fpu &m_fpu;
+    size_t m_write_every;
+    size_t m_count;
+
+    streaming_observer( std::ostream &out , const fpu &f , size_t write_every = 100 )
+    : m_out( out ) , m_fpu( f ) , m_write_every( write_every ) , m_count( 0 ) { }
+
+    template< class State >
+    void operator()( const State &x , double t )
+    {
+        if( ( m_count % m_write_every ) == 0 )
+        {
+            container_type &q = x.first;
+            container_type &p = x.second;
+            container_type energy( q.size() );
+            m_fpu.local_energy( q , p , energy );
+            for( size_t i=0 ; i<q.size() ; ++i )
+            {
+                m_out << t << "\t" << i << "\t" << q[i] << "\t" << p[i] << "\t" << energy[i] << "\n";
+            }
+            m_out << "\n";
+            clog << t << "\t" << accumulate( energy.begin() , energy.end() , 0.0 ) << "\n";
+        }
+        ++m_count;
+    }
+};
+//]
+
+
+
+
+
+
+
+
+int main( int argc , char **argv )
+{
+    //[ fpu_integration
+    const size_t n = 64;
+    container_type q( n , 0.0 ) , p( n , 0.0 );
+
+    for( size_t i=0 ; i<n ; ++i )
+    {
+        p[i] = 0.0;
+        q[i] = 32.0 * sin( double( i + 1 ) / double( n + 1 ) * M_PI );
+    }
+
+
+    const double dt = 0.1;
+
+    typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
+    fpu fpu_instance( 8.0 );
+
+    integrate_const( stepper_type() , fpu_instance ,
+            make_pair( boost::ref( q ) , boost::ref( p ) ) ,
+            0.0 , 1000.0 , dt , streaming_observer( cout , fpu_instance , 10 ) );
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/generation_functions.cpp

@@ -0,0 +1,113 @@
+/*
+ libs/numeric/odeint/examples/stochastic_euler.hpp
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012-2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Stochastic euler stepper example and Ornstein-Uhlenbeck process
+
+ 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)
+ */
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+typedef boost::array< double , 1 > state_type;
+
+using namespace boost::numeric::odeint;
+
+
+//[ generation_functions_own_steppers
+class custom_stepper
+{
+public:
+
+    typedef double value_type;
+    // ...
+};
+
+class custom_controller
+{
+    // ...
+};
+
+class custom_dense_output
+{
+    // ...
+};
+//]
+
+
+//[ generation_functions_get_controller
+namespace boost { namespace numeric { namespace odeint {
+
+template<>
+struct get_controller< custom_stepper >
+{
+    typedef custom_controller type;
+};
+
+} } }
+//]
+
+//[ generation_functions_controller_factory
+namespace boost { namespace numeric { namespace odeint {
+
+template<>
+struct controller_factory< custom_stepper , custom_controller >
+{
+    custom_controller operator()( double abs_tol , double rel_tol , const custom_stepper & ) const
+    {
+        return custom_controller();
+    }
+
+    custom_controller operator()( double abs_tol , double rel_tol , double max_dt ,
+                                  const custom_stepper & ) const
+    {
+        // version with maximal allowed step size max_dt
+        return custom_controller();
+    }
+};
+
+} } }
+//]
+
+int main( int argc , char **argv )
+{
+    {
+        typedef runge_kutta_dopri5< state_type > stepper_type;
+
+        /*
+        //[ generation_functions_syntax_auto
+        auto stepper1 = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() );
+        // or with max step size limit:
+        // auto stepper1 = make_controlled( 1.0e-6 , 1.0e-6 , 0.01, stepper_type() );
+
+        auto stepper2 = make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() );
+        //]
+        */
+
+        //[ generation_functions_syntax_result_of
+        boost::numeric::odeint::result_of::make_controlled< stepper_type >::type stepper3 = make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() );
+        (void)stepper3;
+        boost::numeric::odeint::result_of::make_dense_output< stepper_type >::type stepper4 = make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() );
+        (void)stepper4;
+        //]
+    }
+
+    {
+        /*
+        //[ generation_functions_example_custom_controller
+        auto stepper5 = make_controlled( 1.0e-6 , 1.0e-6 , custom_stepper() );
+        //]
+        */
+
+        boost::numeric::odeint::result_of::make_controlled< custom_stepper >::type stepper5 = make_controlled( 1.0e-6 , 1.0e-6 , custom_stepper() );
+        (void)stepper5;
+    }
+    return 0;
+}

libs/numeric/odeint/examples/gmpxx/Makefile

@@ -0,0 +1,5 @@
+CXXFLAGS = -I${BOOST_ROOT} -O2 -static
+LDFLAGS = -lgmpxx -lgmp 
+
+lorenz_gmpxx: lorenz_gmpxx.cpp
+	${CXX} ${CXXFLAGS} lorenz_gmpxx.cpp ${LDFLAGS} -o lorenz_gmpxx

libs/numeric/odeint/examples/gmpxx/lorenz_gmpxx.cpp

@@ -0,0 +1,83 @@
+/*
+ * lorenz_gmpxx.cpp
+ *
+ * This example demonstrates how odeint can be used with arbitrary precision types.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2011-2012 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)
+ */
+
+
+
+#include <iostream>
+#include <boost/array.hpp>
+
+#include <gmpxx.h>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+//[ gmpxx_lorenz
+typedef mpf_class value_type;
+typedef boost::array< value_type , 3 > state_type;
+
+struct lorenz
+{
+    void operator()( const state_type &x , state_type &dxdt , value_type t ) const
+    {
+        const value_type sigma( 10.0 );
+        const value_type R( 28.0 );
+        const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
+
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+//]
+
+
+
+
+struct streaming_observer
+{
+    std::ostream& m_out;
+
+    streaming_observer( std::ostream &out ) : m_out( out ) { }
+
+    template< class State , class Time >
+    void operator()( const State &x , Time t ) const
+    {
+        m_out << t;
+        for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i] ;
+        m_out << "\n";
+    }
+};
+
+
+
+
+
+
+int main( int argc , char **argv )
+{
+    //[ gmpxx_integration
+    const int precision = 1024;
+    mpf_set_default_prec( precision );
+
+    state_type x = {{ value_type( 10.0 ) , value_type( 10.0 ) , value_type( 10.0 ) }};
+
+    cout.precision( 1000 );
+    integrate_const( runge_kutta4< state_type , value_type >() ,
+            lorenz() , x , value_type( 0.0 ) , value_type( 10.0 ) , value_type( value_type( 1.0 ) / value_type( 10.0 ) ) ,
+            streaming_observer( cout ) );
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/gram_schmidt.hpp

@@ -0,0 +1,89 @@
+/*
+ boost header: numeric/odeint/gram_schmitt.hpp
+
+ Copyright 2011-2013 Karsten Ahnert
+ Copyright 2011 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)
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
+
+#include <boost/throw_exception.hpp>
+#include <iterator>
+#include <algorithm>
+#include <numeric>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+template< class Iterator , class T >
+void normalize( Iterator first , Iterator last , T norm )
+{
+    while( first != last ) *first++ /= norm;
+}
+
+template< class Iterator , class T >
+void substract_vector( Iterator first1 , Iterator last1 ,
+        Iterator first2 , T val )
+{
+    while( first1 != last1 ) *first1++ -= val * ( *first2++ );
+}
+
+template< size_t num_of_lyap , class StateType , class LyapType >
+void gram_schmidt( StateType &x , LyapType &lyap , size_t n )
+{
+    if( !num_of_lyap ) return;
+    if( ptrdiff_t( ( num_of_lyap + 1 ) * n ) != std::distance( x.begin() , x.end() ) )
+        BOOST_THROW_EXCEPTION( std::domain_error( "renormalization() : size of state does not match the number of lyapunov exponents." ) );
+
+    typedef typename StateType::value_type value_type;
+    typedef typename StateType::iterator iterator;
+
+    value_type norm[num_of_lyap];
+    value_type tmp[num_of_lyap];
+    iterator first = x.begin() + n;
+    iterator beg1 = first , end1 = first + n ;
+
+    std::fill( norm , norm+num_of_lyap , 0.0 );
+
+    // normalize first vector
+    norm[0] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
+    normalize( beg1 , end1 , norm[0] );
+
+    beg1 += n;
+    end1 += n;
+
+    for( size_t j=1 ; j<num_of_lyap ; ++j , beg1+=n , end1+=n )
+    {
+        for( size_t k=0 ; k<j ; ++k )
+        {
+            tmp[k] = std::inner_product( beg1 , end1 , first + k*n , 0.0 );
+            //  clog << j << " " << k << " " << tmp[k] << "\n";
+        }
+
+
+
+        for( size_t k=0 ; k<j ; ++k )
+            substract_vector( beg1 , end1 , first + k*n , tmp[k] );
+
+        // normalize j-th vector
+        norm[j] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
+        // clog << j << " " << norm[j] << "\n";
+        normalize( beg1 , end1 , norm[j] );
+    }
+
+    for( size_t j=0 ; j<num_of_lyap ; j++ )
+        lyap[j] += log( norm[j] );
+}
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif //BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED

libs/numeric/odeint/examples/harmonic_oscillator.cpp

@@ -0,0 +1,211 @@
+/*
+ Copyright 2010-2012 Karsten Ahnert
+ Copyright 2011-2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ 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)
+ */
+
+
+#include <iostream>
+#include <vector>
+
+#include <boost/numeric/odeint.hpp>
+
+
+
+//[ rhs_function
+/* The type of container used to hold the state vector */
+typedef std::vector< double > state_type;
+
+const double gam = 0.15;
+
+/* The rhs of x' = f(x) */
+void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
+{
+    dxdt[0] = x[1];
+    dxdt[1] = -x[0] - gam*x[1];
+}
+//]
+
+
+
+
+
+//[ rhs_class
+/* The rhs of x' = f(x) defined as a class */
+class harm_osc {
+
+    double m_gam;
+
+public:
+    harm_osc( double gam ) : m_gam(gam) { }
+
+    void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
+    {
+        dxdt[0] = x[1];
+        dxdt[1] = -x[0] - m_gam*x[1];
+    }
+};
+//]
+
+
+
+
+
+//[ integrate_observer
+struct push_back_state_and_time
+{
+    std::vector< state_type >& m_states;
+    std::vector< double >& m_times;
+
+    push_back_state_and_time( std::vector< state_type > &states , std::vector< double > &times )
+    : m_states( states ) , m_times( times ) { }
+
+    void operator()( const state_type &x , double t )
+    {
+        m_states.push_back( x );
+        m_times.push_back( t );
+    }
+};
+//]
+
+struct write_state
+{
+    void operator()( const state_type &x ) const
+    {
+        std::cout << x[0] << "\t" << x[1] << "\n";
+    }
+};
+
+
+int main(int /* argc */ , char** /* argv */ )
+{
+    using namespace std;
+    using namespace boost::numeric::odeint;
+
+
+    //[ state_initialization
+    state_type x(2);
+    x[0] = 1.0; // start at x=1.0, p=0.0
+    x[1] = 0.0;
+    //]
+
+
+
+    //[ integration
+    size_t steps = integrate( harmonic_oscillator ,
+            x , 0.0 , 10.0 , 0.1 );
+    //]
+
+
+
+    //[ integration_class
+    harm_osc ho(0.15);
+    steps = integrate( ho ,
+            x , 0.0 , 10.0 , 0.1 );
+    //]
+
+
+
+
+
+    //[ integrate_observ
+    vector<state_type> x_vec;
+    vector<double> times;
+
+    steps = integrate( harmonic_oscillator ,
+            x , 0.0 , 10.0 , 0.1 ,
+            push_back_state_and_time( x_vec , times ) );
+
+    /* output */
+    for( size_t i=0; i<=steps; i++ )
+    {
+        cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
+    }
+    //]
+
+
+
+
+
+
+
+    //[ define_const_stepper
+    runge_kutta4< state_type > stepper;
+    integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
+    //]
+
+
+
+
+    //[ integrate_const_loop
+    const double dt = 0.01;
+    for( double t=0.0 ; t<10.0 ; t+= dt )
+        stepper.do_step( harmonic_oscillator , x , t , dt );
+    //]
+
+
+
+
+    //[ define_adapt_stepper
+    typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
+    //]
+
+
+
+    //[ integrate_adapt
+    typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
+    controlled_stepper_type controlled_stepper;
+    integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
+    //]
+
+    {
+    //[integrate_adapt_full
+    double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
+    controlled_stepper_type controlled_stepper( 
+        default_error_checker< double , range_algebra , default_operations >( abs_err , rel_err , a_x , a_dxdt ) );
+    integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
+    //]
+    }
+
+
+    //[integrate_adapt_make_controlled
+    integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) , 
+                        harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
+    //]
+
+
+
+
+    //[integrate_adapt_make_controlled_alternative
+    integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) , 
+                        harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
+    //]
+
+    #ifdef BOOST_NUMERIC_ODEINT_CXX11
+    //[ define_const_stepper_cpp11
+    {
+    runge_kutta4< state_type > stepper;
+    integrate_const( stepper , []( const state_type &x , state_type &dxdt , double t ) {
+            dxdt[0] = x[1]; dxdt[1] = -x[0] - gam*x[1]; }
+        , x , 0.0 , 10.0 , 0.01 );
+    }
+    //]
+    
+    
+    
+    //[ harm_iterator_const_step]
+    std::for_each( make_const_step_time_iterator_begin( stepper , harmonic_oscillator, x , 0.0 , 0.1 , 10.0 ) ,
+                   make_const_step_time_iterator_end( stepper , harmonic_oscillator, x ) ,
+                   []( std::pair< const state_type & , const double & > x ) {
+                       cout << x.second << " " << x.first[0] << " " << x.first[1] << "\n"; } );
+    //]
+    #endif
+    
+    
+
+
+}

libs/numeric/odeint/examples/harmonic_oscillator_units.cpp

@@ -0,0 +1,121 @@
+/*
+ Copyright 2011-2013 Karsten Ahnert
+ Copyright 2011-2013 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)
+ */
+
+
+#include <iostream>
+#include <vector>
+
+
+/* WARNING: Compilation in debug mode might consume enormous memory 
+   (e.g. ~2GB on gcc 4.4 ) 
+*/
+
+// first increase fusion macro variables (now done by odeint itself)
+//#define BOOST_FUSION_INVOKE_MAX_ARITY 15
+//#define BOOST_RESULT_OF_NUM_ARGS 15
+
+//[ units_define_basic_quantities
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/algebra/fusion_algebra.hpp>
+#include <boost/numeric/odeint/algebra/fusion_algebra_dispatcher.hpp>
+
+#include <boost/units/systems/si/length.hpp>
+#include <boost/units/systems/si/time.hpp>
+#include <boost/units/systems/si/velocity.hpp>
+#include <boost/units/systems/si/acceleration.hpp>
+#include <boost/units/systems/si/io.hpp>
+
+#include <boost/fusion/container.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+namespace fusion = boost::fusion;
+namespace units = boost::units;
+namespace si = boost::units::si;
+
+typedef units::quantity< si::time , double > time_type;
+typedef units::quantity< si::length , double > length_type;
+typedef units::quantity< si::velocity , double > velocity_type;
+typedef units::quantity< si::acceleration , double > acceleration_type;
+typedef units::quantity< si::frequency , double > frequency_type;
+
+typedef fusion::vector< length_type , velocity_type > state_type;
+typedef fusion::vector< velocity_type , acceleration_type > deriv_type;
+//]
+
+
+
+//[ units_define_ode
+struct oscillator
+{
+    frequency_type m_omega;
+
+    oscillator( const frequency_type &omega = 1.0 * si::hertz ) : m_omega( omega ) { }
+
+    void operator()( const state_type &x , deriv_type &dxdt , time_type t ) const
+    {
+        fusion::at_c< 0 >( dxdt ) = fusion::at_c< 1 >( x );
+        fusion::at_c< 1 >( dxdt ) = - m_omega * m_omega * fusion::at_c< 0 >( x );
+    }
+};
+//]
+
+
+//[ units_observer
+struct streaming_observer
+{
+    std::ostream& m_out;
+
+    streaming_observer( std::ostream &out ) : m_out( out ) { }
+
+    struct write_element
+    {
+        std::ostream &m_out;
+        write_element( std::ostream &out ) : m_out( out ) { };
+
+        template< class T >
+        void operator()( const T &t ) const
+        {
+            m_out << "\t" << t;
+        }
+    };
+
+    template< class State , class Time >
+    void operator()( const State &x , const Time &t ) const
+    {
+        m_out << t;
+        fusion::for_each( x , write_element( m_out ) );
+        m_out << "\n";
+    }
+};
+//]
+
+
+int main( int argc , char**argv )
+{
+//    typedef dense_output_runge_kutta
+//    <
+//        controlled_runge_kutta
+//        <
+//            runge_kutta_dopri5< state_type , double , deriv_type , time_type , fusion_algebra >
+//        >
+//    > stepper_type;
+
+    //[ units_define_stepper
+    typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type > stepper_type;
+
+    state_type x( 1.0 * si::meter , 0.0 * si::meter_per_second );
+
+    integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , oscillator( 2.0 * si::hertz ) ,
+                     x , 0.0 * si::second , 100.0 * si::second , 0.1 * si::second , streaming_observer( cout ) );
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/heun.cpp

@@ -0,0 +1,170 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/examples/heun.cpp
+
+ [begin_description]
+ Examplary implementation of the method of Heun.
+ [end_description]
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 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)
+ */
+
+#include <iostream>
+
+
+#include <boost/fusion/container/vector.hpp>
+#include <boost/fusion/container/generation/make_vector.hpp>
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+
+
+
+
+
+namespace fusion = boost::fusion;
+
+//[ heun_define_coefficients
+template< class Value = double >
+struct heun_a1 : boost::array< Value , 1 > {
+    heun_a1( void )
+    {
+        (*this)[0] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
+    }
+};
+
+template< class Value = double >
+struct heun_a2 : boost::array< Value , 2 >
+{
+    heun_a2( void )
+    {
+        (*this)[0] = static_cast< Value >( 0 );
+        (*this)[1] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
+    }
+};
+
+
+template< class Value = double >
+struct heun_b : boost::array< Value , 3 >
+{
+    heun_b( void )
+    {
+        (*this)[0] = static_cast<Value>( 1 ) / static_cast<Value>( 4 );
+        (*this)[1] = static_cast<Value>( 0 );
+        (*this)[2] = static_cast<Value>( 3 ) / static_cast<Value>( 4 );
+    }
+};
+
+template< class Value = double >
+struct heun_c : boost::array< Value , 3 >
+{
+    heun_c( void )
+    {
+        (*this)[0] = static_cast< Value >( 0 );
+        (*this)[1] = static_cast< Value >( 1 ) / static_cast< Value >( 3 );
+        (*this)[2] = static_cast< Value >( 2 ) / static_cast< Value >( 3 );
+    }
+};
+//]
+
+
+//[ heun_stepper_definition
+template<
+    class State ,
+    class Value = double ,
+    class Deriv = State ,
+    class Time = Value ,
+    class Algebra = boost::numeric::odeint::range_algebra ,
+    class Operations = boost::numeric::odeint::default_operations ,
+    class Resizer = boost::numeric::odeint::initially_resizer
+>
+class heun : public
+boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
+                                             Algebra , Operations , Resizer >
+{
+
+public:
+
+    typedef boost::numeric::odeint::explicit_generic_rk< 3 , 3 , State , Value , Deriv , Time ,
+                                                         Algebra , Operations , Resizer > stepper_base_type;
+
+    typedef typename stepper_base_type::state_type state_type;
+    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
+    typedef typename stepper_base_type::value_type value_type;
+    typedef typename stepper_base_type::deriv_type deriv_type;
+    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
+    typedef typename stepper_base_type::time_type time_type;
+    typedef typename stepper_base_type::algebra_type algebra_type;
+    typedef typename stepper_base_type::operations_type operations_type;
+    typedef typename stepper_base_type::resizer_type resizer_type;
+    typedef typename stepper_base_type::stepper_type stepper_type;
+
+    heun( const algebra_type &algebra = algebra_type() )
+    : stepper_base_type(
+            fusion::make_vector(
+                heun_a1<Value>() ,
+                heun_a2<Value>() ) ,
+            heun_b<Value>() , heun_c<Value>() , algebra )
+    { }
+};
+//]
+
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+struct lorenz
+{
+    template< class State , class Deriv >
+    void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
+    {
+        typename boost::range_iterator< const State >::type x = boost::begin( x_ );
+        typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );
+
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+
+struct streaming_observer
+{
+    std::ostream &m_out;
+    streaming_observer( std::ostream &out ) : m_out( out ) { }
+    template< typename State , typename Value >
+    void operator()( const State &x , Value t ) const
+    {
+        m_out << t;
+        for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i];
+        m_out << "\n";
+    }
+};
+
+
+
+int main( int argc , char **argv )
+{
+    using namespace std;
+    using namespace boost::numeric::odeint;
+
+
+    //[ heun_example
+    typedef boost::array< double , 3 > state_type;
+    heun< state_type > h;
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+
+    integrate_const( h , lorenz() , x , 0.0 , 100.0 , 0.01 ,
+                     streaming_observer( std::cout ) );
+
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/integrate_times.cpp

@@ -0,0 +1,54 @@
+/* Boost libs/numeric/odeint/examples/integrate_times.cpp
+
+ Copyright 2009-2014 Karsten Ahnert
+ Copyright 2009-2014 Mario Mulansky
+
+ example for the use of integrate_times
+
+ 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)
+ */
+
+
+#include <iostream>
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+/*
+ * simple 1D ODE
+ */
+
+void rhs( const double x , double &dxdt , const double t )
+{
+    dxdt = 3.0/(2.0*t*t) + x/(2.0*t);
+}
+
+void write_cout( const double &x , const double t )
+{
+    cout << t << '\t' << x << endl;
+}
+
+// state_type = double
+typedef runge_kutta_dopri5< double > stepper_type;
+
+const double dt = 0.1;
+
+int main()
+{
+    // create a vector with observation time points
+    std::vector<double> times( 91 );
+    for( size_t i=0 ; i<times.size() ; ++i )
+        times[i] = 1.0 + dt*i;
+
+    double x = 0.0; //initial value x(1) = 0
+    // we can provide the observation time as a boost range (i.e. the vector)
+    integrate_times( make_controlled( 1E-12 , 1E-12 , stepper_type() ) , rhs ,
+                     x , times , dt , write_cout );
+    // or as two iterators
+    //integrate_times( make_controlled( 1E-12 , 1E-12 , stepper_type() ) , rhs ,
+    //                   x , times.begin() , times.end() , dt , write_cout );
+}

libs/numeric/odeint/examples/list_lattice.cpp

@@ -0,0 +1,78 @@
+/*
+ Copyright 2011-2012 Mario Mulansky
+ Copyright 2012-2013 Karsten Ahnert
+
+ 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)
+ */
+
+
+/* example showing how odeint can be used with std::list */
+
+#include <iostream>
+#include <cmath>
+#include <list>
+
+#include <boost/numeric/odeint.hpp>
+
+//[ list_bindings
+typedef std::list< double > state_type;
+
+namespace boost { namespace numeric { namespace odeint {
+
+template< >
+struct is_resizeable< state_type >
+{ // declare resizeability
+    typedef boost::true_type type;
+    const static bool value = type::value;
+};
+
+template< > 
+struct same_size_impl< state_type , state_type >
+{ // define how to check size
+    static bool same_size( const state_type &v1 ,
+                           const state_type &v2 )
+    {
+        return v1.size() == v2.size();
+    }
+};
+
+template< >
+struct resize_impl< state_type , state_type >
+{ // define how to resize
+    static void resize( state_type &v1 ,
+                        const state_type &v2 )
+    {
+        v1.resize( v2.size() );
+    }
+};
+
+} } } 
+//]
+
+void lattice( const state_type &x , state_type &dxdt , const double /* t */ )
+{
+    state_type::const_iterator x_begin = x.begin();
+    state_type::const_iterator x_end = x.end();
+    state_type::iterator dxdt_begin = dxdt.begin();
+    
+    x_end--; // stop one before last
+    while( x_begin != x_end )
+    {
+        *(dxdt_begin++) = std::sin( *(x_begin) - *(x_begin++) );
+    }
+    *dxdt_begin = sin( *x_begin - *(x.begin()) ); // periodic boundary
+}
+
+using namespace boost::numeric::odeint;
+
+int main()
+{
+    const int N = 32;
+    state_type x;
+    for( int i=0 ; i<N ; ++i )
+        x.push_back( 1.0*i/N );
+
+    integrate_const( runge_kutta4< state_type >() , lattice , x , 0.0 , 10.0 , 0.1 );
+}

libs/numeric/odeint/examples/lorenz.cpp

@@ -0,0 +1,31 @@
+#include <iostream>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+typedef boost::array< double , 3 > state_type;
+
+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] = -b * x[2] + x[0] * x[1];
+}
+
+void write_lorenz( const state_type &x , const double t )
+{
+    cout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl;
+}
+
+int main(int argc, char **argv)
+{
+    state_type x = {{ 10.0 , 1.0 , 1.0 }}; // initial conditions
+    integrate( lorenz , x , 0.0 , 25.0 , 0.1 , write_lorenz );
+}

libs/numeric/odeint/examples/lorenz_point.cpp

@@ -0,0 +1,119 @@
+/*
+ * Copyright 2011-2013 Mario Mulansky
+ * Copyright 2012 Karsten Ahnert
+ *
+ * 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)
+ *
+ * Example for the lorenz system with a 3D point type
+*/
+
+#include <iostream>
+#include <cmath>
+
+#include <boost/operators.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+
+//[point3D
+class point3D :
+    boost::additive1< point3D ,
+    boost::additive2< point3D , double ,
+    boost::multiplicative2< point3D , double > > >
+{
+public:
+
+    double x , y , z;
+
+    point3D()
+        : x( 0.0 ) , y( 0.0 ) , z( 0.0 )
+    { }
+
+    point3D( const double val )
+        : x( val ) , y( val ) , z( val )
+    { }
+
+    point3D( const double _x , const double _y , const double _z  )
+        : x( _x ) , y( _y ) , z( _z )
+    { }
+
+    point3D& operator+=( const point3D &p )
+    {
+        x += p.x; y += p.y; z += p.z;
+        return *this;
+    }
+
+    point3D& operator*=( const double a )
+    {
+        x *= a; y *= a; z *= a;
+        return *this;
+    }
+
+};
+//]
+
+//[point3D_abs_div
+// only required for steppers with error control
+point3D operator/( const point3D &p1 , const point3D &p2 )
+{
+    return point3D( p1.x/p2.x , p1.y/p2.y , p1.z/p2.z );
+}
+
+point3D abs( const point3D &p )
+{
+    return point3D( std::abs(p.x) , std::abs(p.y) , std::abs(p.z) );
+}
+//]
+
+//[point3D_norm
+// also only for steppers with error control
+namespace boost { namespace numeric { namespace odeint {
+template<>
+struct vector_space_norm_inf< point3D >
+{
+    typedef double result_type;
+    double operator()( const point3D &p ) const
+    {
+        using std::max;
+        using std::abs;
+        return max( max( abs( p.x ) , abs( p.y ) ) , abs( p.z ) );
+    }
+};
+} } }
+//]
+
+std::ostream& operator<<( std::ostream &out , const point3D &p )
+{
+    out << p.x << " " << p.y << " " << p.z;
+    return out;
+}
+
+//[point3D_main
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const point3D &x , point3D &dxdt , const double t )
+{
+    dxdt.x = sigma * ( x.y - x.x );
+    dxdt.y = R * x.x - x.y - x.x * x.z;
+    dxdt.z = -b * x.z + x.x * x.y;
+}
+
+using namespace boost::numeric::odeint;
+
+int main()
+{
+
+    point3D x( 10.0 , 5.0 , 5.0 );
+    // point type defines it's own operators -> use vector_space_algebra !
+    typedef runge_kutta_dopri5< point3D , double , point3D ,
+                                double , vector_space_algebra > stepper;
+    int steps = integrate_adaptive( make_controlled<stepper>( 1E-10 , 1E-10 ) , lorenz , x ,
+                                    0.0 , 10.0 , 0.1 );
+    std::cout << x << std::endl;
+    std::cout << "steps: " << steps << std::endl;
+}
+//]

libs/numeric/odeint/examples/molecular_dynamics.cpp

@@ -0,0 +1,160 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/examples/molecular_dynamics.cpp
+
+ [begin_description]
+ Molecular dynamics example.
+ [end_description]
+
+ Copyright 2009-2012 Karsten Ahnert
+ Copyright 2009-2012 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)
+ */
+
+#include <boost/numeric/odeint.hpp>
+
+#include <vector>
+#include <iostream>
+#include <random>
+
+using namespace boost::numeric::odeint;
+
+
+
+using namespace std;
+#define tab "\t"
+
+const size_t n1 = 16;
+const size_t n2 = 16;
+
+struct md_system
+{
+    static const size_t n = n1 * n2;
+    typedef std::vector< double > vector_type;
+
+    md_system( double a = 0.0 ,            // strength of harmonic oscillator
+               double gamma = 0.0   ,       // friction
+               double eps = 0.1 ,          // interaction strenght
+               double sigma = 1.0 ,            // interaction radius
+               double xmax = 150.0 , double ymax = 150.0 )
+    : m_a( a ) , m_gamma( gamma ) 
+    , m_eps( eps ) , m_sigma( sigma ) 
+    , m_xmax( xmax ) , m_ymax( ymax )
+    { }
+    
+    static void init_vector_type( vector_type &x ) { x.resize( 2 * n ); }
+    
+    void operator()( vector_type const& x , vector_type const& v , vector_type &a , double t ) const
+    {
+        for( size_t i=0 ; i<n ; ++i )
+        {
+            double diffx = x[i] - 0.5 * m_xmax , diffy = x[i+n] - 0.5 * m_ymax;
+            double r2 = diffx * diffx + diffy * diffy ;
+            double r = std::sqrt( r2 );
+            a[     i ] = - m_a * r * diffx - m_gamma * v[     i ] ;
+            a[ n + i ] = - m_a * r * diffy - m_gamma * v[ n + i ] ;
+        }
+        
+        for( size_t i=0 ; i<n ; ++i )
+        {
+            double xi = x[i] , yi = x[n+i];
+            xi = periodic_bc( xi , m_xmax );
+            yi = periodic_bc( yi , m_ymax );
+            for( size_t j=0 ; j<i ; ++j )
+            {
+                double xj = x[j] , yj = x[n+j];
+                xj = periodic_bc( xj , m_xmax );
+                yj = periodic_bc( yj , m_ymax );
+                
+                double diffx = ( xj - xi ) , diffy = ( yj - yi );
+                double r = sqrt( diffx * diffx + diffy * diffy );
+                double f = lennard_jones( r );
+                a[     i ] += diffx / r * f;
+                a[ n + i ] += diffy / r * f;
+                a[     j ] -= diffx / r * f;
+                a[ n + j ] -= diffy / r * f;
+            }
+        }
+    }
+    
+    void bc( vector_type &x )
+    {
+        for( size_t i=0 ; i<n ; ++i )
+        {
+            x[ i     ] = periodic_bc( x[ i     ] , m_xmax );
+            x[ i + n ] = periodic_bc( x[ i + n ] , m_ymax );
+        }
+    }
+    
+    inline double lennard_jones( double r ) const
+    {
+        double c = m_sigma / r;
+        double c3 = c * c * c;
+        double c6 = c3 * c3;
+        return 4.0 * m_eps * ( -12.0 * c6 * c6 / r + 6.0 * c6 / r );
+    }
+    
+    static inline double periodic_bc( double x , double xmax )
+    {
+        return ( x < 0.0 ) ? x + xmax : ( x > xmax ) ? x - xmax : x ;
+    }
+    
+    double m_a;
+    double m_gamma;
+    double m_eps ;
+    double m_sigma ;
+    double m_xmax , m_ymax;
+};
+
+
+
+
+
+int main( int argc , char *argv[] )
+{
+    const size_t n = md_system::n;
+    typedef md_system::vector_type vector_type;
+    
+    
+    std::mt19937 rng;
+    std::normal_distribution<> dist( 0.0 , 1.0 );
+    
+    vector_type x , v;
+    md_system::init_vector_type( x );
+    md_system::init_vector_type( v );
+    
+    for( size_t i=0 ; i<n1 ; ++i )
+    {
+        for( size_t j=0 ; j<n2 ; ++j )
+        {
+            x[i*n2+j  ] = 5.0 + i * 4.0 ;
+            x[i*n2+j+n] = 5.0 + j * 4.0 ;
+            v[i]   = dist( rng ) ;
+            v[i+n] = dist( rng ) ;
+        }
+    }
+    
+    velocity_verlet< vector_type > stepper;
+    const double dt = 0.025;
+    double t = 0.0;
+    md_system sys;
+    for( size_t oi=0 ; oi<100000 ; ++oi )
+    {
+        for( size_t ii=0 ; ii<100 ; ++ii,t+=dt )
+            stepper.do_step( sys , std::make_pair( std::ref( x ) , std::ref( v ) ) , t , dt );
+        sys.bc( x );
+        
+        std::cout << "set size square" << "\n";
+        std::cout << "unset key" << "\n";
+        std::cout << "p [0:" << sys.m_xmax << "][0:" << sys.m_ymax << "] '-' pt 7 ps 0.5" << "\n";
+        for( size_t i=0 ; i<n ; ++i )
+            std::cout << x[i] << " " << x[i+n] << " " << v[i] << " " << v[i+n] << "\n";
+        std::cout << "e" << std::endl;
+    }
+    
+    
+    return 0;
+}

libs/numeric/odeint/examples/molecular_dynamics_cells.cpp

@@ -0,0 +1,377 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/examples/molecular_dynamics_cells.cpp
+
+ [begin_description]
+ Molecular dynamics example with cells.
+ [end_description]
+
+ Copyright 2009-2012 Karsten Ahnert
+ Copyright 2009-2012 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)
+ */
+
+#include <boost/numeric/odeint.hpp>
+
+#include <cstddef>
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <tuple>
+#include <iostream>
+#include <random>
+
+#include <boost/range/algorithm/for_each.hpp>
+#include <boost/range/algorithm/sort.hpp>
+#include <boost/range/algorithm/unique_copy.hpp>
+#include <boost/range/algorithm_ext/iota.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/iterator/permutation_iterator.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+
+#include "point_type.hpp"
+
+
+
+
+
+
+
+
+struct local_force
+{
+    double m_gamma;        // friction
+    local_force( double gamma = 0.0 ) : m_gamma( gamma ) { }
+    template< typename Point >
+    Point operator()( Point& x , Point& v ) const
+    {
+        return - m_gamma * v;
+    }
+};
+
+
+struct lennard_jones
+{
+    double m_sigma;
+    double m_eps;
+    lennard_jones( double sigma = 1.0 , double eps = 0.1 ) : m_sigma( sigma ) , m_eps( eps ) { }
+    double operator()( double r ) const
+    {
+        double c = m_sigma / r;
+        double c3 = c * c * c;
+        double c6 = c3 * c3;
+        return 4.0 * m_eps * ( -12.0 * c6 * c6 / r + 6.0 * c6 / r );
+    }
+};
+
+template< typename F >
+struct conservative_interaction
+{
+    F m_f;
+    conservative_interaction( F const &f = F() ) : m_f( f ) { }
+    template< typename Point >
+    Point operator()( Point const& x1 , Point const& x2 ) const
+    {
+        Point diff = x1 - x2;
+        double r = abs( diff );
+        double f = m_f( r );
+        return - diff / r * f; 
+    }
+};
+
+template< typename F >
+conservative_interaction< F > make_conservative_interaction( F const &f )
+{
+    return conservative_interaction< F >( f );
+}
+
+
+
+
+
+
+
+// force = interaction( x1 , x2 )
+// force = local_force( x , v )
+template< typename LocalForce , typename Interaction >
+class md_system_bs
+{
+public:
+    
+    typedef std::vector< double > vector_type;
+    typedef point< double , 2 > point_type;
+    typedef point< int , 2 > index_type;
+    typedef std::vector< point_type > point_vector;
+    typedef std::vector< index_type > index_vector;
+    typedef std::vector< size_t > hash_vector;
+    typedef LocalForce local_force_type;
+    typedef Interaction interaction_type;
+
+
+    struct params
+    {
+        size_t n;
+        size_t n_cell_x , n_cell_y , n_cells;
+        double x_max , y_max , cell_size;
+        double eps , sigma;    // interaction strength, interaction radius
+        interaction_type interaction;
+        local_force_type local_force;
+    };
+
+
+    struct cell_functor
+    {
+        params const &m_p;
+
+        cell_functor( params const& p ) : m_p( p ) { }
+        
+        template< typename Tuple >
+        void operator()( Tuple const& t ) const
+        {
+            auto point = boost::get< 0 >( t );
+            size_t i1 = size_t( point[0] / m_p.cell_size ) , i2 = size_t( point[1] / m_p.cell_size );
+            boost::get< 1 >( t ) = index_type( i1 , i2 );
+            boost::get< 2 >( t ) = hash_func( boost::get< 1 >( t ) , m_p );
+        }
+    };
+
+
+
+    struct transform_functor : public std::unary_function< size_t , size_t >
+    {
+        hash_vector const* m_index;
+        transform_functor( hash_vector const& index ) : m_index( &index ) { }
+        size_t operator()( size_t i ) const { return (*m_index)[i]; }
+    };
+
+
+
+    struct interaction_functor
+    {
+        hash_vector const &m_cells_begin;
+        hash_vector const &m_cells_end;
+        hash_vector const &m_order;
+        point_vector const &m_x;
+        point_vector const &m_v;
+        params const &m_p;
+        size_t m_ncellx , m_ncelly;
+        
+        interaction_functor( hash_vector const& cells_begin , hash_vector const& cells_end , hash_vector pos_order ,
+                            point_vector const&x , point_vector const& v , params const &p )
+        : m_cells_begin( cells_begin ) , m_cells_end( cells_end ) , m_order( pos_order ) , m_x( x ) , m_v( v ) ,
+        m_p( p ) { }
+        
+        template< typename Tuple >
+        void operator()( Tuple const &t ) const
+        {
+            point_type x = periodic_bc( boost::get< 0 >( t ) , m_p ) , v = boost::get< 1 >( t );
+            index_type index = boost::get< 3 >( t );
+            size_t pos_hash = boost::get< 4 >( t );
+
+            point_type a = m_p.local_force( x , v );
+
+            for( int i=-1 ; i<=1 ; ++i )
+            {
+                for( int j=-1 ; j<=1 ; ++j )
+                {
+                    index_type cell_index = index + index_type( i , j );
+                    size_t cell_hash = hash_func( cell_index , m_p );
+                    for( size_t ii = m_cells_begin[ cell_hash ] ; ii < m_cells_end[ cell_hash ] ; ++ii )
+                    {
+                        if( m_order[ ii ] == pos_hash ) continue;
+                        point_type x2 = periodic_bc( m_x[ m_order[ ii ] ] , m_p );
+
+                        if( cell_index[0] >= m_p.n_cell_x ) x2[0] += m_p.x_max;
+                        if( cell_index[0] < 0 ) x2[0] -= m_p.x_max;
+                        if( cell_index[1] >= m_p.n_cell_y ) x2[1] += m_p.y_max;
+                        if( cell_index[1] < 0 ) x2[1] -= m_p.y_max;
+
+                        a += m_p.interaction( x , x2 );
+                    }
+                }
+            }
+            boost::get< 2 >( t ) = a;
+        }
+    };
+
+
+
+
+    md_system_bs( size_t n ,
+                  local_force_type const& local_force = local_force_type() ,                  
+                  interaction_type const& interaction = interaction_type() ,
+                  double xmax = 100.0 , double ymax = 100.0 , double cell_size = 2.0 )
+    : m_p() 
+    {
+        m_p.n = n;
+        m_p.x_max = xmax;
+        m_p.y_max = ymax;
+        m_p.interaction = interaction;
+        m_p.local_force = local_force;
+        m_p.n_cell_x = size_t( xmax / cell_size );
+        m_p.n_cell_y = size_t( ymax / cell_size );
+        m_p.n_cells = m_p.n_cell_x * m_p.n_cell_y;
+        m_p.cell_size = cell_size;
+    }
+    
+    void init_point_vector( point_vector &x ) const { x.resize( m_p.n ); }
+    
+    void operator()( point_vector const& x , point_vector const& v , point_vector &a , double t ) const
+    {
+        // init
+        hash_vector pos_hash( m_p.n , 0 );
+        index_vector pos_index( m_p.n );
+        hash_vector pos_order( m_p.n , 0 );
+        hash_vector cells_begin( m_p.n_cells ) , cells_end( m_p.n_cells ) , cell_order( m_p.n_cells );
+
+        boost::iota( pos_order , 0 );
+        boost::iota( cell_order , 0 );
+
+        // calculate grid hash
+        // calcHash( m_dGridParticleHash, m_dGridParticleIndex, dPos, m_numParticles);
+        std::for_each(
+            boost::make_zip_iterator( boost::make_tuple( x.begin() , pos_index.begin() , pos_hash.begin() ) ) ,
+            boost::make_zip_iterator( boost::make_tuple( x.end() , pos_index.end() , pos_hash.end() ) ) ,
+            cell_functor( m_p ) );
+
+//         // sort particles based on hash
+//         // sortParticles(m_dGridParticleHash, m_dGridParticleIndex, m_numParticles);        
+        boost::sort( pos_order , [&]( size_t i1 , size_t i2 ) -> bool {
+            return pos_hash[i1] < pos_hash[i2]; } );
+
+
+        
+        // reorder particle arrays into sorted order and find start and end of each cell
+        std::for_each( cell_order.begin() , cell_order.end() , [&]( size_t i ) {
+            auto pos_begin = boost::make_transform_iterator( pos_order.begin() , transform_functor( pos_hash ) );
+            auto pos_end = boost::make_transform_iterator( pos_order.end() , transform_functor( pos_hash ) );
+            cells_begin[ i ] = std::distance( pos_begin , std::lower_bound( pos_begin , pos_end , i ) );
+            cells_end[ i ] = std::distance( pos_begin , std::upper_bound( pos_begin , pos_end , i ) );
+        } );
+        
+        std::for_each(
+            boost::make_zip_iterator( boost::make_tuple(
+                x.begin() ,
+                v.begin() ,
+                a.begin() ,
+                pos_index.begin() ,
+                boost::counting_iterator< size_t >( 0 )
+            ) ) ,
+            boost::make_zip_iterator( boost::make_tuple(
+                x.end() ,
+                v.end() ,
+                a.end() ,
+                pos_index.end() ,
+                boost::counting_iterator< size_t >( m_p.n )
+            ) ) ,
+            interaction_functor( cells_begin , cells_end , pos_order , x , v , m_p ) );
+    }
+
+    void bc( point_vector &x )
+    {
+        for( size_t i=0 ; i<m_p.n ; ++i )
+        {
+            x[i][0] = periodic_bc( x[ i ][0] , m_p.x_max );
+            x[i][1] = periodic_bc( x[ i ][1] , m_p.y_max );
+        }
+    }
+    
+    static inline double periodic_bc( double x , double xmax )
+    {
+        double tmp = x - xmax * int( x / xmax );
+        return tmp >= 0.0 ? tmp : tmp + xmax;
+    }
+
+
+    static inline point_type periodic_bc( point_type const& x , params const& p ) 
+    {
+        return point_type( periodic_bc( x[0] , p.x_max ) , periodic_bc( x[1] , p.y_max ) );
+    }
+
+
+    static inline int check_interval( int i , int max )
+    {
+        int tmp = i % max;
+        return tmp >= 0 ? tmp : tmp + max;
+    }
+
+
+    static inline size_t hash_func( index_type index , params const & p )
+    {
+        size_t i1 = check_interval( index[0] , p.n_cell_x );
+        size_t i2 = check_interval( index[1] , p.n_cell_y );
+        return i1 * p.n_cell_y + i2;
+    }
+    
+    params m_p;
+};
+
+
+template< typename LocalForce , typename Interaction >
+md_system_bs< LocalForce , Interaction > make_md_system_bs( size_t n , LocalForce const &f1 , Interaction const &f2 ,
+    double xmax = 100.0 , double ymax = 100.0 , double cell_size = 2.0 )
+{
+    return md_system_bs< LocalForce , Interaction >( n , f1 , f2 , xmax , ymax , cell_size );
+}
+
+
+
+
+
+
+using namespace boost::numeric::odeint;
+
+
+
+int main( int argc , char *argv[] )
+{
+    const size_t n1 = 32;
+    const size_t n2 = 32;
+    const size_t n = n1 * n2;
+    auto sys = make_md_system_bs( n , local_force() , make_conservative_interaction( lennard_jones() ) , 100.0 , 100.0 , 2.0 );
+    typedef decltype( sys ) system_type;
+    typedef system_type::point_vector point_vector;
+    
+    std::mt19937 rng;
+    std::normal_distribution<> dist( 0.0 , 1.0 );
+    
+    point_vector x , v;
+    sys.init_point_vector( x );
+    sys.init_point_vector( v );
+    
+    for( size_t i=0 ; i<n1 ; ++i )
+    {
+        for( size_t j=0 ; j<n2 ; ++j )
+        {
+            size_t index = i * n2 + j; 
+            x[index][0] = 10.0 + i * 2.0 ;
+            x[index][1] = 10.0 + j * 2.0 ;
+            v[index][0] = dist( rng ) ;
+            v[index][1] = dist( rng ) ;
+        }
+    }
+    
+    velocity_verlet< point_vector > stepper;
+    const double dt = 0.025;
+    double t = 0.0;
+    // std::cout << "set term x11" << endl;
+    for( size_t oi=0 ; oi<10000 ; ++oi )
+    {
+        for( size_t ii=0 ; ii<50 ; ++ii,t+=dt )
+            stepper.do_step( sys , std::make_pair( std::ref( x ) , std::ref( v ) ) , t , dt );
+        sys.bc( x );
+        
+        std::cout << "set size square" << "\n";
+        std::cout << "unset key" << "\n";
+        std::cout << "p [0:" << sys.m_p.x_max << "][0:" << sys.m_p.y_max << "] '-' pt 7 ps 0.5" << "\n";
+        for( size_t i=0 ; i<n ; ++i )
+            std::cout << x[i][0] << " " << x[i][1] << " " << v[i][0] << " " << v[i][1] << "\n";
+        std::cout << "e" << std::endl;
+    }
+
+    return 0;
+}

libs/numeric/odeint/examples/mpi/Jamfile.v2

@@ -0,0 +1,15 @@
+# Copyright 2011-2013 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# Copyright 2013 Pascal Germroth
+# 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)
+
+project
+    : requirements
+      <define>BOOST_ALL_NO_LIB=1
+      <library>/boost//mpi
+      <library>/boost//timer
+    ;
+
+exe phase_chain : phase_chain.cpp ;

libs/numeric/odeint/examples/mpi/phase_chain.cpp

@@ -0,0 +1,120 @@
+/*
+ * phase_chain.cpp
+ *
+ * Example of MPI parallelization with odeint
+ *
+ * Copyright 2013 Karsten Ahnert
+ * Copyright 2013 Mario Mulansky
+ * Copyright 2013 Pascal Germroth
+ * 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)
+ */
+
+#include <iostream>
+#include <vector>
+#include <boost/random.hpp>
+#include <boost/timer/timer.hpp>
+//[phase_chain_mpi_header
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/mpi/mpi.hpp>
+//]
+
+using namespace std;
+using namespace boost::numeric::odeint;
+using boost::timer::cpu_timer;
+using boost::math::double_constants::pi;
+
+//[phase_chain_state
+typedef mpi_state< vector<double> > state_type;
+//]
+
+//[phase_chain_mpi_rhs
+struct phase_chain_mpi_state
+{
+    phase_chain_mpi_state( double gamma = 0.5 )
+    : m_gamma( gamma ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
+    {
+        const size_t M = x().size();
+        const bool have_left = x.world.rank() > 0,
+                   have_right = x.world.rank() < x.world.size()-1;
+        double x_left, x_right;
+        boost::mpi::request r_left, r_right;
+        if( have_left )
+        {
+            x.world.isend( x.world.rank()-1 , 0 , x().front() ); // send to x_right
+            r_left = x.world.irecv( x.world.rank()-1 , 0 , x_left ); // receive from x().back()
+        }
+        if( have_right )
+        {
+            x.world.isend( x.world.rank()+1 , 0 , x().back() ); // send to x_left
+            r_right = x.world.irecv( x.world.rank()+1 , 0 , x_right ); // receive from x().front()
+        }
+        for(size_t m = 1 ; m < M-1 ; ++m)
+        {
+            dxdt()[m] = coupling_func( x()[m+1] - x()[m] ) +
+                        coupling_func( x()[m-1] - x()[m] );
+        }
+        dxdt()[0] = coupling_func( x()[1] - x()[0] );
+        if( have_left )
+        {
+            r_left.wait();
+            dxdt()[0] += coupling_func( x_left - x().front() );
+        }
+        dxdt()[M-1] = coupling_func( x()[M-2] - x()[M-1] );
+        if( have_right )
+        {
+            r_right.wait();
+            dxdt()[M-1] += coupling_func( x_right - x().back() );
+        }
+    }
+
+    double coupling_func( double x ) const
+    {
+        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
+    }
+
+    double m_gamma;
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+    //[phase_chain_mpi_init
+    boost::mpi::environment env( argc , argv );
+    boost::mpi::communicator world;
+
+    const size_t N = 131101;
+    vector<double> x;
+    if( world.rank() == 0 )
+    {
+        x.resize( N );
+        boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
+        boost::random::mt19937 engine( 0 );
+        generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
+    }
+
+    state_type x_split( world );
+    split( x , x_split );
+    //]
+
+
+    cpu_timer timer;
+    //[phase_chain_mpi_integrate
+    integrate_n_steps( runge_kutta4<state_type>() , phase_chain_mpi_state( 1.2 ) ,
+                       x_split , 0.0 , 0.01 , 100 );
+    unsplit( x_split , x );
+    //]
+
+    if( world.rank() == 0 )
+    {
+        double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
+        std::cerr << run_time << "s" << std::endl;
+        // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
+    }
+
+    return 0;
+}

libs/numeric/odeint/examples/mtl/Jamfile.v2

@@ -0,0 +1,17 @@
+# Copyright 2011-2013 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# 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)
+
+# set your MTL4 directory here
+MTL4_INCLUDE = /home/mario/MTL4 ;
+
+project
+    : requirements 
+      <include>$(MTL4_INCLUDE)
+      <define>BOOST_ALL_NO_LIB=1
+    ;
+    
+exe gauss_packet : gauss_packet.cpp ;
+exe implicit_euler_mtl : implicit_euler_mtl.cpp ;

libs/numeric/odeint/examples/mtl/gauss_packet.cpp

@@ -0,0 +1,141 @@
+/*
+ * gauss_packet.cpp
+ *
+ * Schroedinger equation with potential barrier and periodic boundary conditions
+ * Initial Gauss packet moving to the right
+ *
+ * pipe output into gnuplot to see animation
+ *
+ * Implementation of Hamilton operator via MTL library
+ *
+ * Copyright 2011-2013 Mario Mulansky
+ * Copyright 2011-2012 Karsten Ahnert
+ *
+ * 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)
+ */
+
+
+#include <iostream>
+#include <complex>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/mtl4/mtl4.hpp>
+
+#include <boost/numeric/mtl/mtl.hpp>
+
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef mtl::dense_vector< complex< double > > state_type;
+
+struct hamiltonian {
+
+    typedef mtl::compressed2D< complex< double > > matrix_type;
+    matrix_type m_H;
+
+    hamiltonian( const int N ) : m_H( N , N )
+    {
+        // constructor with zero potential
+        m_H = 0.0;
+        initialize_kinetic_term();
+    }
+
+    //template< mtl::compressed2D< double > >
+    hamiltonian( mtl::compressed2D< double > &V ) : m_H( num_rows( V ) , num_cols( V ) )
+    {
+        // use potential V in hamiltonian
+        m_H = complex<double>( 0.0 , -1.0 ) * V;
+        initialize_kinetic_term();
+    }
+
+    void initialize_kinetic_term( )
+    {
+        const int N = num_rows( m_H );
+        mtl::matrix::inserter< matrix_type , mtl::update_plus< complex<double> > > ins( m_H );
+        const double z = 1.0;
+        // fill diagonal and upper and lower diagonal
+        for( int i = 0 ; i<N ; ++i )
+        {
+            ins[ i ][ (i+1) % N ] << complex< double >( 0.0 , -z );
+            ins[ i ][ i ] << complex< double >( 0.0 , z );
+            ins[ (i+1) % N ][ i ] << complex< double >( 0.0 , -z );
+        }
+    }
+
+    void operator()( const state_type &psi , state_type &dpsidt , const double t )
+    {
+        dpsidt = m_H * psi;
+    }
+
+};
+
+struct write_for_gnuplot
+{
+    size_t m_every , m_count;
+
+    write_for_gnuplot( size_t every = 10 )
+    : m_every( every ) , m_count( 0 ) { }
+
+    void operator()( const state_type &x , double t )
+    {
+        if( ( m_count % m_every ) == 0 )
+        {
+            //clog << t << endl;
+            cout << "p [0:" << mtl::size(x) << "][0:0.02] '-'" << endl;
+            for( size_t i=0 ; i<mtl::size(x) ; ++i )
+            {
+                cout << i << "\t" << norm(x[i]) << "\n";
+            }
+            cout << "e" << endl;
+        }
+
+        ++m_count;
+    }
+};
+
+static const int N = 1024;
+static const int N0 = 256;
+static const double sigma0 = 20;
+static const double k0 = -1.0;
+
+int main( int argc , char** argv )
+{
+    state_type x( N , 0.0 );
+
+    // initialize gauss packet with nonzero velocity
+    for( int i=0 ; i<N ; ++i )
+    {
+        x[i] = exp( -(i-N0)*(i-N0) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , k0*i ) );
+        //x[i] += 2.0*exp( -(i+N0-N)*(i+N0-N) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , -k0*i ) );
+    }
+    x /= mtl::two_norm( x );
+
+    typedef runge_kutta4< state_type > stepper;
+
+    // create potential barrier
+    mtl::compressed2D< double > V( N , N );
+    V = 0.0;
+    {
+        mtl::matrix::inserter< mtl::compressed2D< double > > ins( V );
+        for( int i=0 ; i<N ; ++i )
+        {
+            //ins[i][i] << 1E-4*(i-N/2)*(i-N/2);
+
+            if( i < N/2 )
+                ins[ i ][ i ] << 0.0 ;
+            else
+                ins[ i ][ i ] << 1.0 ;
+
+        }
+    }
+
+    // perform integration, output can be piped to gnuplot
+    integrate_const( stepper() , hamiltonian( V ) , x , 0.0 , 1000.0 , 0.1 , write_for_gnuplot( 10 ) );
+
+    clog << "Norm: " << mtl::two_norm( x ) << endl;
+
+    return 0;
+}

libs/numeric/odeint/examples/mtl/implicit_euler_mtl.cpp

@@ -0,0 +1,324 @@
+/*
+ * Copyright 2012 Karsten Ahnert
+ * Copyright 2012 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)
+ */
+
+
+#include <iostream>
+#include <fstream>
+#include <utility>
+#include "time.h"
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/phoenix/phoenix.hpp>
+#include <boost/numeric/mtl/mtl.hpp>
+
+#include <boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+namespace phoenix = boost::phoenix;
+
+
+
+typedef mtl::dense_vector< double > vec_mtl4;
+typedef mtl::compressed2D< double > mat_mtl4;
+
+typedef boost::numeric::ublas::vector< double > vec_ublas;
+typedef boost::numeric::ublas::matrix< double > mat_ublas;
+
+
+// two systems defined 1 & 2 both are mostly sparse with the number of element variable
+struct system1_mtl4
+{
+
+    void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
+    {
+        int size = mtl::size(x);
+
+        dxdt[ 0 ] = -0.06*x[0];
+
+        for (int i =1; i< size ; ++i){
+
+            dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
+        }
+
+    }
+};
+
+struct jacobi1_mtl4
+{
+    void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
+    {
+        int size = mtl::size(x);
+        mtl::matrix::inserter<mat_mtl4> ins(J);
+
+        ins[0][0]=-0.06;
+
+        for (int i =1; i< size ; ++i)
+        {
+            ins[i][i-1] = + 4.2;
+            ins[i][i] = -4.2*x[i];
+        }
+    }
+};
+
+
+
+struct system1_ublas
+{
+
+    void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
+    {
+        int size = x.size();
+
+        dxdt[ 0 ] = -0.06*x[0];
+
+        for (int i =1; i< size ; ++i){
+
+            dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
+        }
+
+    }
+};
+
+struct jacobi1_ublas
+{
+    void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
+    {
+        int size = x.size();
+// mtl::matrix::inserter<mat_mtl4> ins(J);
+
+        J(0,0)=-0.06;
+
+        for (int i =1; i< size ; ++i){
+//ins[i][0]=120.0*x[i];
+            J(i,i-1) = + 4.2;
+            J(i,i) = -4.2*x[i];
+
+        }
+    }
+};
+
+struct system2_mtl4
+{
+
+    void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
+    {
+        int size = mtl::size(x);
+
+
+        for (int i =0; i< size/5 ; i+=5){
+
+            dxdt[ i ] = -0.5*x[i];
+            dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
+            dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
+            dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
+            dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
+
+        }
+
+    }
+};
+
+struct jacobi2_mtl4
+{
+    void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
+    {
+        int size = mtl::size(x);
+        mtl::matrix::inserter<mat_mtl4> ins(J);
+
+        for (int i =0; i< size/5 ; i+=5){
+
+            ins[ i ][i] = -0.5;
+            ins[i+1][i+1]=25*x[i+2];
+            ins[i+1][i+2] = 25*x[i+1];
+            ins[i+1][i+3] = -740*2*x[i+3];
+            ins[i+1][i] =+4.2e-2;
+
+            ins[i+2][i]= 50*x[i];
+            ins[i+2][i+3]= -740*2*x[i+3];
+            ins[i+3][i+1] = -25*x[i+2];
+            ins[i+3][i+2] = -25*x[i+1];
+            ins[i+3][i+3] = +740*2*x[i+3];
+            ins[i+4][i] = 0.25*x[i+1];
+            ins[i+4][i+1] =0.25*x[i];
+            ins[i+4][i+3]=-44.5;
+
+
+
+        }
+    }
+};
+
+
+
+struct system2_ublas
+{
+
+    void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
+    {
+        int size = x.size();
+        for (int i =0; i< size/5 ; i+=5){
+
+            dxdt[ i ] = -4.2e-2*x[i];
+            dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
+            dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
+            dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
+            dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
+
+        }
+
+    }
+};
+
+struct jacobi2_ublas
+{
+    void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
+    {
+        int size = x.size();
+
+        for (int i =0; i< size/5 ; i+=5){
+
+            J(i ,i) = -4.2e-2;
+            J(i+1,i+1)=25*x[i+2];
+            J(i+1,i+2) = 25*x[i+1];
+            J(i+1,i+3) = -740*2*x[i+3];
+            J(i+1,i) =+4.2e-2;
+
+            J(i+2,i)= 50*x[i];
+            J(i+2,i+3)= -740*2*x[i+3];
+            J(i+3,i+1) = -25*x[i+2];
+            J(i+3,i+2) = -25*x[i+1];
+            J(i+3,i+3) = +740*2*x[i+3];
+            J(i+4,i) = 0.25*x[i+1];
+            J(i+4,i+1) =0.25*x[i];
+            J(i+4,i+3)=-44.5;
+
+
+
+        }
+
+
+    }
+};
+
+
+
+
+void testRidiculouslyMassiveArray( int size )
+{
+    typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
+    typedef boost::numeric::odeint::implicit_euler< double > booststepper;
+
+    vec_mtl4 x(size , 0.0);
+    x[0]=1;
+
+
+    double dt = 0.02;
+    double endtime = 10.0;
+
+    clock_t tstart_mtl4 = clock();
+    size_t num_of_steps_mtl4 = integrate_const(
+        mtl4stepper() ,
+        make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
+        x , 0.0 , endtime , dt  );
+    clock_t tend_mtl4 = clock() ;
+
+    clog << x[0] << endl;
+    clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
+
+    vec_ublas x_ublas(size , 0.0);
+    x_ublas[0]=1;
+
+    clock_t tstart_boost = clock();
+    size_t num_of_steps_ublas = integrate_const(
+        booststepper() ,
+        make_pair( system1_ublas() , jacobi1_ublas() ) ,
+        x_ublas , 0.0 , endtime , dt  );
+    clock_t tend_boost = clock() ;
+    
+    clog << x_ublas[0] << endl;
+    clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
+
+    clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
+    return ;
+}
+
+
+
+void testRidiculouslyMassiveArray2( int size )
+{
+    typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
+    typedef boost::numeric::odeint::implicit_euler< double > booststepper;
+
+
+    vec_mtl4 x(size , 0.0);
+    x[0]=100;
+
+
+    double dt = 0.01;
+    double endtime = 10.0;
+
+    clock_t tstart_mtl4 = clock();
+    size_t num_of_steps_mtl4 = integrate_const(
+        mtl4stepper() ,
+        make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
+        x , 0.0 , endtime , dt );
+
+
+    clock_t tend_mtl4 = clock() ;
+    
+    clog << x[0] << endl;
+    clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
+
+    vec_ublas x_ublas(size , 0.0);
+    x_ublas[0]=100;
+
+    clock_t tstart_boost = clock();
+    size_t num_of_steps_ublas = integrate_const(
+        booststepper() ,
+        make_pair( system1_ublas() , jacobi1_ublas() ) ,
+        x_ublas , 0.0 , endtime , dt  );
+
+
+    clock_t tend_boost = clock() ;
+    
+    clog << x_ublas[0] << endl;
+    clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
+
+    clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
+    return ;
+}
+
+
+
+ 
+int main( int argc , char **argv )
+{
+    std::vector< size_t > length;
+    length.push_back( 8 );
+    length.push_back( 16 );
+    length.push_back( 32 );
+    length.push_back( 64 );
+    length.push_back( 128 );
+    length.push_back( 256 );
+
+    for( size_t i=0 ; i<length.size() ; ++i )
+    {
+        clog << "Testing with size " << length[i] << endl;
+        testRidiculouslyMassiveArray( length[i] );
+    }
+    clog << endl << endl;
+
+    for( size_t i=0 ; i<length.size() ; ++i )
+    {
+        clog << "Testing with size " << length[i] << endl;
+        testRidiculouslyMassiveArray2( length[i] );
+    }
+}

libs/numeric/odeint/examples/multiprecision/Jamfile

@@ -0,0 +1,16 @@
+# Copyright 2011-2013 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# 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)
+
+
+project
+    : requirements 
+      <define>BOOST_ALL_NO_LIB=1
+    : 
+    ;
+    
+
+exe lorenz_mp : lorenz_mp.cpp ;
+exe cmp_precision : cmp_precision.cpp ;

libs/numeric/odeint/examples/multiprecision/cmp_precision.cpp

@@ -0,0 +1,68 @@
+/* Boost libs/numeric/odeint/examples/multiprecision/cmp_precision.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ example comparing double to multiprecision using Boost.Multiprecision
+
+ 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)
+ */
+
+
+#include <iostream>
+#include <boost/numeric/odeint.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::multiprecision::cpp_dec_float_50 mp_50;
+
+/* we solve the simple ODE x' = 3/(2t^2) + x/(2t)
+ * with initial condition x(1) = 0.
+ * Analytic solution is x(t) = sqrt(t) - 1/t
+ */
+
+void rhs_m( const mp_50 x , mp_50 &dxdt , const mp_50 t )
+{   // version for multiprecision
+    dxdt = mp_50(3)/(mp_50(2)*t*t) + x/(mp_50(2)*t);
+}
+
+void rhs_d( const double x , double &dxdt , const double t )
+{   // version for double precision
+    dxdt = 3.0/(2.0*t*t) + x/(2.0*t);
+}
+
+// state_type = mp_50 = deriv_type = time_type = mp_50
+typedef runge_kutta4< mp_50 , mp_50 , mp_50 , mp_50 , vector_space_algebra , default_operations , never_resizer > stepper_type_m;
+
+typedef runge_kutta4< double , double , double , double , vector_space_algebra , default_operations , never_resizer > stepper_type_d;
+
+int main()
+{
+
+    stepper_type_m stepper_m;
+    stepper_type_d stepper_d;
+
+    mp_50 dt_m( 0.5 );
+    double dt_d( 0.5 );
+
+    cout << "dt" << '\t' << "mp" << '\t' << "double" << endl;
+    
+    while( dt_m > 1E-20 )
+    {
+
+        mp_50 x_m = 0; //initial value x(1) = 0
+        stepper_m.do_step( rhs_m , x_m , mp_50( 1 ) , dt_m );
+        double x_d = 0;
+        stepper_d.do_step( rhs_d , x_d , 1.0 , dt_d );        
+
+        cout << dt_m << '\t';
+        cout << abs((x_m - (sqrt(1+dt_m)-mp_50(1)/(1+dt_m)))/x_m) << '\t' ;
+        cout << abs((x_d - (sqrt(1+dt_d)-mp_50(1)/(1+dt_d)))/x_d) << endl ;
+        dt_m /= 2;
+        dt_d /= 2;
+    }
+}

libs/numeric/odeint/examples/multiprecision/lorenz_mp.cpp

@@ -0,0 +1,81 @@
+/*
+ * lorenz_mp.cpp
+ *
+ * This example demonstrates how odeint can be used with boost.multiprecision.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2013 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)
+ */
+
+
+
+#include <iostream>
+
+//[ mp_lorenz_defs
+#include <boost/numeric/odeint.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::multiprecision::cpp_dec_float_50 value_type;
+
+typedef boost::array< value_type , 3 > state_type;
+//]
+
+//[ mp_lorenz_rhs
+struct lorenz
+{
+    void operator()( const state_type &x , state_type &dxdt , value_type t ) const
+    {
+        const value_type sigma( 10 );
+        const value_type R( 28 );
+        const value_type b( value_type( 8 ) / value_type( 3 ) );
+
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+//]
+
+
+
+
+struct streaming_observer
+{
+    std::ostream& m_out;
+
+    streaming_observer( std::ostream &out ) : m_out( out ) { }
+
+    template< class State , class Time >
+    void operator()( const State &x , Time t ) const
+    {
+        m_out << t;
+        for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i] ;
+        m_out << "\n";
+    }
+};
+
+
+
+
+
+
+int main( int argc , char **argv )
+{
+    //[ mp_lorenz_int
+    state_type x = {{ value_type( 10.0 ) , value_type( 10.0 ) , value_type( 10.0 ) }};
+
+    cout.precision( 50 );
+    integrate_const( runge_kutta4< state_type , value_type >() ,
+            lorenz() , x , value_type( 0.0 ) , value_type( 10.0 ) , value_type( value_type( 1.0 ) / value_type( 10.0 ) ) ,
+            streaming_observer( cout ) );
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/my_vector.cpp

@@ -0,0 +1,114 @@
+/*
+ * Copyright 2011-2012 Mario Mulansky
+ * Copyright 2012-2013 Karsten Ahnert
+ *
+ * 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)
+ *
+ * Example for self defined vector type.
+ */
+
+#include <vector>
+
+#include <boost/numeric/odeint.hpp>
+
+//[my_vector
+template< size_t MAX_N >
+class my_vector
+{
+    typedef std::vector< double > vector;
+
+public:
+    typedef vector::iterator iterator;
+    typedef vector::const_iterator const_iterator;
+
+public:
+    my_vector( const size_t N )
+        : m_v( N )
+    { 
+        m_v.reserve( MAX_N );
+    }
+
+    my_vector()
+        : m_v()
+    {
+        m_v.reserve( MAX_N );
+    }
+
+// ... [ implement container interface ]
+//]
+    const double & operator[]( const size_t n ) const
+    { return m_v[n]; }
+
+    double & operator[]( const size_t n )
+    { return m_v[n]; }
+
+    iterator begin()
+    { return m_v.begin(); }
+
+    const_iterator begin() const
+    { return m_v.begin(); }
+
+    iterator end()
+    { return m_v.end(); }
+    
+    const_iterator end() const
+    { return m_v.end(); }
+
+    size_t size() const 
+    { return m_v.size(); }
+
+    void resize( const size_t n )
+    { m_v.resize( n ); }
+
+private:
+    std::vector< double > m_v;
+
+};
+
+//[my_vector_resizeable
+// define my_vector as resizeable
+
+namespace boost { namespace numeric { namespace odeint {
+
+template<size_t N>
+struct is_resizeable< my_vector<N> >
+{
+    typedef boost::true_type type;
+    static const bool value = type::value;
+};
+
+} } }
+//]
+
+
+typedef my_vector<3> state_type;
+
+void lorenz( const state_type &x , state_type &dxdt , const 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] = -b * x[2] + x[0] * x[1];
+}
+
+using namespace boost::numeric::odeint;
+
+int main()
+{
+    state_type x(3);
+    x[0] = 5.0 ; x[1] = 10.0 ; x[2] = 10.0;
+
+    // make sure resizing is ON
+    BOOST_STATIC_ASSERT( is_resizeable<state_type>::value == true );
+
+    // my_vector works with range_algebra as it implements
+    // the required parts of a container interface
+    // no further work is required
+
+    integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 10.0 , 0.1 );
+}

libs/numeric/odeint/examples/nt2/Jamfile.v2

@@ -0,0 +1,34 @@
+#==============================================================================
+#         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+#         Copyright 2014          NumScale SAS
+#
+#          Distributed under the Boost Software License, Version 1.0.
+#                 See accompanying file LICENSE.txt or copy at
+#                     http://www.boost.org/LICENSE_1_0.txt
+#==============================================================================
+
+use-project boost : $(BOOST_ROOT) ;
+
+import os ;
+
+# This must be built using an NT2 installation.
+# NT2_ROOT_PATH should point to the build directory.
+# Currently, cxxflags needs to be set to the required architecture
+# if using avx/avx2, set the environemnt variable NT2_SIMD_FLAGS to the
+# required value for your compiler (i.e. -mavx2 on g++)
+# If using sse2/3/4 in 64 bits, this is set automatically.
+
+local NT2_ROOT_PATH = [ os.environ NT2_ROOT_PATH ] ;
+local NT2_SIMD_FLAGS = [ os.environ NT2_SIMD_FLAGS ] ;
+
+project
+  : requirements
+    <define>BOOST_ALL_NO_LIB=1
+    <include>$(NT2_ROOT_PATH)/include/
+    <link>static
+    <toolset>gcc:<cxxflags>-DBOOST_SIMD_NO_STRICT_ALIASING
+    <toolset>gcc:<cxxflags>-fno-strict-aliasing
+    <cxxflags>$(NT2_SIMD_FLAGS)
+  ;
+
+exe phase_oscillator_ensemble : phase_oscillator_ensemble.cpp ;

libs/numeric/odeint/examples/nt2/phase_oscillator_ensemble.cpp

@@ -0,0 +1,170 @@
+//==============================================================================
+//         Copyright 2011-2014     Karsten Ahnert
+//         Copyright 2011-2014     Mario Mulansky
+//         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+//         Copyright 2014          NumScale SAS
+//
+//          Distributed under the Boost Software License, Version 1.0.
+//                 See accompanying file LICENSE.txt or copy at
+//                     http://www.boost.org/LICENSE_1_0.txt
+//==============================================================================
+
+#include <iostream>
+#include <utility>
+
+#include <boost/numeric/odeint.hpp>
+
+#ifndef M_PI //not there on windows
+#define M_PI 3.141592653589793 //...
+#endif
+
+#include <boost/random.hpp>
+#include <boost/dispatch/meta/as_integer.hpp>
+
+#include <nt2/include/functions/cos.hpp>
+#include <nt2/include/functions/sin.hpp>
+#include <nt2/include/functions/atan2.hpp>
+#include <nt2/table.hpp>
+#include <nt2/include/functions/zeros.hpp>
+#include <nt2/include/functions/sum.hpp>
+#include <nt2/include/functions/mean.hpp>
+#include <nt2/arithmetic/include/functions/hypot.hpp>
+#include <nt2/include/functions/tie.hpp>
+
+#include <boost/numeric/odeint/external/nt2/nt2_algebra_dispatcher.hpp>
+
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+template <typename container_type, typename T>
+pair< T, T > calc_mean_field( const container_type &x )
+
+{
+  T cos_sum = 0.0 , sin_sum = 0.0;
+
+  nt2::tie(cos_sum,sin_sum) = nt2::tie(nt2::mean( nt2::cos(x) ), nt2::mean( nt2::sin(x) ));
+
+  T K = nt2::hypot(sin_sum,cos_sum);
+  T Theta = nt2::atan2( sin_sum , cos_sum );
+
+  return make_pair( K , Theta );
+}
+
+template <typename container_type, typename T>
+struct phase_ensemble
+{
+  typedef typename boost::dispatch::meta::as_integer<T,unsigned>::type int_type;
+  container_type m_omega;
+  T m_epsilon;
+
+  phase_ensemble( const int_type n , T g = 1.0 , T epsilon = 1.0 )
+  : m_epsilon( epsilon )
+  {
+    m_omega = nt2::zeros(nt2::of_size(n), nt2::meta::as_<T>());
+    create_frequencies( g );
+  }
+
+  void create_frequencies( T g )
+  {
+    boost::mt19937 rng;
+    boost::cauchy_distribution<> cauchy( 0.0 , g );
+    boost::variate_generator< boost::mt19937&, boost::cauchy_distribution<> > gen( rng , cauchy );
+    generate( m_omega.begin() , m_omega.end() , gen );
+}
+
+  void set_epsilon( T epsilon ) { m_epsilon = epsilon; }
+
+  T get_epsilon( void ) const { return m_epsilon; }
+
+  void operator()( const container_type &x , container_type &dxdt , T ) const
+  {
+    pair< T, T > mean = calc_mean_field<container_type,T>( x );
+    dxdt = m_omega +  m_epsilon * mean.first * nt2::sin( mean.second - x );
+  }
+};
+
+template<typename T>
+struct statistics_observer
+{
+    typedef typename boost::dispatch::meta::as_integer<T,unsigned>::type int_type;
+    T m_K_mean;
+    int_type m_count;
+
+    statistics_observer( void )
+    : m_K_mean( 0.0 ) , m_count( 0 ) { }
+
+    template< class State >
+    void operator()( const State &x , T t )
+    {
+        pair< T, T > mean = calc_mean_field<State,T>( x );
+        m_K_mean += mean.first;
+        ++m_count;
+    }
+
+    T get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / T( m_count ) : 0.0 ; }
+
+    void reset( void ) { m_K_mean = 0.0; m_count = 0; }
+};
+
+template<typename T>
+struct test_ode_table
+{
+  typedef nt2::table<T> array_type;
+  typedef void experiment_is_immutable;
+
+  typedef typename boost::dispatch::meta::as_integer<T,unsigned>::type int_type;
+
+  test_ode_table ( )
+                 : size_(16384), ensemble( size_ , 1.0 ), unif( 0.0 , 2.0 * M_PI ), gen( rng , unif ), obs()
+  {
+    x.resize(nt2::of_size(size_));
+  }
+
+  void operator()()
+  {
+    for( T epsilon = 0.0 ; epsilon < 5.0 ; epsilon += 0.1 )
+    {
+      ensemble.set_epsilon( epsilon );
+      obs.reset();
+
+      // start with random initial conditions
+      generate( x.begin() , x.end() , gen );
+      // calculate some transients steps
+      integrate_const( runge_kutta4< array_type, T >() , boost::ref( ensemble ) , x , T(0.0) , T(10.0) , dt );
+
+      // integrate and compute the statistics
+      integrate_const( runge_kutta4< array_type, T >() , boost::ref( ensemble ) , x , T(0.0) , T(100.0) , dt , boost::ref( obs ) );
+      cout << epsilon << "\t" << obs.get_K_mean() << endl;
+    }
+  }
+
+  friend std::ostream& operator<<(std::ostream& os, test_ode_table<T> const& p)
+  {
+    return os << "(" << p.size() << ")";
+  }
+
+  std::size_t size() const { return size_; }
+
+  private:
+    std::size_t size_;
+    phase_ensemble<array_type,T> ensemble;
+    boost::uniform_real<> unif;
+    array_type x;
+    boost::mt19937 rng;
+    boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen;
+    statistics_observer<T> obs;
+
+    static const T dt = 0.1;
+};
+
+int main()
+{
+  std::cout<< " With T = [double] \n";
+  test_ode_table<double> test_double;
+  test_double();
+
+  std::cout<< " With T = [float] \n";
+  test_ode_table<float> test_float;
+  test_float();
+}

libs/numeric/odeint/examples/openmp/Jamfile.v2

@@ -0,0 +1,23 @@
+# Copyright 2011-2013 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# Copyright 2013 Pascal Germroth
+# 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)
+
+use-project /boost : $(BOOST_ROOT) ;
+import openmp : * ;
+
+project
+    : requirements
+      <include>..
+      <define>BOOST_ALL_NO_LIB=1
+      <library>/boost//timer
+      [ openmp ]
+    ;
+
+exe lorenz_ensemble : lorenz_ensemble.cpp ;
+exe lorenz_ensemble_simple : lorenz_ensemble_simple.cpp ;
+exe lorenz_ensemble_nested : lorenz_ensemble_nested.cpp ;
+exe phase_chain : phase_chain.cpp ;
+exe phase_chain_omp_state : phase_chain_omp_state.cpp ;

libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp

@@ -0,0 +1,91 @@
+/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Parallelized Lorenz ensembles
+
+ 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)
+ */
+
+#include <omp.h>
+#include <vector>
+#include <iostream>
+#include <iterator>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+#include <boost/lexical_cast.hpp>
+#include "point_type.hpp"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef point<double, 3> point_type;
+typedef vector< point_type > inner_state_type;
+typedef openmp_state<point_type> state_type;
+
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+
+struct sys_func {
+    const vector<double> &R;
+    sys_func( vector<double> &R ) : R(R) {}
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const {
+#       pragma omp parallel for
+        for(size_t j = 0 ; j < x.size() ; j++) {
+            size_t offset = 0;
+            for(size_t i = 0 ; i < j ; i++)
+                offset += x[i].size();
+
+            for(size_t i = 0 ; i < x[j].size() ; i++) {
+                const point_type &xi = x[j][i];
+                point_type &dxdti = dxdt[j][i];
+                dxdti[0] = -sigma * (xi[0] - xi[1]);
+                dxdti[1] = R[offset + i] * xi[0] - xi[1] - xi[0] * xi[2];
+                dxdti[2] = -b * xi[2] + xi[0] * xi[1];
+            }
+        }
+    }
+};
+
+
+int main(int argc, char **argv) {
+    size_t n = 1024;
+    if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]);
+
+    vector<double> R(n);
+    const double Rmin = 0.1, Rmax = 50.0;
+#   pragma omp parallel for
+    for(size_t i = 0 ; i < n ; i++)
+        R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i;
+
+    vector<point_type> inner(n, point_type(10, 10, 10));
+    state_type state;
+    split(inner, state);
+
+    cerr << "openmp_state split " << n << " into";
+    for(size_t i = 0 ; i != state.size() ; i++)
+        cerr << ' ' << state[i].size();
+    cerr << endl;
+
+    typedef runge_kutta4< state_type, double > stepper;
+
+    const double t_max = 10.0, dt = 0.01;
+
+    integrate_const(
+        stepper(),
+        sys_func(R),
+        state,
+        0.0, t_max, dt
+    );
+
+    unsplit(state, inner);
+    std::copy( inner.begin(), inner.end(), ostream_iterator<point_type>(cout, "\n") );
+
+    return 0;
+}

libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp

@@ -0,0 +1,75 @@
+/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Pascal Germroth
+ Copyright 2013 Mario Mulansky
+
+ Parallelized Lorenz ensembles using nested omp algebra
+
+ 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)
+ */
+
+#include <omp.h>
+#include <vector>
+#include <iostream>
+#include <iterator>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+#include <boost/lexical_cast.hpp>
+#include "point_type.hpp"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef point<double, 3> point_type;
+typedef vector< point_type > state_type;
+
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+
+struct sys_func {
+    const vector<double> &R;
+    sys_func( vector<double> &R ) : R(R) {}
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const {
+#       pragma omp parallel for
+        for(size_t i = 0 ; i < x.size() ; i++) {
+            dxdt[i][0] = -sigma * (x[i][0] - x[i][1]);
+            dxdt[i][1] = R[i] * x[i][0] - x[i][1] - x[i][0] * x[i][2];
+            dxdt[i][2] = -b * x[i][2] + x[i][0] * x[i][1];
+        }
+    }
+};
+
+
+int main(int argc, char **argv) {
+    size_t n = 1024;
+    if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]);
+
+    vector<double> R(n);
+    const double Rmin = 0.1, Rmax = 50.0;
+#   pragma omp parallel for
+    for(size_t i = 0 ; i < n ; i++)
+        R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i;
+
+    state_type state( n , point_type(10, 10, 10) );
+
+    typedef runge_kutta4< state_type, double , state_type , double ,
+                          openmp_nested_algebra<vector_space_algebra> > stepper;
+
+    const double t_max = 10.0, dt = 0.01;
+
+    integrate_const(
+        stepper(),
+        sys_func(R),
+        state,
+        0.0, t_max, dt
+    );
+
+    std::copy( state.begin(), state.end(), ostream_iterator<point_type>(cout, "\n") );
+
+    return 0;
+}

libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp

@@ -0,0 +1,79 @@
+/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Parallelized Lorenz ensembles
+
+ 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)
+ */
+
+#include <omp.h>
+#include <vector>
+#include <iostream>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+#include "point_type.hpp"
+
+using namespace std;
+
+typedef vector<double> vector_type;
+typedef point<double, 3> point_type;
+typedef vector<point_type> state_type;
+
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+struct sys_func {
+    const vector_type &R;
+    sys_func( const vector_type &_R ) : R( _R ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const {
+        const size_t n = x.size();
+#       pragma omp parallel for
+        for(size_t i = 0 ; i < n ; i++) {
+            const point_type &xi = x[i];
+            point_type &dxdti = dxdt[i];
+            dxdti[0] = -sigma * (xi[0] - xi[1]);
+            dxdti[1] = R[i] * xi[0] - xi[1] - xi[0] * xi[2];
+            dxdti[2] = -b * xi[2] + xi[0] * xi[1];
+        }
+    }
+};
+
+
+int main() {
+    using namespace boost::numeric::odeint;
+
+    const size_t n = 1024;
+    vector_type R(n);
+    const double Rmin = 0.1, Rmax = 50.0;
+#   pragma omp parallel for
+    for(size_t i = 0 ; i < n ; i++)
+        R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i;
+
+    state_type X(n, point_type(10, 10, 10));
+
+    const double t_max = 10.0, dt = 0.01;
+
+    // Simple stepper with constant step size
+    // typedef runge_kutta4<state_type, double, state_type, double,
+    //                      openmp_range_algebra> stepper;
+
+    // integrate_const(stepper(), sys_func(R), X, 0.0, t_max, dt);
+
+    // Controlled stepper with variable step size
+    typedef runge_kutta_fehlberg78<state_type, double, state_type, double,
+                                   openmp_range_algebra> error_stepper_type;
+    typedef controlled_runge_kutta<error_stepper_type> controlled_stepper_type;
+    controlled_stepper_type controlled_stepper;
+
+    integrate_adaptive(controlled_stepper, sys_func(R), X, 0.0, t_max, dt);
+
+    copy( X.begin(), X.end(), ostream_iterator<point_type>(cout, "\n") );
+
+    return 0;
+}

libs/numeric/odeint/examples/openmp/openmp.jam

@@ -0,0 +1,52 @@
+# Copyright 2013 Karsten Ahnert
+# Copyright 2013 Mario Mulansky
+# Copyright 2013 Pascal Germroth
+# 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)
+
+
+# Only builds target with supported OpenMP enabled toolsets.
+#
+# use as:
+#   exe omp : omp.cpp : [ openmp ] ;
+#
+rule openmp return
+    # default
+    <build>no
+    # GNU C++
+    <toolset>gcc:<cxxflags>-fopenmp
+    <toolset>gcc:<linkflags>-fopenmp
+    <toolset>gcc:<build>yes
+    # Microsoft Visual C++
+    <toolset>msvc:<cxxflags>/openmp
+    <toolset>msvc:<linkflags>/openmp
+    <toolset>msvc:<build>yes
+    # Intel C++
+    <toolset>intel-linux:<cxxflags>-openmp
+    <toolset>intel-linux:<linkflags>-openmp
+    <toolset>intel-linux:<build>yes
+    <toolset>intel-win:<cxxflags>-Qopenmp
+    <toolset>intel-win:<linkflags>-Qopenmp
+    <toolset>intel-win:<build>yes
+    # HP aC++
+    <toolset>acc:<cxxflags>+Oopenmp
+    <toolset>acc:<linkflags>+Oopenmp
+    <toolset>acc:<build>yes
+    # Sun Studio
+    <toolset>sun:<cxxflags>-xopenmp
+    <toolset>sun:<linkflags>-xopenmp
+    <toolset>sun:<build>yes
+    # IBM XL
+    <toolset>vacpp:<cxxflags>-qsmp=omp
+    <toolset>vacpp:<linkflags>-qsmp=omp
+    <toolset>vacpp:<build>yes
+    # PG++
+    <toolset>pgi:<cxxflags>-mp
+    <toolset>pgi:<linkflags>-mp
+    <toolset>pgi:<build>yes
+    # Pathscale
+    <toolset>pathscale:<cxxflags>-mp
+    <toolset>pathscale:<linkflags>-mp
+    <toolset>pathscale:<build>yes
+    ;

libs/numeric/odeint/examples/openmp/phase_chain.cpp

@@ -0,0 +1,95 @@
+/*
+ * phase_chain.cpp
+ *
+ * Example of OMP parallelization with odeint
+ *
+ * Copyright 2013 Karsten Ahnert
+ * Copyright 2013 Mario Mulansky
+ * Copyright 2013 Pascal Germroth
+ * 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)
+ */
+
+#include <iostream>
+#include <vector>
+#include <boost/random.hpp>
+#include <boost/timer/timer.hpp>
+//[phase_chain_openmp_header
+#include <omp.h>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+//]
+
+using namespace std;
+using namespace boost::numeric::odeint;
+using boost::timer::cpu_timer;
+using boost::math::double_constants::pi;
+
+//[phase_chain_vector_state
+typedef std::vector< double > state_type;
+//]
+
+//[phase_chain_rhs
+struct phase_chain
+{
+    phase_chain( double gamma = 0.5 )
+    : m_gamma( gamma ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
+    {
+        const size_t N = x.size();
+        #pragma omp parallel for schedule(runtime)
+        for(size_t i = 1 ; i < N - 1 ; ++i)
+        {
+            dxdt[i] = coupling_func( x[i+1] - x[i] ) +
+                      coupling_func( x[i-1] - x[i] );
+        }
+        dxdt[0  ] = coupling_func( x[1  ] - x[0  ] );
+        dxdt[N-1] = coupling_func( x[N-2] - x[N-1] );
+    }
+
+    double coupling_func( double x ) const
+    {
+        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
+    }
+
+    double m_gamma;
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+    //[phase_chain_init
+    size_t N = 131101;
+    state_type x( N );
+    boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
+    boost::random::mt19937 engine( 0 );
+    generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
+    //]
+
+    //[phase_chain_stepper
+    typedef runge_kutta4<
+                      state_type , double ,
+                      state_type , double ,
+                      openmp_range_algebra
+                    > stepper_type;
+    //]
+
+    //[phase_chain_scheduling
+    int chunk_size = N/omp_get_max_threads();
+    omp_set_schedule( omp_sched_static , chunk_size );
+    //]
+
+    cpu_timer timer;
+    //[phase_chain_integrate
+    integrate_n_steps( stepper_type() , phase_chain( 1.2 ) ,
+                       x , 0.0 , 0.01 , 100 );
+    //]
+    double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
+    std::cerr << run_time << "s" << std::endl;
+    // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
+
+    return 0;
+}

libs/numeric/odeint/examples/openmp/phase_chain_omp_state.cpp

@@ -0,0 +1,98 @@
+/*
+ * phase_chain_omp_state.cpp
+ *
+ * Example of OMP parallelization with odeint
+ *
+ * Copyright 2013 Karsten Ahnert
+ * Copyright 2013 Mario Mulansky
+ * Copyright 2013 Pascal Germroth
+ * 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)
+ */
+
+#include <iostream>
+#include <vector>
+#include <boost/random.hpp>
+#include <boost/timer/timer.hpp>
+#include <omp.h>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+using boost::timer::cpu_timer;
+using boost::math::double_constants::pi;
+
+typedef openmp_state<double> state_type;
+
+//[phase_chain_state_rhs
+struct phase_chain_omp_state
+{
+    phase_chain_omp_state( double gamma = 0.5 )
+    : m_gamma( gamma ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
+    {
+        const size_t N = x.size();
+        #pragma omp parallel for schedule(runtime)
+        for(size_t n = 0 ; n < N ; ++n)
+        {
+            const size_t M = x[n].size();
+            for(size_t m = 1 ; m < M-1 ; ++m)
+            {
+                dxdt[n][m] = coupling_func( x[n][m+1] - x[n][m] ) +
+                             coupling_func( x[n][m-1] - x[n][m] );
+            }
+            dxdt[n][0] = coupling_func( x[n][1] - x[n][0] );
+            if( n > 0 )
+            {
+                dxdt[n][0] += coupling_func( x[n-1].back() - x[n].front() );
+            }
+            dxdt[n][M-1] = coupling_func( x[n][M-2] - x[n][M-1] );
+            if( n < N-1 )
+            {
+                dxdt[n][M-1] += coupling_func( x[n+1].front() - x[n].back() );
+            }
+        }
+    }
+
+    double coupling_func( double x ) const
+    {
+        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
+    }
+
+    double m_gamma;
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+    //[phase_chain_state_init
+    const size_t N = 131101;
+    vector<double> x( N );
+    boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
+    boost::random::mt19937 engine( 0 );
+    generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
+    const size_t blocks = omp_get_max_threads();
+    state_type x_split( blocks );
+    split( x , x_split );
+    //]
+
+
+    cpu_timer timer;
+    //[phase_chain_state_integrate
+    integrate_n_steps( runge_kutta4<state_type>() , phase_chain_omp_state( 1.2 ) ,
+                       x_split , 0.0 , 0.01 , 100 );
+    unsplit( x_split , x );
+    //]
+
+    double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
+    std::cerr << run_time << "s" << std::endl;
+    // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
+
+    return 0;
+}

libs/numeric/odeint/examples/phase_oscillator_ensemble.cpp

@@ -0,0 +1,151 @@
+/*
+ * phase_oscillator_ensemble.cpp
+ *
+ * Demonstrates the phase transition from an unsynchronized to an synchronized state.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2011-2012 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)
+ *
+ */
+
+#include <iostream>
+#include <utility>
+
+#include <boost/numeric/odeint.hpp>
+
+#ifndef M_PI //not there on windows
+#define M_PI 3.141592653589793 //...
+#endif
+
+#include <boost/random.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+//[ phase_oscillator_ensemble_system_function
+typedef vector< double > container_type;
+
+
+pair< double , double > calc_mean_field( const container_type &x )
+{
+    size_t n = x.size();
+    double cos_sum = 0.0 , sin_sum = 0.0;
+    for( size_t i=0 ; i<n ; ++i )
+    {
+        cos_sum += cos( x[i] );
+        sin_sum += sin( x[i] );
+    }
+    cos_sum /= double( n );
+    sin_sum /= double( n );
+
+    double K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum );
+    double Theta = atan2( sin_sum , cos_sum );
+
+    return make_pair( K , Theta );
+}
+
+
+struct phase_ensemble
+{
+    container_type m_omega;
+    double m_epsilon;
+
+    phase_ensemble( const size_t n , double g = 1.0 , double epsilon = 1.0 )
+    : m_omega( n , 0.0 ) , m_epsilon( epsilon )
+    {
+        create_frequencies( g );
+    }
+
+    void create_frequencies( double g )
+    {
+        boost::mt19937 rng;
+        boost::cauchy_distribution<> cauchy( 0.0 , g );
+        boost::variate_generator< boost::mt19937&, boost::cauchy_distribution<> > gen( rng , cauchy );
+        generate( m_omega.begin() , m_omega.end() , gen );
+    }
+
+    void set_epsilon( double epsilon ) { m_epsilon = epsilon; }
+
+    double get_epsilon( void ) const { return m_epsilon; }
+
+    void operator()( const container_type &x , container_type &dxdt , double /* t */ ) const
+    {
+        pair< double , double > mean = calc_mean_field( x );
+        for( size_t i=0 ; i<x.size() ; ++i )
+            dxdt[i] = m_omega[i] + m_epsilon * mean.first * sin( mean.second - x[i] );
+    }
+};
+//]
+
+
+
+//[ phase_oscillator_ensemble_observer
+struct statistics_observer
+{
+    double m_K_mean;
+    size_t m_count;
+
+    statistics_observer( void )
+    : m_K_mean( 0.0 ) , m_count( 0 ) { }
+
+    template< class State >
+    void operator()( const State &x , double t )
+    {
+        pair< double , double > mean = calc_mean_field( x );
+        m_K_mean += mean.first;
+        ++m_count;
+    }
+
+    double get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / double( m_count ) : 0.0 ; }
+
+    void reset( void ) { m_K_mean = 0.0; m_count = 0; }
+};
+//]
+
+
+
+
+
+
+
+
+int main( int argc , char **argv )
+{
+    //[ phase_oscillator_ensemble_integration
+    const size_t n = 16384;
+    const double dt = 0.1;
+
+    container_type x( n );
+
+    boost::mt19937 rng;
+    boost::uniform_real<> unif( 0.0 , 2.0 * M_PI );
+    boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
+
+    // gamma = 1, the phase transition occurs at epsilon = 2
+    phase_ensemble ensemble( n , 1.0 );
+    statistics_observer obs;
+
+    for( double epsilon = 0.0 ; epsilon < 5.0 ; epsilon += 0.1 )
+    {
+        ensemble.set_epsilon( epsilon );
+        obs.reset();
+
+        // start with random initial conditions
+        generate( x.begin() , x.end() , gen );
+
+        // calculate some transients steps
+        integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 10.0 , dt );
+
+        // integrate and compute the statistics
+        integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 100.0 , dt , boost::ref( obs ) );
+        cout << epsilon << "\t" << obs.get_K_mean() << endl;
+    }
+
+
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/point_type.hpp

@@ -0,0 +1,177 @@
+/* Boost libs/numeric/odeint/examples/point_type.hpp
+
+ Copyright 2010-2012 Karsten Ahnert
+ Copyright 2011 Mario Mulansky
+
+ solar system example for Hamiltonian stepper
+
+ 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 POINT_TYPE_HPP_INCLUDED
+#define POINT_TYPE_HPP_INCLUDED
+
+
+#include <boost/operators.hpp>
+#include <ostream>
+
+
+//[ point_type
+/*the point type */
+template< class T , size_t Dim >
+class point :
+    boost::additive1< point< T , Dim > ,
+    boost::additive2< point< T , Dim  > , T ,
+    boost::multiplicative2< point< T , Dim > , T
+    > > >
+    {
+    public:
+
+        const static size_t dim = Dim;
+        typedef T value_type;
+        typedef point< value_type , dim > point_type;
+
+        // ...
+        // constructors
+        //<-
+        point( void )
+        {
+            for( size_t i=0 ; i<dim ; ++i ) m_val[i] = 0.0;
+        }
+
+        point( value_type val )
+        {
+            for( size_t i=0 ; i<dim ; ++i ) m_val[i] = val;
+        }
+
+        point( value_type x , value_type y , value_type z = 0.0 )
+        {
+            if( dim > 0 ) m_val[0] = x;
+            if( dim > 1 ) m_val[1] = y;
+            if( dim > 2 ) m_val[2] = z;
+        }
+        //->
+
+        // ...
+        // operators
+        //<-
+        T operator[]( size_t i ) const { return m_val[i]; }
+        T& operator[]( size_t i ) { return m_val[i]; }
+
+        point_type& operator+=( const point_type& p )
+        {
+            for( size_t i=0 ; i<dim ; ++i )
+                m_val[i] += p[i];
+            return *this;
+        }
+
+        point_type& operator-=( const point_type& p )
+        {
+            for( size_t i=0 ; i<dim ; ++i )
+                m_val[i] -= p[i];
+            return *this;
+        }
+
+        point_type& operator+=( const value_type& val )
+        {
+            for( size_t i=0 ; i<dim ; ++i )
+                m_val[i] += val;
+            return *this;
+        }
+
+        point_type& operator-=( const value_type& val )
+        {
+            for( size_t i=0 ; i<dim ; ++i )
+                m_val[i] -= val;
+            return *this;
+        }
+
+        point_type& operator*=( const value_type &val )
+        {
+            for( size_t i=0 ; i<dim ; ++i )
+                m_val[i] *= val;
+            return *this;
+        }
+
+        point_type& operator/=( const value_type &val )
+        {
+            for( size_t i=0 ; i<dim ; ++i )
+                m_val[i] /= val;
+            return *this;
+        }
+
+        //->
+
+    private:
+
+        T m_val[dim];
+    };
+
+    //...
+    // more operators
+    //]
+
+    //
+    // the - operator
+    //
+    template< class T , size_t Dim >
+    point< T , Dim > operator-( const point< T , Dim > &p )
+    {
+        point< T , Dim > tmp;
+        for( size_t i=0 ; i<Dim ; ++i ) tmp[i] = -p[i];
+        return tmp;
+    }
+
+    //
+    // scalar product
+    //
+    template< class T , size_t Dim >
+    T scalar_prod( const point< T , Dim > &p1 , const point< T , Dim > &p2 )
+    {
+        T tmp = 0.0;
+        for( size_t i=0 ; i<Dim ; ++i ) tmp += p1[i] * p2[i];
+        return tmp;
+    }
+
+
+
+    //
+    // norm
+    //
+    template< class T , size_t Dim >
+    T norm( const point< T , Dim > &p1 )
+    {
+        return scalar_prod( p1 , p1 );
+    }
+
+
+
+
+    //
+    // absolute value
+    //
+    template< class T , size_t Dim >
+    T abs( const point< T , Dim > &p1 )
+    {
+        return sqrt( norm( p1 ) );
+    }
+
+
+
+
+    //
+    // output operator
+    //
+    template< class T , size_t Dim >
+    std::ostream& operator<<( std::ostream &out , const point< T , Dim > &p )
+    {
+        if( Dim > 0 ) out << p[0];
+        for( size_t i=1 ; i<Dim ; ++i ) out << " " << p[i];
+        return out;
+    }
+
+
+
+#endif //POINT_TYPE_HPP_INCLUDED

libs/numeric/odeint/examples/quadmath/Jamfile.v2

@@ -0,0 +1,16 @@
+# Copyright 2011 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# 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)
+
+
+project
+    : requirements 
+      <define>BOOST_ALL_NO_LIB=1
+    : 
+    ;
+    
+lib quadmath : : <name>quadmath <link>shared ;
+
+exe black_hole : black_hole.cpp quadmath : <cxxflags>-std=c++0x ;

libs/numeric/odeint/examples/quadmath/black_hole.cpp

@@ -0,0 +1,151 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/examples/black_hole.cpp
+
+ [begin_description]
+ This example shows how the __float128 from gcc libquadmath can be used with odeint.
+ [end_description]
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 Lee Hodgkinson
+ Copyright 2012 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)
+ */
+
+#include <cstdlib>
+#include <cmath>
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <algorithm>
+#include <cassert>
+#include <vector>
+#include <complex>
+
+extern "C" {
+#include <quadmath.h>
+}
+
+const __float128 zero =strtoflt128 ("0.0", NULL);
+
+namespace std {
+
+    inline __float128 abs( __float128 x )
+    {
+        return fabsq( x );
+    }
+
+    inline __float128 sqrt( __float128 x )
+    {
+        return sqrtq( x );
+    }
+
+    inline __float128 pow( __float128 x , __float128 y )
+    {
+        return powq( x , y );
+    }
+
+    inline __float128 abs( std::complex< __float128 > x )
+    {
+        return sqrtq( x.real() * x.real() + x.imag() * x.imag() );
+    }
+
+    inline std::complex< __float128 > pow( std::complex< __float128> x , __float128 y )
+    {
+        __float128 r = pow( abs(x) , y );
+        __float128 phi = atanq( x.imag() / x.real() );
+        return std::complex< __float128 >( r * cosq( y * phi ) , r * sinq( y * phi ) );
+    }
+}
+
+inline std::ostream& operator<< (std::ostream& os, const __float128& f) {
+
+    char* y = new char[1000];
+    quadmath_snprintf(y, 1000, "%.30Qg", f) ;
+    os.precision(30);
+    os<<y;
+    delete[] y;
+    return os;
+}
+
+
+#include <boost/array.hpp>
+#include <boost/range/algorithm.hpp>
+#include <boost/range/adaptor/filtered.hpp>
+#include <boost/range/numeric.hpp>
+#include <boost/numeric/odeint.hpp>
+
+
+
+using namespace boost::numeric::odeint;
+using namespace std;
+
+typedef __float128 my_float;
+typedef std::vector< std::complex < my_float > > state_type;
+
+struct radMod
+{
+    my_float m_om;
+    my_float m_l;
+
+    radMod( my_float om , my_float l )
+        : m_om( om ) , m_l( l ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , my_float r ) const
+    {
+
+        dxdt[0] = x[1];
+        dxdt[1] = -(2*(r-1)/(r*(r-2)))*x[1]-((m_om*m_om*r*r/((r-2)*(r-2)))-(m_l*(m_l+1)/(r*(r-2))))*x[0];
+    }
+};
+
+
+
+
+
+
+
+int main( int argc , char **argv )
+{
+
+
+    state_type x(2);
+
+    my_float re0 = strtoflt128( "-0.00008944230755601224204687038354994353820468" , NULL );
+    my_float im0 = strtoflt128( "0.00004472229441850588228136889483397204368247" , NULL );
+    my_float re1 = strtoflt128( "-4.464175354293244250869336196695966076150E-6 " , NULL );
+    my_float im1 = strtoflt128( "-8.950483248390306670770345406051469584488E-6" , NULL );
+
+    x[0] = complex< my_float >( re0 ,im0 );
+    x[1] = complex< my_float >( re1 ,im1 );
+
+    const my_float dt =strtoflt128 ("-0.001", NULL);
+    const my_float start =strtoflt128 ("10000.0", NULL);
+    const my_float end =strtoflt128 ("9990.0", NULL);
+    const my_float omega =strtoflt128 ("2.0", NULL);
+    const my_float ell =strtoflt128 ("1.0", NULL);
+
+
+
+    my_float abs_err = strtoflt128( "1.0E-15" , NULL ) , rel_err = strtoflt128( "1.0E-10" , NULL );
+    my_float a_x = strtoflt128( "1.0" , NULL ) , a_dxdt = strtoflt128( "1.0" , NULL );
+
+    typedef runge_kutta_dopri5< state_type, my_float > dopri5_type;
+    typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
+    typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
+    
+    dense_output_dopri5_type dopri5( controlled_dopri5_type( default_error_checker< my_float >( abs_err , rel_err , a_x , a_dxdt ) ) );
+
+    std::for_each( make_adaptive_time_iterator_begin(dopri5 , radMod(omega , ell) , x , start , end , dt) ,
+                   make_adaptive_time_iterator_end(dopri5 , radMod(omega , ell) , x ) ,
+                   []( const std::pair< state_type&, my_float > &x ) {
+                       std::cout << x.second << ", " << x.first[0].real() << "\n"; }
+        );
+
+
+
+    return 0;
+}

libs/numeric/odeint/examples/resizing_lattice.cpp

@@ -0,0 +1,169 @@
+/*
+ * resizing_lattice.cpp
+ *
+ * Demonstrates the usage of resizing of the state type during integration.
+ * Examplary system is a strongly nonlinear, disordered Hamiltonian lattice
+ * where the spreading of energy is investigated
+ *
+ * Copyright 2011-2012 Mario Mulansky
+ * Copyright 2012-2013 Karsten Ahnert
+ * 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)
+ *
+ */
+
+#include <iostream>
+#include <utility>
+
+#include <boost/numeric/odeint.hpp>
+
+#include <boost/ref.hpp>
+#include <boost/random.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+//[ resizing_lattice_system_class
+typedef vector< double > coord_type;
+typedef pair< coord_type , coord_type > state_type;
+
+struct compacton_lattice
+{
+    const int m_max_N;
+    const double m_beta;
+    int m_pot_start_index;
+    vector< double > m_pot;
+
+    compacton_lattice( int max_N , double beta , int pot_start_index )
+        : m_max_N( max_N ) , m_beta( beta ) , m_pot_start_index( pot_start_index ) , m_pot( max_N )
+    {
+        srand( time( NULL ) );
+        // fill random potential with iid values from [0,1]
+        boost::mt19937 rng;
+        boost::uniform_real<> unif( 0.0 , 1.0 );
+        boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
+        generate( m_pot.begin() , m_pot.end() , gen );
+    }
+
+    void operator()( const coord_type &q , coord_type &dpdt )
+    {
+        // calculate dpdt = -dH/dq of this hamiltonian system
+        // dp_i/dt = - V_i * q_i^3 - beta*(q_i - q_{i-1})^3 + beta*(q_{i+1} - q_i)^3
+        const int N = q.size();
+        double diff = q[0] - q[N-1];
+        for( int i=0 ; i<N ; ++i )
+        {
+            dpdt[i] = - m_pot[m_pot_start_index+i] * q[i]*q[i]*q[i] -
+                    m_beta * diff*diff*diff;
+            diff = q[(i+1) % N] - q[i];
+            dpdt[i] += m_beta * diff*diff*diff;
+        }
+    }
+
+    void energy_distribution( const coord_type &q , const coord_type &p , coord_type &energies )
+    {
+        // computes the energy per lattice site normalized by total energy
+        const size_t N = q.size();
+        double en = 0.0;
+        for( size_t i=0 ; i<N ; i++ )
+        {
+            const double diff = q[(i+1) % N] - q[i];
+            energies[i] = p[i]*p[i]/2.0
+                + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i]/4.0
+                + m_beta/4.0 * diff*diff*diff*diff;
+            en += energies[i];
+        }
+        en = 1.0/en;
+        for( size_t i=0 ; i<N ; i++ )
+        {
+            energies[i] *= en;
+        }
+    }
+
+    double energy( const coord_type &q , const coord_type &p )
+    {
+        // calculates the total energy of the excitation
+        const size_t N = q.size();
+        double en = 0.0;
+        for( size_t i=0 ; i<N ; i++ )
+        {
+            const double diff = q[(i+1) % N] - q[i];
+            en += p[i]*p[i]/2.0
+                + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i] / 4.0
+                + m_beta/4.0 * diff*diff*diff*diff;
+        }
+        return en;
+    }
+
+    void change_pot_start( const int delta )
+    {
+        m_pot_start_index += delta;
+    }
+};
+//]
+
+//[ resizing_lattice_resize_function
+void do_resize( coord_type &q , coord_type &p , coord_type &distr , const int N )
+{
+    q.resize( N );
+    p.resize( N );
+    distr.resize( N );
+}
+//]
+
+const int max_N = 1024;
+const double beta = 1.0;
+
+int main()
+{
+    //[ resizing_lattice_initialize
+    //start with 60 sites
+    const int N_start = 60;
+    coord_type q( N_start , 0.0 );
+    q.reserve( max_N );
+    coord_type p( N_start , 0.0 );
+    p.reserve( max_N );
+    // start with uniform momentum distribution over 20 sites
+    fill( p.begin()+20 , p.end()-20 , 1.0/sqrt(20.0) );
+
+    coord_type distr( N_start , 0.0 );
+    distr.reserve( max_N );
+
+    // create the system
+    compacton_lattice lattice( max_N , beta , (max_N-N_start)/2 );
+
+    //create the stepper, note that we use an always_resizer because state size might change during steps
+    typedef symplectic_rkn_sb3a_mclachlan< coord_type , coord_type , double , coord_type , coord_type , double ,
+            range_algebra , default_operations , always_resizer > hamiltonian_stepper;
+    hamiltonian_stepper stepper;
+    hamiltonian_stepper::state_type state = make_pair( q , p );
+    //]
+
+    //[ resizing_lattice_steps_loop
+    double t = 0.0;
+    const double dt = 0.1;
+    const int steps = 10000;
+    for( int step = 0 ; step < steps ; ++step )
+    {
+        stepper.do_step( boost::ref(lattice) , state , t , dt );
+        lattice.energy_distribution( state.first , state.second , distr );
+        if( distr[10] > 1E-150 )
+        {
+            do_resize( state.first , state.second , distr , state.first.size()+20 );
+            rotate( state.first.begin() , state.first.end()-20 , state.first.end() );
+            rotate( state.second.begin() , state.second.end()-20 , state.second.end() );
+            lattice.change_pot_start( -20 );
+            cout << t << ": resized left to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
+        }
+        if( distr[distr.size()-10] > 1E-150 )
+        {
+            do_resize( state.first , state.second , distr , state.first.size()+20 );
+            cout << t << ": resized right to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
+        }
+        t += dt;
+    }
+    //]
+
+    cout << "final lattice size: " << distr.size() << ", final energy: " << lattice.energy( state.first , state.second ) << endl;
+}

libs/numeric/odeint/examples/simple1d.cpp

@@ -0,0 +1,44 @@
+/* Boost libs/numeric/odeint/examples/simple1d.cpp
+
+ Copyright 2012-2013 Mario Mulansky
+ Copyright 2012 Karsten Ahnert
+
+ example for a simple one-dimensional 1st order ODE
+
+ 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)
+ */
+
+
+#include <iostream>
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+/* we solve the simple ODE x' = 3/(2t^2) + x/(2t)
+ * with initial condition x(1) = 0.
+ * Analytic solution is x(t) = sqrt(t) - 1/t
+ */
+
+void rhs( const double x , double &dxdt , const double t )
+{
+    dxdt = 3.0/(2.0*t*t) + x/(2.0*t);
+}
+
+void write_cout( const double &x , const double t )
+{
+    cout << t << '\t' << x << endl;
+}
+
+// state_type = double
+typedef runge_kutta_dopri5< double > stepper_type;
+
+int main()
+{
+    double x = 0.0; //initial value x(1) = 0
+    // use dopri5 with stepsize control and allowed errors 10^-12, integrate t=1...10
+    integrate_adaptive( make_controlled( 1E-12 , 1E-12 , stepper_type() ) , rhs , x , 1.0 , 10.0 , 0.1 , write_cout );
+}

libs/numeric/odeint/examples/solar_system.cpp

@@ -0,0 +1,201 @@
+/* Boost libs/numeric/odeint/examples/solar_system.cpp
+
+ Copyright 2010-2012 Karsten Ahnert
+ Copyright 2011 Mario Mulansky
+
+ Solar system example for Hamiltonian stepper
+
+ 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)
+ */
+
+
+#include <iostream>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+#include "point_type.hpp"
+
+//[ container_type_definition
+// we simulate 5 planets and the sun
+const size_t n = 6;
+
+typedef point< double , 3 > point_type;
+typedef boost::array< point_type , n > container_type;
+typedef boost::array< double , n > mass_type;
+//]
+
+
+
+
+
+
+
+//[ coordinate_function
+const double gravitational_constant = 2.95912208286e-4;
+
+struct solar_system_coor
+{
+    const mass_type &m_masses;
+
+    solar_system_coor( const mass_type &masses ) : m_masses( masses ) { }
+
+    void operator()( const container_type &p , container_type &dqdt ) const
+    {
+        for( size_t i=0 ; i<n ; ++i )
+            dqdt[i] = p[i] / m_masses[i];
+    }
+};
+//]
+
+
+//[ momentum_function
+struct solar_system_momentum
+{
+    const mass_type &m_masses;
+
+    solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }
+
+    void operator()( const container_type &q , container_type &dpdt ) const
+    {
+        const size_t n = q.size();
+        for( size_t i=0 ; i<n ; ++i )
+        {
+            dpdt[i] = 0.0;
+            for( size_t j=0 ; j<i ; ++j )
+            {
+                point_type diff = q[j] - q[i];
+                double d = abs( diff );
+                diff *= ( gravitational_constant * m_masses[i] * m_masses[j] / d / d / d );
+                dpdt[i] += diff;
+                dpdt[j] -= diff;
+
+            }
+        }
+    }
+};
+//]
+
+
+
+
+
+
+
+//[ some_helpers
+point_type center_of_mass( const container_type &x , const mass_type &m )
+{
+    double overall_mass = 0.0;
+    point_type mean( 0.0 );
+    for( size_t i=0 ; i<x.size() ; ++i )
+    {
+        overall_mass += m[i];
+        mean += m[i] * x[i];
+    }
+    if( !x.empty() ) mean /= overall_mass;
+    return mean;
+}
+
+
+double energy( const container_type &q , const container_type &p , const mass_type &masses )
+{
+    const size_t n = q.size();
+    double en = 0.0;
+    for( size_t i=0 ; i<n ; ++i )
+    {
+        en += 0.5 * norm( p[i] ) / masses[i];
+        for( size_t j=0 ; j<i ; ++j )
+        {
+            double diff = abs( q[i] - q[j] );
+            en -= gravitational_constant * masses[j] * masses[i] / diff;
+        }
+    }
+    return en;
+}
+//]
+
+
+//[ streaming_observer
+struct streaming_observer
+{
+    std::ostream& m_out;
+
+    streaming_observer( std::ostream &out ) : m_out( out ) { }
+
+    template< class State >
+    void operator()( const State &x , double t ) const
+    {
+        container_type &q = x.first;
+        m_out << t;
+        for( size_t i=0 ; i<q.size() ; ++i ) m_out << "\t" << q[i];
+        m_out << "\n";
+    }
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+
+    using namespace std;
+    using namespace boost::numeric::odeint;
+
+    mass_type masses = {{
+            1.00000597682 ,      // sun
+            0.000954786104043 ,  // jupiter
+            0.000285583733151 ,  // saturn
+            0.0000437273164546 , // uranus
+            0.0000517759138449 , // neptune
+            1.0 / ( 1.3e8 )      // pluto
+    }};
+
+    container_type q = {{
+            point_type( 0.0 , 0.0 , 0.0 ) ,                        // sun
+            point_type( -3.5023653 , -3.8169847 , -1.5507963 ) ,   // jupiter
+            point_type( 9.0755314 , -3.0458353 , -1.6483708 ) ,    // saturn
+            point_type( 8.3101420 , -16.2901086 , -7.2521278 ) ,   // uranus
+            point_type( 11.4707666 , -25.7294829 , -10.8169456 ) , // neptune
+            point_type( -15.5387357 , -25.2225594 , -3.1902382 )   // pluto
+    }};
+
+    container_type p = {{
+            point_type( 0.0 , 0.0 , 0.0 ) ,                        // sun
+            point_type( 0.00565429 , -0.00412490 , -0.00190589 ) , // jupiter
+            point_type( 0.00168318 , 0.00483525 , 0.00192462 ) ,   // saturn
+            point_type( 0.00354178 , 0.00137102 , 0.00055029 ) ,   // uranus
+            point_type( 0.00288930 , 0.00114527 , 0.00039677 ) ,   // neptune
+            point_type( 0.00276725 , -0.00170702 , -0.00136504 )   // pluto
+    }};
+
+    point_type qmean = center_of_mass( q , masses );
+    point_type pmean = center_of_mass( p , masses );
+    for( size_t i=0 ; i<n ; ++i )
+    {
+        q[i] -= qmean ;
+        p[i] -= pmean;
+    }
+
+    for( size_t i=0 ; i<n ; ++i ) p[i] *= masses[i];
+
+    //[ integration_solar_system
+    typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
+    const double dt = 100.0;
+    
+    integrate_const(
+            stepper_type() ,
+            make_pair( solar_system_coor( masses ) , solar_system_momentum( masses ) ) ,
+            make_pair( boost::ref( q ) , boost::ref( p ) ) ,
+            0.0 , 200000.0 , dt , streaming_observer( cout ) );
+    //]
+
+
+    return 0;
+}
+
+
+/*
+Plot with gnuplot:
+p "solar_system.dat" u 2:4 w l,"solar_system.dat" u 5:7 w l,"solar_system.dat" u 8:10 w l,"solar_system.dat" u 11:13 w l,"solar_system.dat" u 14:16 w l,"solar_system.dat" u 17:19 w l
+ */

libs/numeric/odeint/examples/stepper_details.cpp

@@ -0,0 +1,200 @@
+/*
+ * stepper_details.cpp
+ *
+ * This example demonstrates some details about the steppers in odeint.
+ *
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2012 Mario Mulansky
+ * Copyright 2013 Pascal Germroth
+ *
+ * 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)
+ */
+
+#include <iostream>
+#include <boost/array.hpp>
+#include <boost/bind.hpp>
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+const size_t N = 3;
+
+typedef boost::array< double , N > state_type;
+
+//[ system_function_structure
+void sys( const state_type & /*x*/ , state_type & /*dxdt*/ , const double /*t*/ )
+{
+    // ...
+}
+//]
+
+void sys1( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
+{
+}
+
+void sys2( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ )
+{
+}
+
+
+//[ symplectic_stepper_detail_system_function
+typedef boost::array< double , 1 > vector_type;
+
+
+struct harm_osc_f1
+{
+    void operator()( const vector_type &p , vector_type &dqdt )
+    {
+        dqdt[0] = p[0];
+    }
+};
+
+struct harm_osc_f2
+{
+    void operator()( const vector_type &q , vector_type &dpdt )
+    {
+        dpdt[0] = -q[0];
+    }
+};
+//]
+
+//[ symplectic_stepper_detail_system_class
+struct harm_osc
+{
+    void f1( const vector_type &p , vector_type &dqdt ) const
+    {
+        dqdt[0] = p[0];
+    }
+
+    void f2( const vector_type &q , vector_type &dpdt ) const
+    {
+        dpdt[0] = -q[0];
+    }
+};
+//]
+
+int main( int argc , char **argv )
+{
+    using namespace std;
+
+    // Explicit stepper example
+    {
+        double t( 0.0 ) , dt( 0.1 );
+        state_type in , out , dxdtin , inout;
+        //[ explicit_stepper_detail_example
+        runge_kutta4< state_type > rk;
+        rk.do_step( sys1 , inout , t , dt );               // In-place transformation of inout
+        rk.do_step( sys2 , inout , t , dt );               // call with different system: Ok
+        rk.do_step( sys1 , in , t , out , dt );            // Out-of-place transformation
+        rk.do_step( sys1 , inout , dxdtin , t , dt );      // In-place tranformation of inout
+        rk.do_step( sys1 , in , dxdtin , t , out , dt );   // Out-of-place transformation
+        //]
+    }
+
+
+
+    // FSAL stepper example
+    {
+        double t( 0.0 ) , dt( 0.1 );
+        state_type in , in2 , in3 , out , dxdtin , dxdtout , inout , dxdtinout;
+        //[ fsal_stepper_detail_example
+        runge_kutta_dopri5< state_type > rk;
+        rk.do_step( sys1 , in , t , out , dt );
+        rk.do_step( sys2 , in , t , out , dt );         // DONT do this, sys1 is assumed
+
+        rk.do_step( sys2 , in2 , t , out , dt );
+        rk.do_step( sys2 , in3 , t , out , dt );        // DONT do this, in2 is assumed
+
+        rk.do_step( sys1 , inout , dxdtinout , t , dt );
+        rk.do_step( sys2 , inout , dxdtinout , t , dt );           // Ok, internal derivative is not used, dxdtinout is updated
+
+        rk.do_step( sys1 , in , dxdtin , t , out , dxdtout , dt );
+        rk.do_step( sys2 , in , dxdtin , t , out , dxdtout , dt ); // Ok, internal derivative is not used
+        //]
+    }
+
+
+    // Symplectic harmonic oscillator example
+    {
+        double t( 0.0 ) , dt( 0.1 );
+        //[ symplectic_stepper_detail_example
+        pair< vector_type , vector_type > x;
+        x.first[0] = 1.0; x.second[0] = 0.0;
+        symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
+        rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , x , t , dt );
+        //]
+
+        //[ symplectic_stepper_detail_system_class_example
+        harm_osc h;
+        rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) ,
+                x , t , dt );
+        //]
+    }
+
+    // Simplified harmonic oscillator example
+    {
+        double t = 0.0, dt = 0.1;
+        //[ simplified_symplectic_stepper_example
+        pair< vector_type , vector_type > x;
+        x.first[0] = 1.0; x.second[0] = 0.0;
+        symplectic_rkn_sb3a_mclachlan< vector_type > rkn;
+        rkn.do_step( harm_osc_f1() , x , t , dt );
+        //]
+
+        vector_type q = {{ 1.0 }} , p = {{ 0.0 }};
+        //[ symplectic_stepper_detail_ref_usage
+        rkn.do_step( harm_osc_f1() , make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt );
+        rkn.do_step( harm_osc_f1() , q , p , t , dt );
+        rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , q , p , t , dt );
+        //]
+    }
+    
+    // adams_bashforth_moulton stepper example
+    {
+        double t = 0.0 , dt = 0.1;
+        state_type inout;
+        //[ multistep_detail_example
+        adams_bashforth_moulton< 5 , state_type > abm;
+        abm.initialize( sys , inout , t , dt );
+        abm.do_step( sys , inout , t , dt );
+        //]
+        
+        //[ multistep_detail_own_stepper_initialization
+        abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt );
+        //]
+    }
+
+
+
+    // dense output stepper examples
+    {
+        double t = 0.0 , dt = 0.1;
+        state_type in;
+        //[ dense_output_detail_example
+        dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense;
+        dense.initialize( in , t , dt );
+        pair< double , double > times = dense.do_step( sys );
+        (void)times;
+        //]
+
+        state_type inout;
+        double t_start = 0.0 , t_end = 1.0;
+        //[ dense_output_detail_generation1
+        typedef boost::numeric::odeint::result_of::make_dense_output<
+            runge_kutta_dopri5< state_type > >::type dense_stepper_type;
+        dense_stepper_type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
+        (void)dense2;
+        //]
+
+        //[ dense_output_detail_generation2
+        integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt );
+        //]
+    }
+
+
+    return 0;
+}

libs/numeric/odeint/examples/stiff_system.cpp

@@ -0,0 +1,118 @@
+/*
+ * rosenbrock4.cpp
+ *
+ * Copyright 2010-2012 Mario Mulansky
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2012 Andreas Angelopoulos
+ *
+ * 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)
+ */
+
+#include <iostream>
+#include <fstream>
+#include <utility>
+
+#include <boost/numeric/odeint.hpp>
+
+#include <boost/phoenix/core.hpp>
+
+#include <boost/phoenix/core.hpp>
+#include <boost/phoenix/operator.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+namespace phoenix = boost::phoenix;
+
+
+
+//[ stiff_system_definition
+typedef boost::numeric::ublas::vector< double > vector_type;
+typedef boost::numeric::ublas::matrix< double > matrix_type;
+
+struct stiff_system
+{
+    void operator()( const vector_type &x , vector_type &dxdt , double /* t */ )
+    {
+        dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 * x[ 1 ];
+        dxdt[ 1 ] = x[ 0 ];
+    }
+};
+
+struct stiff_system_jacobi
+{
+    void operator()( const vector_type & /* x */ , matrix_type &J , const double & /* t */ , vector_type &dfdt )
+    {
+        J( 0 , 0 ) = -101.0;
+        J( 0 , 1 ) = -100.0;
+        J( 1 , 0 ) = 1.0;
+        J( 1 , 1 ) = 0.0;
+        dfdt[0] = 0.0;
+        dfdt[1] = 0.0;
+    }
+};
+//]
+
+
+
+/*
+//[ stiff_system_alternative_definition
+typedef boost::numeric::ublas::vector< double > vector_type;
+typedef boost::numeric::ublas::matrix< double > matrix_type;
+
+struct stiff_system
+{
+    template< class State >
+    void operator()( const State &x , State &dxdt , double t )
+    {
+        ...
+    }
+};
+
+struct stiff_system_jacobi
+{
+    template< class State , class Matrix >
+    void operator()( const State &x , Matrix &J , const double &t , State &dfdt )
+    {
+        ...
+    }
+};
+//]
+ */
+
+
+
+int main( int argc , char **argv )
+{
+//    typedef rosenbrock4< double > stepper_type;
+//    typedef rosenbrock4_controller< stepper_type > controlled_stepper_type;
+//    typedef rosenbrock4_dense_output< controlled_stepper_type > dense_output_type;
+    //[ integrate_stiff_system
+    vector_type x( 2 , 1.0 );
+
+    size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< double > >( 1.0e-6 , 1.0e-6 ) ,
+            make_pair( stiff_system() , stiff_system_jacobi() ) ,
+            x , 0.0 , 50.0 , 0.01 ,
+            cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
+    //]
+    clog << num_of_steps << endl;
+
+
+
+//    typedef runge_kutta_dopri5< vector_type > dopri5_type;
+//    typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
+//    typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
+    //[ integrate_stiff_system_alternative
+
+    vector_type x2( 2 , 1.0 );
+
+    size_t num_of_steps2 = integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) ,
+            stiff_system() , x2 , 0.0 , 50.0 , 0.01 ,
+            cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << "\n" );
+    //]
+    clog << num_of_steps2 << endl;
+
+
+    return 0;
+}

libs/numeric/odeint/examples/stochastic_euler.cpp

@@ -0,0 +1,147 @@
+/*
+ libs/numeric/odeint/examples/stochastic_euler.hpp
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 Mario Mulansky
+
+ Stochastic euler stepper example and Ornstein-Uhlenbeck process
+
+ 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)
+ */
+
+
+#include <vector>
+#include <iostream>
+#include <boost/random.hpp>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+
+/*
+//[ stochastic_euler_class_definition
+template< size_t N > class stochastic_euler
+{
+public:
+
+    typedef boost::array< double , N > state_type;
+    typedef boost::array< double , N > deriv_type;
+    typedef double value_type;
+    typedef double time_type;
+    typedef unsigned short order_type;
+    typedef boost::numeric::odeint::stepper_tag stepper_category;
+
+    static order_type order( void ) { return 1; }
+
+    // ...
+};
+//]
+*/
+
+
+/*
+//[ stochastic_euler_do_step
+template< size_t N > class stochastic_euler
+{
+public:
+
+    // ...
+
+    template< class System >
+    void do_step( System system , state_type &x , time_type t , time_type dt ) const
+    {
+        deriv_type det , stoch ;
+        system.first( x , det );
+        system.second( x , stoch );
+        for( size_t i=0 ; i<x.size() ; ++i )
+            x[i] += dt * det[i] + sqrt( dt ) * stoch[i];
+    }
+};
+//]
+*/
+
+
+
+
+//[ stochastic_euler_class
+template< size_t N >
+class stochastic_euler
+{
+public:
+
+    typedef boost::array< double , N > state_type;
+    typedef boost::array< double , N > deriv_type;
+    typedef double value_type;
+    typedef double time_type;
+    typedef unsigned short order_type;
+
+    typedef boost::numeric::odeint::stepper_tag stepper_category;
+
+    static order_type order( void ) { return 1; }
+
+    template< class System >
+    void do_step( System system , state_type &x , time_type t , time_type dt ) const
+    {
+        deriv_type det , stoch ;
+        system.first( x , det );
+        system.second( x , stoch );
+        for( size_t i=0 ; i<x.size() ; ++i )
+            x[i] += dt * det[i] + sqrt( dt ) * stoch[i];
+    }
+};
+//]
+
+
+
+//[ stochastic_euler_ornstein_uhlenbeck_def
+const static size_t N = 1;
+typedef boost::array< double , N > state_type;
+
+struct ornstein_det
+{
+    void operator()( const state_type &x , state_type &dxdt ) const
+    {
+        dxdt[0] = -x[0];
+    }
+};
+
+struct ornstein_stoch
+{
+    boost::mt19937 &m_rng;
+    boost::normal_distribution<> m_dist;
+
+  ornstein_stoch( boost::mt19937 &rng , double sigma ) : m_rng( rng ) , m_dist( 0.0 , sigma ) { }
+
+    void operator()( const state_type &x , state_type &dxdt )
+    {
+        dxdt[0] = m_dist( m_rng );
+    }
+};
+//]
+
+struct streaming_observer
+{
+    template< class State >
+    void operator()( const State &x , double t ) const
+    {
+        std::cout << t << "\t" << x[0] << "\n";
+    }
+};
+
+
+int main( int argc , char **argv )
+{
+    using namespace std;
+    using namespace boost::numeric::odeint;
+
+    //[ ornstein_uhlenbeck_main
+    boost::mt19937 rng;
+    double dt = 0.1;
+    state_type x = {{ 1.0 }};
+    integrate_const( stochastic_euler< N >() , make_pair( ornstein_det() , ornstein_stoch( rng , 1.0 ) ),
+            x , 0.0 , 10.0 , dt , streaming_observer() );
+    //]
+    return 0;
+}

libs/numeric/odeint/examples/stuart_landau.cpp

@@ -0,0 +1,87 @@
+/*
+ * stuart_landau.cpp
+ *
+ * This example demonstrates how one can use odeint can be used with state types consisting of complex variables.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2011-2013 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)
+ */
+
+#include <iostream>
+#include <complex>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+//[ stuart_landau_system_function
+typedef complex< double > state_type;
+
+struct stuart_landau
+{
+    double m_eta;
+    double m_alpha;
+
+    stuart_landau( double eta = 1.0 , double alpha = 1.0 )
+    : m_eta( eta ) , m_alpha( alpha ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const
+    {
+        const complex< double > I( 0.0 , 1.0 );
+        dxdt = ( 1.0 + m_eta * I ) * x - ( 1.0 + m_alpha * I ) * norm( x ) * x;
+    }
+};
+//]
+
+
+/*
+//[ stuart_landau_system_function_alternative
+double eta = 1.0;
+double alpha = 1.0;
+
+void stuart_landau( const state_type &x , state_type &dxdt , double t )
+{
+    const complex< double > I( 0.0 , 1.0 );
+    dxdt = ( 1.0 + m_eta * I ) * x - ( 1.0 + m_alpha * I ) * norm( x ) * x;
+}
+//]
+*/
+
+
+struct streaming_observer
+{
+    std::ostream& m_out;
+
+    streaming_observer( std::ostream &out ) : m_out( out ) { }
+
+    template< class State >
+    void operator()( const State &x , double t ) const
+    {
+        m_out << t;
+        m_out << "\t" << x.real() << "\t" << x.imag() ;
+        m_out << "\n";
+    }
+};
+
+
+
+
+int main( int argc , char **argv )
+{
+    //[ stuart_landau_integration
+    state_type x = complex< double >( 1.0 , 0.0 );
+
+    const double dt = 0.1;
+
+    typedef runge_kutta4< state_type > stepper_type;
+
+    integrate_const( stepper_type() , stuart_landau( 2.0 , 1.0 ) , x , 0.0 , 10.0 , dt , streaming_observer( cout ) );
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/thrust/Makefile

@@ -0,0 +1,34 @@
+# Copyright 2011-2014 Mario Mulansky
+# Copyright 2011-2012 Karsten Ahnert
+#
+# 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)
+
+# make sure BOOST_ROOT is pointing to your boost directory
+# otherwise, set it here:
+# BOOST_ROOT = /path/to/boost
+
+# path to the cuda installation
+CUDA_ROOT = /usr/local/cuda
+# target architecture
+ARCH = sm_13
+
+NVCC = $(CUDA_ROOT)/bin/nvcc
+
+INCLUDES += -I../../include/ -I$(BOOST_ROOT)
+
+NVCCFLAGS = -O3 $(INCLUDES) -arch $(ARCH)
+
+%.o : %.cu
+	$(NVCC) $(NVCCFLAGS) -c $< -o $@
+
+% : %.o
+	$(NVCC) $(NVCCFLAGS) -o $@ $<
+
+
+all : phase_oscillator_chain phase_oscillator_ensemble lorenz_parameters relaxation
+
+
+clean :
+	-rm *~ *.o phase_oscillator_chain phase_oscillator_ensemble lorenz_parameters relaxation 

libs/numeric/odeint/examples/thrust/lorenz_parameters.cu

@@ -0,0 +1,296 @@
+/*
+ Copyright 2011-2012 Karsten Ahnert
+ Copyright 2011-2013 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)
+ */
+
+#include <iostream>
+#include <cmath>
+#include <utility>
+
+
+#include <thrust/device_vector.h>
+#include <thrust/reduce.h>
+#include <thrust/functional.h>
+
+#include <boost/numeric/odeint.hpp>
+
+#include <boost/numeric/odeint/external/thrust/thrust.hpp>
+
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/variate_generator.hpp>
+
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+//change this to float if your device does not support double computation
+typedef double value_type;
+
+//change this to host_vector< ... > of you want to run on CPU
+typedef thrust::device_vector< value_type > state_type;
+typedef thrust::device_vector< size_t > index_vector_type;
+// typedef thrust::host_vector< value_type > state_type;
+// typedef thrust::host_vector< size_t > index_vector_type;
+
+
+const value_type sigma = 10.0;
+const value_type b = 8.0 / 3.0;
+
+
+//[ thrust_lorenz_parameters_define_simple_system
+struct lorenz_system
+{
+    struct lorenz_functor
+    {
+        template< class T >
+        __host__ __device__
+        void operator()( T t ) const
+        {
+            // unpack the parameter we want to vary and the Lorenz variables
+            value_type R = thrust::get< 3 >( t );
+            value_type x = thrust::get< 0 >( t );
+            value_type y = thrust::get< 1 >( t );
+            value_type z = thrust::get< 2 >( t );
+            thrust::get< 4 >( t ) = sigma * ( y - x );
+            thrust::get< 5 >( t ) = R * x - y - x * z;
+            thrust::get< 6 >( t ) = -b * z + x * y ;
+
+        }
+    };
+
+    lorenz_system( size_t N , const state_type &beta )
+    : m_N( N ) , m_beta( beta ) { }
+
+    template< class State , class Deriv >
+    void operator()(  const State &x , Deriv &dxdt , value_type t ) const
+    {
+        thrust::for_each(
+                thrust::make_zip_iterator( thrust::make_tuple(
+                        boost::begin( x ) ,
+                        boost::begin( x ) + m_N ,
+                        boost::begin( x ) + 2 * m_N ,
+                        m_beta.begin() ,
+                        boost::begin( dxdt ) ,
+                        boost::begin( dxdt ) + m_N ,
+                        boost::begin( dxdt ) + 2 * m_N  ) ) ,
+                thrust::make_zip_iterator( thrust::make_tuple(
+                        boost::begin( x ) + m_N ,
+                        boost::begin( x ) + 2 * m_N ,
+                        boost::begin( x ) + 3 * m_N ,
+                        m_beta.begin() ,
+                        boost::begin( dxdt ) + m_N ,
+                        boost::begin( dxdt ) + 2 * m_N ,
+                        boost::begin( dxdt ) + 3 * m_N  ) ) ,
+                lorenz_functor() );
+    }
+    size_t m_N;
+    const state_type &m_beta;
+};
+//]
+
+struct lorenz_perturbation_system
+{
+    struct lorenz_perturbation_functor
+    {
+        template< class T >
+        __host__ __device__
+        void operator()( T t ) const
+        {
+            value_type R = thrust::get< 1 >( t );
+            value_type x = thrust::get< 0 >( thrust::get< 0 >( t ) );
+            value_type y = thrust::get< 1 >( thrust::get< 0 >( t ) );
+            value_type z = thrust::get< 2 >( thrust::get< 0 >( t ) );
+            value_type dx = thrust::get< 3 >( thrust::get< 0 >( t ) );
+            value_type dy = thrust::get< 4 >( thrust::get< 0 >( t ) );
+            value_type dz = thrust::get< 5 >( thrust::get< 0 >( t ) );
+            thrust::get< 0 >( thrust::get< 2 >( t ) ) = sigma * ( y - x );
+            thrust::get< 1 >( thrust::get< 2 >( t ) ) = R * x - y - x * z;
+            thrust::get< 2 >( thrust::get< 2 >( t ) ) = -b * z + x * y ;
+            thrust::get< 3 >( thrust::get< 2 >( t ) ) = sigma * ( dy - dx );
+            thrust::get< 4 >( thrust::get< 2 >( t ) ) = ( R - z ) * dx - dy - x * dz;
+            thrust::get< 5 >( thrust::get< 2 >( t ) ) = y * dx + x * dy - b * dz;
+        }
+    };
+
+    lorenz_perturbation_system( size_t N , const state_type &beta )
+    : m_N( N ) , m_beta( beta ) { }
+
+    template< class State , class Deriv >
+    void operator()(  const State &x , Deriv &dxdt , value_type t ) const
+    {
+        thrust::for_each(
+                thrust::make_zip_iterator( thrust::make_tuple(
+                        thrust::make_zip_iterator( thrust::make_tuple(
+                                boost::begin( x ) ,
+                                boost::begin( x ) + m_N ,
+                                boost::begin( x ) + 2 * m_N ,
+                                boost::begin( x ) + 3 * m_N ,
+                                boost::begin( x ) + 4 * m_N ,
+                                boost::begin( x ) + 5 * m_N ) ) ,
+                        m_beta.begin() ,
+                        thrust::make_zip_iterator( thrust::make_tuple(
+                                boost::begin( dxdt ) ,
+                                boost::begin( dxdt ) + m_N ,
+                                boost::begin( dxdt ) + 2 * m_N ,
+                                boost::begin( dxdt ) + 3 * m_N ,
+                                boost::begin( dxdt ) + 4 * m_N ,
+                                boost::begin( dxdt ) + 5 * m_N ) )
+                ) ) ,
+                thrust::make_zip_iterator( thrust::make_tuple(
+                        thrust::make_zip_iterator( thrust::make_tuple(
+                                boost::begin( x ) + m_N ,
+                                boost::begin( x ) + 2 * m_N ,
+                                boost::begin( x ) + 3 * m_N ,
+                                boost::begin( x ) + 4 * m_N ,
+                                boost::begin( x ) + 5 * m_N ,
+                                boost::begin( x ) + 6 * m_N ) ) ,
+                        m_beta.begin() ,
+                        thrust::make_zip_iterator( thrust::make_tuple(
+                                boost::begin( dxdt ) + m_N ,
+                                boost::begin( dxdt ) + 2 * m_N ,
+                                boost::begin( dxdt ) + 3 * m_N ,
+                                boost::begin( dxdt ) + 4 * m_N ,
+                                boost::begin( dxdt ) + 5 * m_N ,
+                                boost::begin( dxdt ) + 6 * m_N  ) )
+                ) ) ,
+                lorenz_perturbation_functor() );
+    }
+
+    size_t m_N;
+    const state_type &m_beta;
+};
+
+struct lyap_observer
+{
+    //[thrust_lorenz_parameters_observer_functor
+    struct lyap_functor
+    {
+        template< class T >
+        __host__ __device__
+        void operator()( T t ) const
+        {
+            value_type &dx = thrust::get< 0 >( t );
+            value_type &dy = thrust::get< 1 >( t );
+            value_type &dz = thrust::get< 2 >( t );
+            value_type norm = sqrt( dx * dx + dy * dy + dz * dz );
+            dx /= norm;
+            dy /= norm;
+            dz /= norm;
+            thrust::get< 3 >( t ) += log( norm );
+        }
+    };
+    //]
+
+    lyap_observer( size_t N , size_t every = 100 )
+    : m_N( N ) , m_lyap( N ) , m_every( every ) , m_count( 0 )
+    {
+        thrust::fill( m_lyap.begin() , m_lyap.end() , 0.0 );
+    }
+
+    template< class Lyap >
+    void fill_lyap( Lyap &lyap )
+    {
+        thrust::copy( m_lyap.begin() , m_lyap.end() , lyap.begin() );
+        for( size_t i=0 ; i<lyap.size() ; ++i )
+            lyap[i] /= m_t_overall;
+    }
+
+
+    template< class State >
+    void operator()( State &x , value_type t )
+    {
+        if( ( m_count != 0 ) && ( ( m_count % m_every ) == 0 ) )
+        {
+            thrust::for_each(
+                    thrust::make_zip_iterator( thrust::make_tuple(
+                            boost::begin( x ) + 3 * m_N ,
+                            boost::begin( x ) + 4 * m_N ,
+                            boost::begin( x ) + 5 * m_N ,
+                            m_lyap.begin() ) ) ,
+                    thrust::make_zip_iterator( thrust::make_tuple(
+                            boost::begin( x ) + 4 * m_N ,
+                            boost::begin( x ) + 5 * m_N ,
+                            boost::begin( x ) + 6 * m_N ,
+                            m_lyap.end() ) ) ,
+                    lyap_functor() );
+            clog << t << "\n";
+        }
+        ++m_count;
+        m_t_overall = t;
+    }
+
+    size_t m_N;
+    state_type m_lyap;
+    size_t m_every;
+    size_t m_count;
+    value_type m_t_overall;
+};
+
+const size_t N = 1024*2;
+const value_type dt = 0.01;
+
+
+int main( int arc , char* argv[] )
+{
+    int driver_version , runtime_version;
+    cudaDriverGetVersion( &driver_version );
+    cudaRuntimeGetVersion ( &runtime_version );
+    cout << driver_version << "\t" << runtime_version << endl;
+
+
+    //[ thrust_lorenz_parameters_define_beta
+    vector< value_type > beta_host( N );
+    const value_type beta_min = 0.0 , beta_max = 56.0;
+    for( size_t i=0 ; i<N ; ++i )
+        beta_host[i] = beta_min + value_type( i ) * ( beta_max - beta_min ) / value_type( N - 1 );
+
+    state_type beta = beta_host;
+    //]
+
+    //[ thrust_lorenz_parameters_integration
+    state_type x( 6 * N );
+
+    // initialize x,y,z
+    thrust::fill( x.begin() , x.begin() + 3 * N , 10.0 );
+
+    // initial dx
+    thrust::fill( x.begin() + 3 * N , x.begin() + 4 * N , 1.0 );
+
+    // initialize dy,dz
+    thrust::fill( x.begin() + 4 * N , x.end() , 0.0 );
+
+
+    // create error stepper, can be used with make_controlled or make_dense_output
+    typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type > stepper_type;
+
+
+    lorenz_system lorenz( N , beta );
+    lorenz_perturbation_system lorenz_perturbation( N , beta );
+    lyap_observer obs( N , 1 );
+
+    // calculate transients
+    integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz , std::make_pair( x.begin() , x.begin() + 3 * N ) , 0.0 , 10.0 , dt );
+
+    // calculate the Lyapunov exponents -- the main loop
+    double t = 0.0;
+    while( t < 10000.0 )
+    {
+        integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz_perturbation , x , t , t + 1.0 , 0.1 );
+        t += 1.0;
+        obs( x , t );
+    }
+
+    vector< value_type > lyap( N );
+    obs.fill_lyap( lyap );
+
+    for( size_t i=0 ; i<N ; ++i )
+        cout << beta_host[i] << "\t" << lyap[i] << "\n";
+    //]
+
+    return 0;
+}

libs/numeric/odeint/examples/thrust/phase_oscillator_chain.cu

@@ -0,0 +1,156 @@
+/*
+ Copyright 2011-2013 Mario Mulansky
+ Copyright 2011 Karsten Ahnert
+
+ 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)
+ */
+
+/*
+ * This example shows how to use odeint on CUDA devices with thrust.
+ * Note that we require at least Version 3.2 of the nVidia CUDA SDK
+ * and the thrust library should be installed in the CUDA include
+ * folder.
+ *
+ * As example we use a chain of phase oscillators with nearest neighbour
+ * coupling, as described in:
+ *
+ * Avis H. Cohen, Philip J. Holmes and Richard H. Rand:
+ * JOURNAL OF MATHEMATICAL BIOLOGY Volume 13, Number 3, 345-369,
+ *
+ */
+
+#include <iostream>
+#include <cmath>
+
+#include <thrust/device_vector.h>
+#include <thrust/iterator/permutation_iterator.h>
+#include <thrust/iterator/counting_iterator.h>
+
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/integrate/integrate_const.hpp>
+#include <boost/numeric/odeint/external/thrust/thrust.hpp>
+
+using namespace std;
+
+using namespace boost::numeric::odeint;
+
+
+//change this to float if your device does not support double computation
+typedef double value_type;
+
+
+//[ thrust_phase_chain_system
+//change this to host_vector< ... > if you want to run on CPU
+typedef thrust::device_vector< value_type > state_type;
+typedef thrust::device_vector< size_t > index_vector_type;
+//typedef thrust::host_vector< value_type > state_type;
+//typedef thrust::host_vector< size_t > index_vector_type;
+
+//<-
+/*
+ * This implements the rhs of the dynamical equation:
+ * \phi'_0 = \omega_0 + sin( \phi_1 - \phi_0 )
+ * \phi'_i  = \omega_i + sin( \phi_i+1 - \phi_i ) + sin( \phi_i - \phi_i-1 )
+ * \phi'_N-1 = \omega_N-1 + sin( \phi_N-1 - \phi_N-2 )
+ */
+//->
+class phase_oscillators
+{
+
+public:
+
+    struct sys_functor
+    {
+        template< class Tuple >
+        __host__ __device__
+        void operator()( Tuple t )  // this functor works on tuples of values
+        {
+            // first, unpack the tuple into value, neighbors and omega
+            const value_type phi = thrust::get<0>(t);
+            const value_type phi_left = thrust::get<1>(t);  // left neighbor
+            const value_type phi_right = thrust::get<2>(t); // right neighbor
+            const value_type omega = thrust::get<3>(t);
+            // the dynamical equation
+            thrust::get<4>(t) = omega + sin( phi_right - phi ) + sin( phi - phi_left );
+        }
+    };
+
+    phase_oscillators( const state_type &omega )
+        : m_omega( omega ) , m_N( omega.size() ) , m_prev( omega.size() ) , m_next( omega.size() )
+    {
+        // build indices pointing to left and right neighbours
+        thrust::counting_iterator<size_t> c( 0 );
+        thrust::copy( c , c+m_N-1 , m_prev.begin()+1 );
+        m_prev[0] = 0; // m_prev = { 0 , 0 , 1 , 2 , 3 , ... , N-1 }
+
+        thrust::copy( c+1 , c+m_N , m_next.begin() );
+        m_next[m_N-1] = m_N-1; // m_next = { 1 , 2 , 3 , ... , N-1 , N-1 }
+    }
+
+    void operator() ( const state_type &x , state_type &dxdt , const value_type dt )
+    {
+        thrust::for_each(
+                thrust::make_zip_iterator(
+                        thrust::make_tuple(
+                                x.begin() ,
+                                thrust::make_permutation_iterator( x.begin() , m_prev.begin() ) ,
+                                thrust::make_permutation_iterator( x.begin() , m_next.begin() ) ,
+                                m_omega.begin() ,
+                                dxdt.begin()
+                                ) ),
+                thrust::make_zip_iterator(
+                        thrust::make_tuple(
+                                x.end() ,
+                                thrust::make_permutation_iterator( x.begin() , m_prev.end() ) ,
+                                thrust::make_permutation_iterator( x.begin() , m_next.end() ) ,
+                                m_omega.end() ,
+                                dxdt.end()) ) ,
+                sys_functor()
+                );
+    }
+
+private:
+
+    const state_type &m_omega;
+    const size_t m_N;
+    index_vector_type m_prev;
+    index_vector_type m_next;
+};
+//]
+
+const size_t N = 32768;
+const value_type pi = 3.1415926535897932384626433832795029;
+const value_type epsilon = 6.0 / ( N * N ); // should be < 8/N^2 to see phase locking
+const value_type dt = 0.1;
+
+int main( int arc , char* argv[] )
+{
+    //[ thrust_phase_chain_integration
+    // create initial conditions and omegas on host:
+    vector< value_type > x_host( N );
+    vector< value_type > omega_host( N );
+    for( size_t i=0 ; i<N ; ++i )
+    {
+        x_host[i] = 2.0 * pi * drand48();
+        omega_host[i] = ( N - i ) * epsilon; // decreasing frequencies
+    }
+
+    // copy to device
+    state_type x = x_host;
+    state_type omega = omega_host;
+
+    // create stepper
+    runge_kutta4< state_type , value_type , state_type , value_type > stepper;
+
+    // create phase oscillator system function
+    phase_oscillators sys( omega );
+
+    // integrate
+    integrate_const( stepper , sys , x , 0.0 , 10.0 , dt );
+
+    thrust::copy( x.begin() , x.end() , std::ostream_iterator< value_type >( std::cout , "\n" ) );
+    std::cout << std::endl;
+    //]
+}

libs/numeric/odeint/examples/thrust/phase_oscillator_ensemble.cu

@@ -0,0 +1,280 @@
+/*
+ The example how the phase_oscillator ensemble can be implemented using CUDA and thrust
+
+ Copyright 2011-2013 Mario Mulansky
+ Copyright 2011 Karsten Ahnert
+
+ 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)
+ */
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <utility>
+
+#include <thrust/device_vector.h>
+#include <thrust/reduce.h>
+#include <thrust/functional.h>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/thrust/thrust.hpp>
+
+#include <boost/timer.hpp>
+#include <boost/random/cauchy_distribution.hpp>
+
+using namespace std;
+
+using namespace boost::numeric::odeint;
+
+/*
+ * Sorry for that dirty hack, but nvcc has large problems with boost::random.
+ *
+ * Nevertheless we need the cauchy distribution from boost::random, and therefore
+ * we need a generator. Here it is:
+ */
+struct drand48_generator
+{
+    typedef double result_type;
+    result_type operator()( void ) const { return drand48(); }
+    result_type min( void ) const { return 0.0; }
+    result_type max( void ) const { return 1.0; }
+};
+
+//[ thrust_phase_ensemble_state_type
+//change this to float if your device does not support double computation
+typedef double value_type;
+
+//change this to host_vector< ... > of you want to run on CPU
+typedef thrust::device_vector< value_type > state_type;
+// typedef thrust::host_vector< value_type > state_type;
+//]
+
+
+//[ thrust_phase_ensemble_mean_field_calculator
+struct mean_field_calculator
+{
+    struct sin_functor : public thrust::unary_function< value_type , value_type >
+    {
+        __host__ __device__
+        value_type operator()( value_type x) const
+        {
+            return sin( x );
+        }
+    };
+
+    struct cos_functor : public thrust::unary_function< value_type , value_type >
+    {
+        __host__ __device__
+        value_type operator()( value_type x) const
+        {
+            return cos( x );
+        }
+    };
+
+    static std::pair< value_type , value_type > get_mean( const state_type &x )
+    {
+        //[ thrust_phase_ensemble_sin_sum
+        value_type sin_sum = thrust::reduce(
+                thrust::make_transform_iterator( x.begin() , sin_functor() ) ,
+                thrust::make_transform_iterator( x.end() , sin_functor() ) );
+        //]
+        value_type cos_sum = thrust::reduce(
+                thrust::make_transform_iterator( x.begin() , cos_functor() ) ,
+                thrust::make_transform_iterator( x.end() , cos_functor() ) );
+
+        cos_sum /= value_type( x.size() );
+        sin_sum /= value_type( x.size() );
+
+        value_type K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum );
+        value_type Theta = atan2( sin_sum , cos_sum );
+
+        return std::make_pair( K , Theta );
+    }
+};
+//]
+
+
+
+//[ thrust_phase_ensemble_sys_function
+class phase_oscillator_ensemble
+{
+
+public:
+
+    struct sys_functor
+    {
+        value_type m_K , m_Theta , m_epsilon;
+
+        sys_functor( value_type K , value_type Theta , value_type epsilon )
+        : m_K( K ) , m_Theta( Theta ) , m_epsilon( epsilon ) { }
+
+        template< class Tuple >
+        __host__ __device__
+        void operator()( Tuple t )
+        {
+            thrust::get<2>(t) = thrust::get<1>(t) + m_epsilon * m_K * sin( m_Theta - thrust::get<0>(t) );
+        }
+    };
+
+    // ...
+    //<-
+    phase_oscillator_ensemble( size_t N , value_type g = 1.0 , value_type epsilon = 1.0 )
+        : m_omega() , m_N( N ) , m_epsilon( epsilon )
+    {
+        create_frequencies( g );
+    }
+
+    void create_frequencies( value_type g )
+    {
+        boost::cauchy_distribution< value_type > cauchy( 0.0 , g );
+//        boost::variate_generator< boost::mt19937&, boost::cauchy_distribution< value_type > > gen( rng , cauchy );
+        drand48_generator d48;
+        vector< value_type > omega( m_N );
+        for( size_t i=0 ; i<m_N ; ++i )
+            omega[i] = cauchy( d48 );
+//        generate( omega.begin() , omega.end() , gen );
+        m_omega = omega;
+    }
+
+    void set_epsilon( value_type epsilon ) { m_epsilon = epsilon; }
+
+    value_type get_epsilon( void ) const { return m_epsilon; }
+    //->
+
+    void operator() ( const state_type &x , state_type &dxdt , const value_type dt ) const
+    {
+        std::pair< value_type , value_type > mean_field = mean_field_calculator::get_mean( x );
+
+        thrust::for_each(
+                thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_omega.begin() , dxdt.begin() ) ),
+                thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_omega.end() , dxdt.end()) ) ,
+                sys_functor( mean_field.first , mean_field.second , m_epsilon )
+                );
+    }
+
+    // ...
+    //<-
+private:
+
+    state_type m_omega;
+    const size_t m_N;
+    value_type m_epsilon;
+    //->
+};
+//]
+
+
+//[ thrust_phase_ensemble_observer
+struct statistics_observer
+{
+    value_type m_K_mean;
+    size_t m_count;
+
+    statistics_observer( void )
+    : m_K_mean( 0.0 ) , m_count( 0 ) { }
+
+    template< class State >
+    void operator()( const State &x , value_type t )
+    {
+        std::pair< value_type , value_type > mean = mean_field_calculator::get_mean( x );
+        m_K_mean += mean.first;
+        ++m_count;
+    }
+
+    value_type get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / value_type( m_count ) : 0.0 ; }
+
+    void reset( void ) { m_K_mean = 0.0; m_count = 0; }
+};
+//]
+
+
+
+// const size_t N = 16384 * 128;
+const size_t N = 16384;
+const value_type pi = 3.1415926535897932384626433832795029;
+const value_type dt = 0.1;
+const value_type d_epsilon = 0.1;
+const value_type epsilon_min = 0.0;
+const value_type epsilon_max = 5.0;
+const value_type t_transients = 10.0;
+const value_type t_max = 100.0;
+
+int main( int arc , char* argv[] )
+{
+    // initial conditions on host
+    vector< value_type > x_host( N );
+    for( size_t i=0 ; i<N ; ++i ) x_host[i] = 2.0 * pi * drand48();
+
+    //[ thrust_phase_ensemble_system_instance
+    phase_oscillator_ensemble ensemble( N , 1.0 );
+    //]
+
+
+
+    boost::timer timer;
+    boost::timer timer_local;
+    double dopri5_time = 0.0 , rk4_time = 0.0;
+    {
+        //[thrust_phase_ensemble_define_dopri5
+        typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type > stepper_type;
+        //]
+
+        ofstream fout( "phase_ensemble_dopri5.dat" );
+        timer.restart();
+        for( value_type epsilon = epsilon_min ; epsilon < epsilon_max ; epsilon += d_epsilon )
+        {
+            ensemble.set_epsilon( epsilon );
+            statistics_observer obs;
+            state_type x = x_host;
+
+            timer_local.restart();
+
+            // calculate some transients steps
+            //[ thrust_phase_ensemble_integration
+            size_t steps1 = integrate_const( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , boost::ref( ensemble ) , x , 0.0 , t_transients , dt );
+            //]
+
+            // integrate and compute the statistics
+            size_t steps2 = integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , boost::ref( ensemble ) , x , 0.0 , t_max , dt , boost::ref( obs ) );
+
+            fout << epsilon << "\t" << obs.get_K_mean() << endl;
+            cout << "Dopri5 : " << epsilon << "\t" << obs.get_K_mean() << "\t" << timer_local.elapsed() << "\t" << steps1 << "\t" << steps2 << endl;
+        }
+        dopri5_time = timer.elapsed();
+    }
+
+
+
+    {
+        //[ thrust_phase_ensemble_define_rk4
+        typedef runge_kutta4< state_type , value_type , state_type , value_type > stepper_type;
+        //]
+
+        ofstream fout( "phase_ensemble_rk4.dat" );
+        timer.restart();
+        for( value_type epsilon = epsilon_min ; epsilon < epsilon_max ; epsilon += d_epsilon )
+        {
+            ensemble.set_epsilon( epsilon );
+            statistics_observer obs;
+            state_type x = x_host;
+
+            timer_local.restart();
+
+            // calculate some transients steps
+            size_t steps1 = integrate_const( stepper_type() , boost::ref( ensemble ) , x , 0.0 , t_transients , dt );
+
+            // integrate and compute the statistics
+            size_t steps2 = integrate_const( stepper_type() , boost::ref( ensemble ) , x , 0.0 , t_max , dt , boost::ref( obs ) );
+            fout << epsilon << "\t" << obs.get_K_mean() << endl;
+            cout << "RK4     : " << epsilon << "\t" << obs.get_K_mean() << "\t" << timer_local.elapsed() << "\t" << steps1 << "\t" << steps2 << endl;
+        }
+        rk4_time = timer.elapsed();
+    }
+
+    cout << "Dopri 5 : " << dopri5_time << " s\n";
+    cout << "RK4     : " << rk4_time << "\n";
+
+    return 0;
+}

libs/numeric/odeint/examples/thrust/relaxation.cu

@@ -0,0 +1,81 @@
+/*
+ Copyright 2011-2012 Karsten Ahnert
+ Copyright 2013 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)
+ */
+
+
+/*
+ * Solves many relaxation equations dxdt = - a * x in parallel and for different values of a.
+ * The relaxation equations are completely uncoupled.
+ */
+
+#include <thrust/device_vector.h>
+
+#include <boost/ref.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/thrust/thrust.hpp>
+
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+// change to float if your GPU does not support doubles
+typedef double value_type;
+typedef thrust::device_vector< value_type > state_type;
+typedef runge_kutta4< state_type , value_type , state_type , value_type > stepper_type;
+
+struct relaxation
+{
+    struct relaxation_functor
+    {
+        template< class T >
+        __host__ __device__
+        void operator()( T t ) const
+        {
+            // unpack the parameter we want to vary and the Lorenz variables
+            value_type a = thrust::get< 1 >( t );
+            value_type x = thrust::get< 0 >( t );
+            thrust::get< 2 >( t ) = -a * x;
+        }
+    };
+
+    relaxation( size_t N , const state_type &a )
+    : m_N( N ) , m_a( a ) { }
+
+    void operator()(  const state_type &x , state_type &dxdt , value_type t ) const
+    {
+        thrust::for_each(
+            thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_a.begin() , dxdt.begin() ) ) ,
+            thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_a.end() , dxdt.end() ) ) ,
+            relaxation_functor() );
+    }
+
+    size_t m_N;
+    const state_type &m_a;
+};
+
+const size_t N = 1024 * 1024;
+const value_type dt = 0.01;
+
+int main( int arc , char* argv[] )
+{
+    // initialize the relaxation constants a
+    vector< value_type > a_host( N );
+    for( size_t i=0 ; i<N ; ++i ) a_host[i] = drand48();
+    state_type a = a_host;
+
+    // initialize the intial state x
+    state_type x( N );
+    thrust::fill( x.begin() , x.end() , 1.0 );
+
+    // integrate
+    relaxation relax( N , a );
+    integrate_const( stepper_type() , boost::ref( relax ) , x , 0.0 , 10.0 , dt );
+
+    return 0;
+}

libs/numeric/odeint/examples/two_dimensional_phase_lattice.cpp

@@ -0,0 +1,158 @@
+/*
+ * two_dimensional_phase_lattice.cpp
+ *
+ * This example show how one can use matrices as state types in odeint.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2011-2013 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)
+ */
+
+#include <iostream>
+#include <map>
+#include <string>
+#include <fstream>
+
+#ifndef M_PI //not there on windows
+#define M_PI 3.1415927 //...
+#endif
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+//[ two_dimensional_phase_lattice_definition
+typedef boost::numeric::ublas::matrix< double > state_type;
+
+struct two_dimensional_phase_lattice
+{
+    two_dimensional_phase_lattice( double gamma = 0.5 )
+    : m_gamma( gamma ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
+    {
+        size_t size1 = x.size1() , size2 = x.size2();
+
+        for( size_t i=1 ; i<size1-1 ; ++i )
+        {
+            for( size_t j=1 ; j<size2-1 ; ++j )
+            {
+                dxdt( i , j ) =
+                        coupling_func( x( i + 1 , j ) - x( i , j ) ) +
+                        coupling_func( x( i - 1 , j ) - x( i , j ) ) +
+                        coupling_func( x( i , j + 1 ) - x( i , j ) ) +
+                        coupling_func( x( i , j - 1 ) - x( i , j ) );
+            }
+        }
+
+        for( size_t i=0 ; i<x.size1() ; ++i ) dxdt( i , 0 ) = dxdt( i , x.size2() -1 ) = 0.0;
+        for( size_t j=0 ; j<x.size2() ; ++j ) dxdt( 0 , j ) = dxdt( x.size1() -1 , j ) = 0.0;
+    }
+
+    double coupling_func( double x ) const
+    {
+        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
+    }
+
+    double m_gamma;
+};
+//]
+
+
+struct write_for_gnuplot
+{
+    size_t m_every , m_count;
+
+    write_for_gnuplot( size_t every = 10 )
+    : m_every( every ) , m_count( 0 ) { }
+
+    void operator()( const state_type &x , double t )
+    {
+        if( ( m_count % m_every ) == 0 )
+        {
+            clog << t << endl;
+            cout << "sp '-'" << endl;
+            for( size_t i=0 ; i<x.size1() ; ++i )
+            {
+                for( size_t j=0 ; j<x.size2() ; ++j )
+                {
+                    cout << i << "\t" << j << "\t" << sin( x( i , j ) ) << "\n";
+                }
+                cout << "\n";
+            }
+            cout << "e" << endl;
+        }
+
+        ++m_count;
+    }
+};
+
+class write_snapshots
+{
+public:
+
+    typedef std::map< size_t , std::string > map_type;
+
+    write_snapshots( void ) : m_count( 0 ) { }
+
+    void operator()( const state_type &x , double t )
+    {
+        map< size_t , string >::const_iterator it = m_snapshots.find( m_count );
+        if( it != m_snapshots.end() )
+        {
+            ofstream fout( it->second.c_str() );
+            for( size_t i=0 ; i<x.size1() ; ++i )
+            {
+                for( size_t j=0 ; j<x.size2() ; ++j )
+                {
+                    fout << i << "\t" << j << "\t" << x( i , j ) << "\t" << sin( x( i , j ) ) << "\n";
+                }
+                fout << "\n";
+            }
+        }
+        ++m_count;
+    }
+
+    map_type& snapshots( void ) { return m_snapshots; }
+    const map_type& snapshots( void ) const { return m_snapshots; }
+
+private:
+
+    size_t m_count;
+    map_type m_snapshots;
+};
+
+
+int main( int argc , char **argv )
+{
+    size_t size1 = 128 , size2 = 128;
+    state_type x( size1 , size2 , 0.0 );
+
+    for( size_t i=(size1/2-10) ; i<(size1/2+10) ; ++i )
+        for( size_t j=(size2/2-10) ; j<(size2/2+10) ; ++j )
+            x( i , j ) = static_cast<double>( rand() ) / RAND_MAX * 2.0 * M_PI;
+
+    write_snapshots snapshots;
+    snapshots.snapshots().insert( make_pair( size_t( 0 ) , string( "lat_0000.dat" ) ) );
+    snapshots.snapshots().insert( make_pair( size_t( 100 ) , string( "lat_0100.dat" ) ) );
+    snapshots.snapshots().insert( make_pair( size_t( 1000 ) , string( "lat_1000.dat" ) ) );
+    observer_collection< state_type , double > obs;
+    obs.observers().push_back( write_for_gnuplot( 10 ) );
+    obs.observers().push_back( snapshots );
+
+    cout << "set term x11" << endl;
+    cout << "set pm3d map" << endl;
+
+    integrate_const( runge_kutta4<state_type>() , two_dimensional_phase_lattice( 1.2 ) ,
+                     x , 0.0 , 1001.0 , 0.1 , boost::ref( obs ) );
+
+    // controlled steppers work only after ublas bugfix
+    //integrate_const( make_dense_output< runge_kutta_dopri5< state_type > >( 1E-6 , 1E-6 ) , two_dimensional_phase_lattice( 1.2 ) ,
+    //        x , 0.0 , 1001.0 , 0.1 , boost::ref( obs ) );
+
+
+    return 0;
+}

libs/numeric/odeint/examples/ublas/Jamfile.v2

@@ -0,0 +1,13 @@
+# Copyright 2011 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# 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)
+
+
+project
+    : requirements 
+      <define>BOOST_ALL_NO_LIB=1
+    ;
+    
+exe lorenz_ublas : lorenz_ublas.cpp ;

libs/numeric/odeint/examples/ublas/lorenz_ublas.cpp

@@ -0,0 +1,40 @@
+/*
+ * Copyright 2011-2013 Mario Mulansky
+ * Copyright 2012 Karsten Ahnert
+ *
+ * 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)
+ */
+
+
+#include <iostream>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
+typedef boost::numeric::ublas::vector< double > state_type;
+
+void lorenz( const state_type &x , state_type &dxdt , const 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] = -b * x[2] + x[0] * x[1];
+}
+
+using namespace boost::numeric::odeint;
+
+//[ublas_main
+int main()
+{
+    state_type x(3);
+    x[0] = 10.0; x[1] = 5.0 ; x[2] = 0.0;
+    typedef runge_kutta_dopri5< state_type > stepper;
+    integrate_const( make_dense_output< stepper >( 1E-6 , 1E-6 ) , lorenz , x ,
+                     0.0 , 10.0 , 0.1 );
+}
+//]

libs/numeric/odeint/examples/van_der_pol_stiff.cpp

@@ -0,0 +1,95 @@
+/*
+ * van_der_pol_stiff.cpp
+ *  
+ * Created on: Dec 12, 2011
+ *
+ * Copyright 2012 Karsten Ahnert
+ * Copyright 2012-2013 Rajeev Singh
+ * Copyright 2012-2013 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)
+ */
+
+#include <iostream>
+#include <fstream>
+#include <utility>
+
+#include <boost/numeric/odeint.hpp>
+
+#include <boost/phoenix/core.hpp>
+#include <boost/phoenix/operator.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+namespace phoenix = boost::phoenix;
+
+const double mu = 1000.0;
+
+
+typedef boost::numeric::ublas::vector< double > vector_type;
+typedef boost::numeric::ublas::matrix< double > matrix_type;
+
+struct vdp_stiff
+{
+    void operator()( const vector_type &x , vector_type &dxdt , double t )
+    {
+        dxdt[0] = x[1];
+        dxdt[1] = -x[0] - mu * x[1] * (x[0]*x[0]-1.0);
+    }
+};
+
+struct vdp_stiff_jacobi
+{
+    void operator()( const vector_type &x , matrix_type &J , const double &t , vector_type &dfdt )
+    {
+        J(0, 0) = 0.0;
+        J(0, 1) = 1.0;
+        J(1, 0) = -1.0 - 2.0*mu * x[0] * x[1];
+        J(1, 1) = -mu * ( x[0] * x[0] - 1.0);
+
+        dfdt[0] = 0.0;
+        dfdt[1] = 0.0;
+    }
+};
+
+
+int main( int argc , char **argv )
+{
+    //[ integrate_stiff_system
+    vector_type x( 2 );
+    /* initialize random seed: */
+    srand ( time(NULL) );
+
+    // initial conditions
+    for (int i=0; i<2; i++)
+        x[i] = 1.0; //(1.0 * rand()) / RAND_MAX;
+
+    size_t num_of_steps = integrate_const( make_dense_output< rosenbrock4< double > >( 1.0e-6 , 1.0e-6 ) ,
+            make_pair( vdp_stiff() , vdp_stiff_jacobi() ) ,
+            x , 0.0 , 1000.0 , 1.0
+            , cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << " " << phoenix::arg_names::arg1[1] << "\n"
+            );
+    //]
+    clog << num_of_steps << endl;
+
+
+
+    //[ integrate_stiff_system_alternative
+
+    vector_type x2( 2 );
+    // initial conditions
+    for (int i=0; i<2; i++)
+        x2[i] = 1.0; //(1.0 * rand()) / RAND_MAX;
+
+    //size_t num_of_steps2 = integrate_const( make_dense_output< runge_kutta_dopri5< vector_type > >( 1.0e-6 , 1.0e-6 ) ,
+    //        vdp_stiff() , x2 , 0.0 , 1000.0 , 1.0
+    //        , cout << phoenix::arg_names::arg2 << " " << phoenix::arg_names::arg1[0] << " " << phoenix::arg_names::arg1[1] << "\n"
+    //        );
+    //]
+    //clog << num_of_steps2 << endl;
+
+
+    return 0;
+}

libs/numeric/odeint/examples/vexcl/Jamfile.v2

@@ -0,0 +1,32 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2013 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 boost ;
+import os ;
+
+boost.use-project ;
+
+
+# change these lines to fit you configuration
+local HOME = [ os.environ HOME ] ;
+local VEXCL_INCLUDE = [ os.environ VEXCL_ROOT ] ;
+OPENCL_INCLUDE = /usr/local/cuda/include ;
+
+
+
+lib opencl : : <name>OpenCL ;
+
+project : requirements
+        <implicit-dependency>/boost//headers
+        <include>$(VEXCL_INCLUDE)
+        <include>$(OPENCL_INCLUDE)
+        <toolset>gcc:<cxxflags>-std=c++0x
+        <library>/boost//system/
+        ;
+
+exe lorenz_ensemble : lorenz_ensemble.cpp opencl ;

libs/numeric/odeint/examples/vexcl/lorenz_ensemble.cpp

@@ -0,0 +1,86 @@
+/*
+ * Copyright 2012 Karsten Ahnert
+ * Copyright 2013 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)
+ */
+
+#include <iostream>
+#include <vector>
+
+#include <vexcl/vexcl.hpp>
+
+#include <boost/numeric/odeint.hpp>
+//[ vexcl_includes
+#include <boost/numeric/odeint/external/vexcl/vexcl.hpp>
+//]
+
+namespace odeint = boost::numeric::odeint;
+
+//[ vexcl_state_types
+typedef vex::vector< double >    vector_type;
+typedef vex::multivector< double, 3 > state_type;
+//]
+
+
+//[ vexcl_system
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+struct sys_func
+{
+    const vector_type &R;
+
+    sys_func( const vector_type &_R ) : R( _R ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const
+    {
+        dxdt(0) = -sigma * ( x(0) - x(1) );
+        dxdt(1) = R * x(0) - x(1) - x(0) * x(2);
+        dxdt(2) = - b * x(2) + x(0) * x(1);
+    }
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+    using namespace std;
+    using namespace odeint;
+
+    //[ vexcl_main
+    // setup the opencl context
+    vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
+    std::cout << ctx << std::endl;
+
+    // set up number of system, time step and integration time
+    const size_t n = 1024 * 1024;
+    const double dt = 0.01;
+    const double t_max = 1000.0;
+
+    // initialize R
+    double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
+    std::vector<double> x( n * 3 ) , r( n );
+    for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
+    vector_type R( ctx.queue() , r );
+
+    // initialize the state of the lorenz ensemble
+    state_type X(ctx.queue(), n);
+    X(0) = 10.0;
+    X(1) = 10.0;
+    X(2) = 10.0;
+
+    // create a stepper
+    runge_kutta4< state_type > stepper;
+
+    // solve the system
+    integrate_const( stepper , sys_func( R ) , X , 0.0 , t_max , dt );
+    //]
+
+    std::vector< double > res( 3 * n );
+    vex::copy( X(0) , res );
+    for( size_t i=0 ; i<n ; ++i )
+        cout << r[i] << "\t" << res[i] << "\t" << "\n";
+}