manual merge

.travis.yml

@@ -1,5 +1,7 @@
 language: cpp
-
+os: 
+  - linux
+  - osx
 compiler:
   - gcc
   - clang
@@ -8,6 +10,12 @@ env:
   - CXXSTD=''
   - CXXSTD='cxxflags="-std=c++0x"'
 
+# For now disable gcc on osx as g++4.8 is not yet available
+matrix:
+  exclude:
+    - os: osx
+      compiler: gcc
+
 before_install:
  - wget http://sourceforge.net/projects/boost/files/boost/1.55.0/boost_1_55_0.tar.bz2/download -O /tmp/boost.tar.bz2
  - tar jxf /tmp/boost.tar.bz2
@@ -19,9 +27,11 @@ before_install:
  - cd $BOOST_ROOT
  - ./bootstrap.sh
  - cd $TRAVIS_BUILD_DIR
-
+ - if [ "$TRAVIS_OS_NAME" = "osx" ] && [ "$CC" = "gcc" ]; then export CC=gcc-4.8; fi
+ - $CC --version
 
 script:
- - $BOOST_ROOT/b2 toolset=$CC$GCCVER $CXXSTD
+ - $BOOST_ROOT/b2 toolset=$CC $CXXSTD
 # build in c++11 mode only with gcc
 # - if [ "$CXX" = "g++" ]; then $BOOST_ROOT/b2 -a toolset=$CC$GCCVER cxxflags="-std=c++0x"; fi
+

doc/Jamfile.v2

@@ -38,7 +38,7 @@ if --enable-index in  [ modules.peek : ARGV ]
         # PDF native index support is probably better for PDFs as then you actually get page numbers.
 
         <auto-index-script>odeint.idx # Specifies the name of the index script to load. 
-        <auto-index-prefix>../../../.. 
+        <auto-index-prefix>../include
         
         # Inform Quickbook that there is to be an index(es).
         <quickbook-define>enable_index

doc/details_state_types_algebras_operations.qbk

@@ -124,7 +124,7 @@ If we wouldn't specialize the `is_resizeable` template, the code would still
 compile but odeint would not adjust the size of temporary internal instances
 of my_vector and hence try to fill zero-sized vectors resulting in
 segmentation faults!
-The full example can be found in [github_link libs/numeric/odeint/examples/my_vector.cpp my_vector.cpp]
+The full example can be found in [github_link examples/my_vector.cpp my_vector.cpp]
 
 [endsect]
 
@@ -152,7 +152,7 @@ The following code shows the required template specializations:
 
 With these definitions odeint knows how to resize `std::list`s and so they can
 be used as state types.
-A complete example can be found in [github_link libs/numeric/odeint/examples/list_lattice.cpp list_lattice.cpp].
+A complete example can be found in [github_link examples/list_lattice.cpp list_lattice.cpp].
 
 [endsect]
 
@@ -407,7 +407,7 @@ The following code shows the corresponding definitions:
 Note again, that we haven't supported the requirements for controlled steppers,
 but only for simple Runge-Kutta methods.
 You can find the full example in [github_link
-libs/numeric/odeint/examples/ublas/lorenz_ublas.cpp lorenz_ublas.cpp].
+examples/ublas/lorenz_ublas.cpp lorenz_ublas.cpp].
 
 [endsect]
 /]
@@ -453,7 +453,7 @@ template argument list:
 [point3D_main]
 
 The whole example can be found in [github_link
-libs/numeric/odeint/examples/lorenz_point.cpp lorenz_point.cpp]
+examples/lorenz_point.cpp lorenz_point.cpp]
 
 [note For the most `state_types`, odeint is able to automatically determine
 the correct algebra and operations. But if you want to use your own `state_type`, as in this

doc/examples_table.qbk

@@ -15,128 +15,128 @@
 [table Examples Overview
   [[File] [Brief Description]]
 
-  [[[github_link libs/numeric/odeint/examples/bind_member_functions.cpp bind_member_functions.cpp]]
+  [[[github_link examples/bind_member_functions.cpp bind_member_functions.cpp]]
   [This examples shows how member functions can be used as system functions in odeint.]]
 
-  [[[github_link libs/numeric/odeint/examples/bind_member_functions.cpp bind_member_functions_cpp11.cpp]]
+  [[[github_link examples/bind_member_functions.cpp bind_member_functions_cpp11.cpp]]
   [This examples shows how member functions can be used as system functions in odeint with `std::bind` in C++11.]]
 
-  [[[github_link libs/numeric/odeint/examples/bulirsch_stoer.cpp bulirsch_stoer.cpp]]
+  [[[github_link examples/bulirsch_stoer.cpp bulirsch_stoer.cpp]]
   [Shows the usage of the Bulirsch-Stoer method.]]
 
-  [[[github_link libs/numeric/odeint/examples/chaotic_system.cpp chaotic_system.cpp]]
+  [[[github_link examples/chaotic_system.cpp chaotic_system.cpp]]
   [The chaotic system examples integrates the Lorenz system and calculates the Lyapunov exponents.]]
 
-  [[[github_link libs/numeric/odeint/examples/elliptic_functions.cpp elliptic_functions.cpp]]
+  [[[github_link examples/elliptic_functions.cpp elliptic_functions.cpp]]
   [Example calculating the elliptic functions using Bulirsch-Stoer and Runge-Kutta-Dopri5 Steppers with dense output.]] 
  
-  [[[github_link libs/numeric/odeint/examples/fpu.cpp fpu.cpp]]
+  [[[github_link examples/fpu.cpp fpu.cpp]]
   [The Fermi-Pasta-Ulam (FPU) example shows how odeint can be used to integrate lattice systems.]]
 
-  [[[github_link libs/numeric/odeint/examples/generation_functions.cpp generation_functions.cpp]]
+  [[[github_link examples/generation_functions.cpp generation_functions.cpp]]
   [Shows skeletal code on how to implement own factory functions.]]
 
-  [[[github_link libs/numeric/odeint/examples/harmonic_oscillator.cpp harmonic_oscillator.cpp]]
+  [[[github_link examples/harmonic_oscillator.cpp harmonic_oscillator.cpp]]
   [The harmonic oscillator examples gives a brief introduction to odeint and shows the usage of the classical Runge-Kutta-solvers.]]
 
-  [[[github_link libs/numeric/odeint/examples/harmonic_oscillator_units.cpp harmonic_oscillator_units.cpp]]
+  [[[github_link examples/harmonic_oscillator_units.cpp harmonic_oscillator_units.cpp]]
   [This examples shows how __boost_units can be used with odeint.]]
 
-  [[[github_link libs/numeric/odeint/examples/heun.cpp heun.cpp]]
+  [[[github_link examples/heun.cpp heun.cpp]]
   [The Heun example shows how an custom Runge-Kutta stepper can be created with odeint generic Runge-Kutta method.]]
 
-  [[[github_link libs/numeric/odeint/examples/list_lattice.cpp list_lattice.cpp]]
+  [[[github_link examples/list_lattice.cpp list_lattice.cpp]]
   [Example of a phase lattice integration using `std::list` as state type.]]
 
-  [[[github_link libs/numeric/odeint/examples/lorenz_point.cpp lorenz_point.cpp]]
+  [[[github_link examples/lorenz_point.cpp lorenz_point.cpp]]
   [Alternative way of integrating lorenz by using a self defined point3d data type as state type.]]
 
-  [[[github_link libs/numeric/odeint/examples/my_vector.cpp my_vector.cpp]]
+  [[[github_link examples/my_vector.cpp my_vector.cpp]]
   [Simple example showing how to get odeint to work with a self-defined vector type.]]
 
-  [[[github_link libs/numeric/odeint/examples/phase_oscillator_ensemble.cpp phase_oscillator_ensemble.cpp]]
+  [[[github_link examples/phase_oscillator_ensemble.cpp phase_oscillator_ensemble.cpp]]
   [The phase oscillator ensemble example shows how globally coupled oscillators can be analyzed and how statistical measures can be computed during integration.]]
 
-  [[[github_link libs/numeric/odeint/examples/resizing_lattice.cpp resizing_lattice.cpp]]
+  [[[github_link examples/resizing_lattice.cpp resizing_lattice.cpp]]
   [Shows the strength of odeint's memory management by simulating a Hamiltonian system on an expanding lattice.]]
 
-  [[[github_link libs/numeric/odeint/examples/simple1d.cpp simple1d.cpp]]
+  [[[github_link examples/simple1d.cpp simple1d.cpp]]
   [Integrating a simple, one-dimensional ODE showing the usage of integrate- and generate-functions.]]
 
-  [[[github_link libs/numeric/odeint/examples/solar_system.cpp solar_system.cpp]]
+  [[[github_link examples/solar_system.cpp solar_system.cpp]]
   [The solar system example shows the usage of the symplectic solvers.]]
 
-  [[[github_link libs/numeric/odeint/examples/stepper_details.cpp stepper_details.cpp]]
+  [[[github_link examples/stepper_details.cpp stepper_details.cpp]]
   [Trivial example showing the usability of the several stepper classes.]]
 
-  [[[github_link libs/numeric/odeint/examples/stiff_system.cpp stiff_system.cpp]]
+  [[[github_link examples/stiff_system.cpp stiff_system.cpp]]
   [The stiff system example shows the usage of the stiff solvers using the Jacobian of the system function.]]
 
-  [[[github_link libs/numeric/odeint/examples/stochastic_euler.cpp stochastic_euler.cpp]]
+  [[[github_link examples/stochastic_euler.cpp stochastic_euler.cpp]]
   [Implementation of a custom stepper - the stochastic euler -  for solving stochastic differential equations.]]
 
-  [[[github_link libs/numeric/odeint/examples/stuart_landau.cpp stuart_landau.cpp]]
+  [[[github_link examples/stuart_landau.cpp stuart_landau.cpp]]
   [The Stuart-Landau example shows how odeint can be used with complex state types.]]
 
-  [[[github_link libs/numeric/odeint/examples/two_dimensional_phase_lattice.cpp two_dimensional_phase_lattice.cpp]]
+  [[[github_link examples/two_dimensional_phase_lattice.cpp two_dimensional_phase_lattice.cpp]]
   [The 2D phase oscillator example shows how a two-dimensional lattice works with odeint and how matrix types can be used as state types in odeint.]]
 
-  [[[github_link libs/numeric/odeint/examples/van_der_pol_stiff.cpp van_der_pol_stiff.cpp]]
+  [[[github_link examples/van_der_pol_stiff.cpp van_der_pol_stiff.cpp]]
   [This stiff system example again shows the usage of the stiff solvers by integrating the van der Pol oscillator.]]
 
-  [[[github_link libs/numeric/odeint/examples/gmpxx/lorenz_gmpxx.cpp gmpxx/lorenz_gmpxx.cpp]]
+  [[[github_link examples/gmpxx/lorenz_gmpxx.cpp gmpxx/lorenz_gmpxx.cpp]]
   [This examples integrates the Lorenz system by means of an arbitrary precision type.]]
 
-  [[[github_link libs/numeric/odeint/examples/mtl/gauss_packet.cpp mtl/gauss_packet.cpp]]
+  [[[github_link examples/mtl/gauss_packet.cpp mtl/gauss_packet.cpp]]
   [The MTL-Gauss-packet example shows how the MTL can be easily used with odeint.]]
 
-  [[[github_link libs/numeric/odeint/examples/mtl/implicit_euler_mtl.cpp mtl/implicit_euler_mtl.cpp]]
+  [[[github_link examples/mtl/implicit_euler_mtl.cpp mtl/implicit_euler_mtl.cpp]]
   [This examples shows the usage of the MTL implicit Euler method with a sparse matrix type.]]
 
-  [[[github_link libs/numeric/odeint/examples/thrust/phase_oscillator_ensemble.cu thrust/phase_oscillator_ensemble.cu]]
+  [[[github_link examples/thrust/phase_oscillator_ensemble.cu thrust/phase_oscillator_ensemble.cu]]
   [The Thrust phase oscillator ensemble example shows how globally coupled oscillators can be analyzed with Thrust and CUDA, employing the power of modern graphic devices.]]
 
-  [[[github_link libs/numeric/odeint/examples/thrust/phase_oscillator_chain.cu thrust/phase_oscillator_chain.cu]]
+  [[[github_link examples/thrust/phase_oscillator_chain.cu thrust/phase_oscillator_chain.cu]]
   [The Thrust phase oscillator chain example shows how chains of nearest neighbor coupled oscillators can be integrated with Thrust and odeint.]]
 
-  [[[github_link libs/numeric/odeint/examples/thrust/lorenz_parameters.cu thrust/lorenz_parameters.cu]]
+  [[[github_link examples/thrust/lorenz_parameters.cu thrust/lorenz_parameters.cu]]
   [The Lorenz parameters examples show how ensembles of ordinary differential equations can be solved by means of Thrust to study
   the dependence of an ODE on some parameters.]]
 
-  [[[github_link libs/numeric/odeint/examples/thrust/relaxation.cu thrust/relaxation.cu]]
+  [[[github_link examples/thrust/relaxation.cu thrust/relaxation.cu]]
   [Another examples for the usage of Thrust.]]
 
-  [[[github_link libs/numeric/odeint/examples/ublas/lorenz_ublas.cpp ublas/lorenz_ublas.cpp]]
+  [[[github_link examples/ublas/lorenz_ublas.cpp ublas/lorenz_ublas.cpp]]
   [This example shows how the ublas vector types can be used with odeint.]]
 
-  [[[github_link libs/numeric/odeint/examples/vexcl/lorenz_ensemble.cpp vexcl/lorenz_ensemble.cpp]]
+  [[[github_link examples/vexcl/lorenz_ensemble.cpp vexcl/lorenz_ensemble.cpp]]
   [This example shows how the VexCL - a framework for OpenCL computation - can be used with odeint.]]
 
-  [[[github_link libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp openmp/lorenz_ensemble_simple.cpp]]
+  [[[github_link examples/openmp/lorenz_ensemble_simple.cpp openmp/lorenz_ensemble_simple.cpp]]
   [OpenMP Lorenz attractor parameter study with continuous data.]]
 
-  [[[github_link libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp openmp/lorenz_ensemble.cpp]]
+  [[[github_link examples/openmp/lorenz_ensemble.cpp openmp/lorenz_ensemble.cpp]]
   [OpenMP Lorenz attractor parameter study with split data.]]
 
-  [[[github_link libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp openmp/lorenz_ensemble_nested.cpp]]
+  [[[github_link examples/openmp/lorenz_ensemble.cpp openmp/lorenz_ensemble_nested.cpp]]
   [OpenMP Lorenz attractor parameter study with nested `vector_space_algebra`.]]
 
-  [[[github_link libs/numeric/odeint/examples/openmp/phase_chain.cpp openmp/phase_chain.cpp]]
+  [[[github_link examples/openmp/phase_chain.cpp openmp/phase_chain.cpp]]
   [OpenMP nearest neighbour coupled phase chain with continuous state.]]
 
-  [[[github_link libs/numeric/odeint/examples/openmp/phase_chain_omp_state.cpp openmp/phase_chain_omp_state.cpp]]
+  [[[github_link examples/openmp/phase_chain_omp_state.cpp openmp/phase_chain_omp_state.cpp]]
   [OpenMP nearest neighbour coupled phase chain with split state.]]
 
-  [[[github_link libs/numeric/odeint/examples/mpi/phase_chain.cpp mpi/phase_chain.cpp]]
+  [[[github_link examples/mpi/phase_chain.cpp mpi/phase_chain.cpp]]
   [MPI nearest neighbour coupled phase chain.]]
 
-  [[[github_link libs/numeric/odeint/examples/2d_lattice/spreading.cpp 2d_lattice/spreading.cpp]]
+  [[[github_link examples/2d_lattice/spreading.cpp 2d_lattice/spreading.cpp]]
   [This examples shows how a `vector< vector< T > >` can be used a state type for odeint and how a resizing mechanism of this state can be implemented.]]
 
-  [[[github_link libs/numeric/odeint/examples/quadmath/black_hole.cpp quadmath/black_hole.cpp]]
+  [[[github_link examples/quadmath/black_hole.cpp quadmath/black_hole.cpp]]
   [This examples shows how gcc libquadmath can be used with odeint. It provides a high precision floating point type which is adapted to odeint in this example.]]
   
-  [[[github_link libs/numeric/odeint/examples/molecular_dynamics.cpp molecular_dynamics.cpp]]
+  [[[github_link examples/molecular_dynamics.cpp molecular_dynamics.cpp]]
   [A very basic molecular dynamics simulation with the Velocity-Verlet method.]]
 
 ]

doc/getting_started.qbk

@@ -133,7 +133,7 @@ to pass this container to the integration function:
 
 That is all. You can use functional libraries like __boost_lambda or __boost_phoenix to ease the creation of observer functions.
 
-The full cpp file for this example can be found here: [github_link libs/numeric/odeint/examples/harmonic_oscillator.cpp harmonic_oscillator.cpp]
+The full cpp file for this example can be found here: [github_link examples/harmonic_oscillator.cpp harmonic_oscillator.cpp]
 
 
 

doc/odeint.qbk

@@ -181,7 +181,7 @@
 [template sub[x]'''<subscript>'''[x]'''</subscript>''']
 [template subl[x]'''<subscript>'''__space[x]'''</subscript>''']
 
-[template github_link[url text]'''<ulink url="https://github.com/headmyshoulder/odeint-v2/tree/master/'''[url]'''" target="_blank">'''[text]'''</ulink>''']
+[template github_link[url text]'''<ulink url="https://github.com/headmyshoulder/odeint-v2/blob/master/'''[url]'''" target="_blank">'''[text]'''</ulink>''']
 [/ [template github_link[url text]'''<ulink url="../../../../../'''[url]'''" target="_blank">'''[text]'''</ulink>''']]
 
 

doc/tutorial_chaotic_system.qbk

@@ -157,6 +157,6 @@ Having integrated a sufficient number of transients steps we are now able to cal
 
 [lyapunov_full_code]
 
-The full code can be found here:  [github_link libs/numeric/odeint/examples/chaotic_system.cpp chaotic_system.cpp]
+The full code can be found here:  [github_link examples/chaotic_system.cpp chaotic_system.cpp]
 
 [endsect]

doc/tutorial_harmonic_oscillator.qbk

@@ -117,7 +117,7 @@ odeint supports iterators for solving ODEs. That is, you instantiate a pair of i
 
 [endsect]
 
-The full source file for this example can be found here: [github_link libs/numeric/odeint/examples/harmonic_oscillator.cpp harmonic_oscillator.cpp]
+The full source file for this example can be found here: [github_link examples/harmonic_oscillator.cpp harmonic_oscillator.cpp]
 
 
 [endsect]

doc/tutorial_parallel.qbk

@@ -70,7 +70,7 @@ Note, that you can specify the OpenMP scheduling by calling `omp_set_schedule`
 in the beginning of your program:
 [phase_chain_scheduling]
 
-See [github_link libs/numeric/odeint/examples/openmp/phase_chain.cpp
+See [github_link examples/openmp/phase_chain.cpp
 openmp/phase_chain.cpp] for the complete example.
 
 [heading Split state]
@@ -109,7 +109,7 @@ back together into a single vector.
 and supports any model of Random Access Range as the outer, parallel state type,
 and will use the given algebra on its elements.]
 
-See [github_link libs/numeric/odeint/examples/openmp/phase_chain_omp_state.cpp
+See [github_link examples/openmp/phase_chain_omp_state.cpp
 openmp/phase_chain_omp_state.cpp] for the complete example.
 
 [endsect]
@@ -165,7 +165,7 @@ guarded by any barriers either, so if you don't manually place any (for example
 in parameter studies cases where the elements are completely independent) you
 might see the nodes diverging, returning from this call at different times.]
 
-See [github_link libs/numeric/odeint/examples/mpi/phase_chain.cpp
+See [github_link examples/mpi/phase_chain.cpp
 mpi/phase_chain.cpp] for the complete example.
 
 [endsect]

doc/tutorial_solar_system.qbk

@@ -102,7 +102,7 @@ These integration routine was used to produce the above sketch of the solar syst
 
 [tip You can use C++11 lambda to create the observers]
 
-The full example can be found here: [github_link libs/numeric/odeint/examples/solar_system.cpp solar_system.cpp]
+The full example can be found here: [github_link examples/solar_system.cpp solar_system.cpp]
 
 [endsect]
 

doc/tutorial_special_topics.qbk

@@ -46,7 +46,7 @@ found in the section __adapt_state_types.
 [stuart_landau_integration]
 
 The full cpp file for the Stuart-Landau example can be found here [github_link
-libs/numeric/odeint/examples/stuart_landau.cpp stuart_landau.cpp]
+examples/stuart_landau.cpp stuart_landau.cpp]
 
 [endsect]
 
@@ -75,7 +75,7 @@ The observer uses a reference to the system object to calculate the local energi
 
 [fpu_observer]
 
-The full cpp file for this FPU example can be found here [github_link libs/numeric/odeint/examples/fpu.cpp fpu.cpp]
+The full cpp file for this FPU example can be found here [github_link examples/fpu.cpp fpu.cpp]
 
 [endsect]
 
@@ -103,7 +103,7 @@ Now, we do several integrations for different values of ['__epsilon] and record
 
 [phase_oscillator_ensemble_integration]
 
-The full cpp file for this example can be found here [github_link libs/numeric/odeint/examples/phase_oscillator_ensemble.cpp phase_oscillator_ensemble.cpp]
+The full cpp file for this example can be found here [github_link examples/phase_oscillator_ensemble.cpp phase_oscillator_ensemble.cpp]
 
 [endsect]
 
@@ -159,7 +159,7 @@ It is quite easy but the compilation time might take very long. Furthermore, the
 memory consuming. For example the unit test for the usage of __boost_units in odeint take up to 4 GB
 of memory at compilation.]
 
-The full cpp file for this example can be found here [github_link libs/numeric/odeint/examples/harmonic_oscillator_units.cpp harmonic_oscillator_units.cpp].
+The full cpp file for this example can be found here [github_link examples/harmonic_oscillator_units.cpp harmonic_oscillator_units.cpp].
 
 [endsect] 
 
@@ -179,7 +179,7 @@ In principle this is all. Please note, that the above code is far from being opt
 
 [$phase_lattice_2d_0000.jpg] [$phase_lattice_2d_0100.jpg] [$phase_lattice_2d_1000.jpg]
 
-The full cpp for this example can be found here [github_link libs/numeric/odeint/examples/two_dimensional_phase_lattice.cpp two_dimensional_phase_lattice.cpp].
+The full cpp for this example can be found here [github_link examples/two_dimensional_phase_lattice.cpp two_dimensional_phase_lattice.cpp].
 
 [endsect]
 
@@ -224,13 +224,13 @@ The actual integration then is straight forward:
 
 [mp_lorenz_int]
 
-The full example can be found at [github_link libs/numeric/odeint/examples/multiprecision/lorenz_mp.cpp lorenz_mp.cpp].
+The full example can be found at [github_link examples/multiprecision/lorenz_mp.cpp lorenz_mp.cpp].
 Another example that compares the accuracy of the high precision type with
-standard double can be found at [github_link libs/numeric/odeint/examples/multiprecision/cmp_precision.cpp cmp_precision.cpp].
+standard double can be found at [github_link examples/multiprecision/cmp_precision.cpp cmp_precision.cpp].
 
 Furthermore, odeint can also be run with other multiprecision libraries,
 e.g. [@http://gmplib.org/ gmp].
-An example for this is given in [github_link libs/numeric/odeint/examples/gmpxx/lorenz_gmpxx.cpp lorenz_gmpxx.cpp].
+An example for this is given in [github_link examples/gmpxx/lorenz_gmpxx.cpp lorenz_gmpxx.cpp].
 
 [endsect]
 
@@ -267,7 +267,7 @@ The `do_resize` function simply calls `vector.resize` of `q` , `p` and `distr`.
 
 [resizing_lattice_resize_function]
 
-The full example can be found in [github_link libs/numeric/odeint/examples/resizing_lattice.cpp resizing_lattice.cpp]
+The full example can be found in [github_link examples/resizing_lattice.cpp resizing_lattice.cpp]
 
 [endsect]
 

doc/tutorial_stiff_systems.qbk

@@ -56,7 +56,7 @@ During the integration 71 steps have been done. Comparing to a classical Runge-K
 
 Note, that we have used __boost_phoenix, a great functional programming library, to create and compose the observer.
 
-The full example can be found here: [github_link libs/numeric/odeint/examples/stiff_system.cpp stiff_system.cpp]
+The full example can be found here: [github_link examples/stiff_system.cpp stiff_system.cpp]
 
 
 [endsect]

doc/tutorial_thrust_cuda.qbk

@@ -107,7 +107,7 @@ Then, it is straightforward to integrate the phase ensemble by creating an insta
 
 We have to use `boost::ref` here in order to pass the rhs class as reference and not by value. This ensures that the natural frequencies of each oscillator are not copied when calling `integrate_const`. In the full example the performance and results of the Runge-Kutta-4 and the Dopri5 solver are compared.
 
-The full example can be found at [github_link libs/numeric/odeint/examples/thrust/phase_oscillator_ensemble.cu phase_oscillator_example.cu].
+The full example can be found at [github_link examples/thrust/phase_oscillator_ensemble.cu phase_oscillator_example.cu].
 
 [endsect]
 
@@ -130,7 +130,7 @@ Now we put everything together. We create random initial conditions and decreasi
 
 [thrust_phase_chain_integration]
 
-The full example can be found at [github_link libs/numeric/odeint/examples/thrust/phase_oscillator_chain.cu phase_oscillator_chain.cu].
+The full example can be found at [github_link examples/thrust/phase_oscillator_chain.cu phase_oscillator_chain.cu].
 
 [endsect]
 
@@ -153,7 +153,7 @@ The next thing we have to implement is the Lorenz system without perturbations.
 
 As `state_type` a `thrust::device_vector` or a __boost_range of a `device_vector` is used. The length of the state is ['3N] where ['N] is the number of systems. The system is encoded into this vector such that all ['x] components come first, then every ['y] components and finally every ['z] components. Implementing the device function is then a simple task, you only have to decompose the tuple originating from the zip iterators.
 
-Besides the system without perturbations we furthermore need to calculate the system including linearized equations governing the time evolution of small perturbations. Using the method from above this is straightforward, with a small difficulty that Thrust's tuples have a maximal arity of 10. But this is only a small problem since we can create a zip iterator packed with zip iterators. So the top level zip iterator contains one zip iterator for the state, one normal iterator for the parameter, and one zip iterator for the derivative. Accessing the elements of this tuple in the system function is then straightforward, you unpack the tuple with `thrust::get<>()`. We will not show the code here, it is to large. It can be found [github_link libs/numeric/odeint/examples/thrust/lorenz_parameters.cu here] and is easy to understand.
+Besides the system without perturbations we furthermore need to calculate the system including linearized equations governing the time evolution of small perturbations. Using the method from above this is straightforward, with a small difficulty that Thrust's tuples have a maximal arity of 10. But this is only a small problem since we can create a zip iterator packed with zip iterators. So the top level zip iterator contains one zip iterator for the state, one normal iterator for the parameter, and one zip iterator for the derivative. Accessing the elements of this tuple in the system function is then straightforward, you unpack the tuple with `thrust::get<>()`. We will not show the code here, it is to large. It can be found [github_link examples/thrust/lorenz_parameters.cu here] and is easy to understand.
 
 Furthermore, we need an observer which determines the norm of the perturbations, normalizes them and averages the logarithm of the norm. The device functor which is used within this observer is defined
 
@@ -165,7 +165,7 @@ Now we complete the whole code to calculate the Lyapunov exponents. First, we ha
 
 [thrust_lorenz_parameters_integration]
 
-The full example can be found at [github_link libs/numeric/odeint/examples/thrust/lorenz_parameters.cu lorenz_parameters.cu].
+The full example can be found at [github_link examples/thrust/lorenz_parameters.cu lorenz_parameters.cu].
 
 [endsect]
 

examples/Jamfile.v2

@@ -8,7 +8,6 @@
 
 project
     : requirements 
-      <include>../../../..
       <define>BOOST_ALL_NO_LIB=1
     : 
     ;

examples/molecular_dynamics_cells.cpp

@@ -22,6 +22,7 @@
 #include <algorithm>
 #include <tuple>
 #include <iostream>
+#include <random>
 
 #include <boost/range/algorithm/for_each.hpp>
 #include <boost/range/algorithm/sort.hpp>

examples/mpi/Jamfile.v2

@@ -7,7 +7,6 @@
 
 project
     : requirements
-      <include>../../../../..
       <define>BOOST_ALL_NO_LIB=1
       <library>/boost//mpi
       <library>/boost//timer

examples/mtl/Jamfile.v2

@@ -9,7 +9,6 @@ MTL4_INCLUDE = /home/mario/MTL4 ;
 
 project
     : requirements 
-      <include>../../../../..
       <include>$(MTL4_INCLUDE)
       <define>BOOST_ALL_NO_LIB=1
     ;

examples/multiprecision/Jamfile

@@ -7,7 +7,6 @@
 
 project
     : requirements 
-      <include>../../../../..
       <define>BOOST_ALL_NO_LIB=1
     : 
     ;

examples/nt2/Jamfile.v2

@@ -24,7 +24,6 @@ local NT2_SIMD_FLAGS = [ os.environ NT2_SIMD_FLAGS ] ;
 project
   : requirements
     <define>BOOST_ALL_NO_LIB=1
-    <include>../../../../..
     <include>$(NT2_ROOT_PATH)/include/
     <link>static
     <toolset>gcc:<cxxflags>-DBOOST_SIMD_NO_STRICT_ALIASING

examples/openmp/Jamfile.v2

@@ -10,7 +10,6 @@ import openmp : * ;
 
 project
     : requirements
-      <include>../../../../..
       <include>..
       <define>BOOST_ALL_NO_LIB=1
       <library>/boost//timer

examples/quadmath/Jamfile.v2

@@ -7,7 +7,6 @@
 
 project
     : requirements 
-      <include>../../../..
       <define>BOOST_ALL_NO_LIB=1
     : 
     ;

examples/stochastic_euler.cpp

@@ -109,10 +109,10 @@ struct ornstein_det
 
 struct ornstein_stoch
 {
-    boost::mt19937 m_rng;
+    boost::mt19937 &m_rng;
     boost::normal_distribution<> m_dist;
 
-    ornstein_stoch( double sigma ) : m_rng() , m_dist( 0.0 , sigma ) { }
+  ornstein_stoch( boost::mt19937 &rng , double sigma ) : m_rng( rng ) , m_dist( 0.0 , sigma ) { }
 
     void operator()( const state_type &x , state_type &dxdt )
     {
@@ -137,9 +137,10 @@ int main( int argc , char **argv )
     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( 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;

examples/thrust/Makefile

@@ -1,66 +1,34 @@
-# Copyright 2011-2013 Mario Mulansky
+# 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
 
-
-# CUDA_ROOT = /home/karsten/boost/cuda4.1/cuda/
+# path to the cuda installation
 CUDA_ROOT = /usr/local/cuda
+# target architecture
 ARCH = sm_13
 
-CC = gcc
-CXX = g++
 NVCC = $(CUDA_ROOT)/bin/nvcc
-# NVCC = g++
 
-INCLUDES += -I$(BOOST_ROOT) -I$(THRUST_ROOT) -I$(CUDA_ROOT)/include -I../../../../..
+INCLUDES += -I../../include/ -I$(BOOST_ROOT)
 
-NVCCFLAGS = -O3 $(INCLUDES) -arch $(ARCH) --compiler-bindir=/opt/gcc4.6.2/bin/  -Xcompiler -fopenmp -DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP
-# NVCCFLAGS = -O3 $(INCLUDES) -arch $(ARCH) --compiler-bindir=/usr/bin/g++-4.3 -Xcompiler -fopenmp -DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP 
+NVCCFLAGS = -O3 $(INCLUDES) -arch $(ARCH)
 
-# NVCCFLAGS = -O3 $(INCLUDES) -arch $(ARCH) --compiler-bindir=/usr/bin/g++-4.3
-
-#--compiler-bindir=/usr/bin/g++-4.4 
-#-Xcompiler -fopenmp -DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP 
-# NVCCFLAGS = -O3 $(INCLUDES) -arch $(ARCH) --compiler-bindir=/usr/bin/g++-4.3 
-#-Xcompiler -fopenmp -DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP 
-
-
-LDLIBS = -lstdc++ -lm -lcudart -lgomp
-LDFLAGS = -L$(CUDA_ROOT)/lib64
-
-%.co : %.cu
+%.o : %.cu
 	$(NVCC) $(NVCCFLAGS) -c $< -o $@
 
+% : %.o
+	$(NVCC) $(NVCCFLAGS) -o $@ $<
+
 
 all : phase_oscillator_chain phase_oscillator_ensemble lorenz_parameters relaxation
 
 
-phase_oscillator_chain.co : phase_oscillator_chain.cu
-phase_oscillator_chain : phase_oscillator_chain.co
-	$(CC) -o phase_oscillator_chain $(LDFLAGS) $(LDLIBS) phase_oscillator_chain.co
-
-
-
-phase_oscillator_ensemble.co : phase_oscillator_ensemble.cu
-phase_oscillator_ensemble : phase_oscillator_ensemble.co
-	$(CC) -o phase_oscillator_ensemble $(LDFLAGS) $(LDLIBS) phase_oscillator_ensemble.co
-
-
-
-lorenz_parameters : lorenz_parameters.co
-	$(CC) -o lorenz_parameters $(LDFLAGS) $(LDLIBS) lorenz_parameters.co
-lorenz_parameters.co : lorenz_parameters.cu
-
-
-relaxation : relaxation.co
-	$(CC) -o relaxation $(LDFLAGS) $(LDLIBS) relaxation.co
-relaxation.co : relaxation.cu
-
-
-
 clean :
-	-rm *~ *.o *.co phase_oscillator_chain phase_oscillator_ensemble lorenz_parameters relaxation 
+	-rm *~ *.o phase_oscillator_chain phase_oscillator_ensemble lorenz_parameters relaxation 

examples/ublas/Jamfile.v2

@@ -7,7 +7,6 @@
 
 project
     : requirements 
-      <include>../../../../..
       <define>BOOST_ALL_NO_LIB=1
     ;
     

examples/vexcl/Jamfile.v2

@@ -23,7 +23,6 @@ lib opencl : : <name>OpenCL ;
 
 project : requirements
         <library>/boost//headers
-        <include>../../../../..
         <include>$(VEXCL_INCLUDE)
         <include>$(OPENCL_INCLUDE)
         <toolset>gcc:<cxxflags>-std=c++0x

include/boost/numeric/odeint/algebra/algebra_dispatcher.hpp

@@ -17,8 +17,9 @@
 #ifndef BOOST_NUMERIC_ODEINT_ALGEBRA_ALGEBRA_DISPATCHER_HPP_INCLUDED
 #define BOOST_NUMERIC_ODEINT_ALGEBRA_ALGEBRA_DISPATCHER_HPP_INCLUDED
 
-#include <complex>
+#include <boost/numeric/odeint/config.hpp>
 
+#include <complex>
 #include <boost/type_traits/is_floating_point.hpp>
 
 #include <boost/numeric/ublas/vector.hpp>
@@ -38,14 +39,14 @@ namespace odeint {
 template< class StateType , class Enabler = void >
 struct algebra_dispatcher_sfinae
 {
-    //  range_algebra is the standard algebra^
+    // range_algebra is the standard algebra
     typedef range_algebra algebra_type;
 };
 
 template< class StateType >
 struct algebra_dispatcher : algebra_dispatcher_sfinae< StateType > { };
 
-//specialize for array
+// specialize for array
 template< class T , size_t N >
 struct algebra_dispatcher< boost::array< T , N > >
 {
@@ -84,4 +85,26 @@ struct algebra_dispatcher< boost::numeric::ublas::matrix< T , L , A > >
 }
 }
 
+#ifdef BOOST_NUMERIC_ODEINT_CXX11
+
+// c++11 mode: specialization for std::array if available
+
+#include <array>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+    
+// specialize for std::array
+template< class T , size_t N >
+struct algebra_dispatcher< std::array< T , N > >
+{
+    typedef array_algebra algebra_type;
+};
+
+} } }
+
+#endif
+
+
 #endif

include/boost/numeric/odeint/algebra/array_algebra.hpp

@@ -3,7 +3,11 @@
  boost/numeric/odeint/algebra/array_algebra.hpp
 
  [begin_description]
- Algebra for boost::array. Highly specialized for odeint. Const arguments are introduce to work with odeint.
+ Algebra for Arrays. Highly specialized for odeint. Const arguments are
+ introduce to work with odeint.
+ The Array algebra can be used for Array structures with two template
+ parameters: 
+ Array<T, N>
  [end_description]
 
  Copyright 2011-2013 Mario Mulansky
@@ -29,234 +33,251 @@ namespace odeint {
 
 struct array_algebra
 {
-    template< typename T , size_t dim , class Op >
-    static void for_each1( boost::array< T , dim > &s1 , Op op )
+    //template< typename T , size_t dim , class Op >
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each1( Array< T, dim > &s1, Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] );
     }
 
-
-    template< typename T1 , typename T2 , size_t dim , class Op >
-    static void for_each2( boost::array< T1 , dim > &s1 ,
-            const boost::array< T2 , dim > &s2 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each2( Array< T, dim > &s1, const Array< T, dim > &s2,
+                           Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each3( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each3( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] );
     }
 
     /* different const signature - required for the scale_sum_swap2 operation */
-    template< typename T , size_t dim , class Op >
-    static void for_each3( boost::array< T , dim > &s1 ,
-            boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each3( Array< T , dim > &s1 ,
+            Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each4( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each4( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each5( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each5( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each6( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each6( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each7( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each7( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each8( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each8( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each9( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each9( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each10( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each10( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each11( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each11( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each12( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each12( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each13( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 ,
-            const boost::array< T , dim > &s13 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each13( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 ,
+            const Array< T , dim > &s13 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each14( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 ,
-            const boost::array< T , dim > &s13 ,
-            const boost::array< T , dim > &s14 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each14( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 ,
+            const Array< T , dim > &s13 ,
+            const Array< T , dim > &s14 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] , s14[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each15( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 ,
-            const boost::array< T , dim > &s13 ,
-            const boost::array< T , dim > &s14 ,
-            const boost::array< T , dim > &s15 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each15( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 ,
+            const Array< T , dim > &s13 ,
+            const Array< T , dim > &s14 ,
+            const Array< T , dim > &s15 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] , s14[i] , s15[i] );
     }
 
 
-    template< typename T , size_t dim >
-    static typename norm_result_type< boost::array< T , dim > >::type norm_inf( const boost::array< T , dim > &s )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim>
+    static typename norm_result_type< Array< T , dim > >::type norm_inf( const Array< T , dim > &s )
     {
         BOOST_USING_STD_MAX();
         using std::abs;
-        typedef typename norm_result_type< boost::array< T , dim > >::type result_type;
+        typedef typename norm_result_type< Array< T , dim > >::type result_type;
         result_type init = static_cast< result_type >( 0 );
         for( size_t i=0 ; i<dim ; ++i )
             init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type >(abs(s[i])) );

include/boost/numeric/odeint/algebra/multi_array_algebra.hpp

@@ -113,13 +113,13 @@ struct multi_array_algebra
     template< class S1 , class S2 , class S3 , class S4 , class S5 , class S6 ,class S7 , class S8 , class S9 , class S10 , class S11 , class S12 , class S13 , class S14 , class Op >
     static void for_each14( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , S5 &s5 , S6 &s6 , S7 &s7 , S8 &s8 , S9 &s9 , S10 &s10 , S11 &s11 , S12 &s12 , S13 &s13 , S14 &s14 , Op op )
     {
-        detail::for_each14( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s14.data() , op );
+        detail::for_each14( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , op );
     }
 
     template< class S1 , class S2 , class S3 , class S4 , class S5 , class S6 ,class S7 , class S8 , class S9 , class S10 , class S11 , class S12 , class S13 , class S14 , class S15 , class Op >
     static void for_each15( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , S5 &s5 , S6 &s6 , S7 &s7 , S8 &s8 , S9 &s9 , S10 &s10 , S11 &s11 , S12 &s12 , S13 &s13 , S14 &s14 , S15 &s15 , Op op )
     {
-        detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s14.data() , s15.data() , op );
+        detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , s15.data() , op );
     }
     
     template< typename S >

include/boost/numeric/odeint/external/eigen/eigen.hpp

@@ -0,0 +1,27 @@
+/*
+  [auto_generated]
+  boost/numeric/odeint/external/eigen/eigen.hpp
+
+  [begin_description]
+  tba.
+  [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)
+*/
+
+
+#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED
+
+
+#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen_resize.hpp>
+
+
+#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED

include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp

@@ -71,6 +71,19 @@ operator/(const Eigen::MatrixBase<D1> &x1, const Eigen::MatrixBase<D2> &x2) {
     return x1.cwiseQuotient(x2);
 }
 
+
+template< typename D >
+inline const 
+typename Eigen::CwiseUnaryOp<
+    typename Eigen::internal::scalar_abs_op<
+        typename Eigen::internal::traits< D >::Scalar > ,
+        const D >
+abs( const Eigen::MatrixBase< D > &m ) {
+    return m.cwiseAbs();
+}
+
+
+
 } // end Eigen namespace
 
 

include/boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp

@@ -0,0 +1,49 @@
+/*
+  [auto_generated]
+  boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp
+
+  [begin_description]
+  tba.
+  [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)
+*/
+
+
+#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+template< class Derived >
+struct algebra_dispatcher_sfinae< Derived ,
+                      typename boost::enable_if< typename boost::is_base_of< Eigen::MatrixBase< Derived > , Derived >::type >::type >
+{ 
+    typedef vector_space_algebra algebra_type;
+};
+
+
+template < class Derived  >
+struct algebra_dispatcher_sfinae< Derived ,
+                      typename boost::enable_if< typename boost::is_base_of< Eigen::ArrayBase< Derived > , Derived >::type >::type >
+{ 
+    typedef vector_space_algebra algebra_type;
+};
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED

include/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp

@@ -24,6 +24,7 @@
 #include <boost/numeric/odeint/external/thrust/thrust_algebra.hpp>
 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
 
+// specializations for the standard thrust containers
 
 namespace boost {
 namespace numeric {
@@ -48,5 +49,59 @@ struct algebra_dispatcher< thrust::device_vector< T , A > >
 } // namespace boost
 
 
+// add support for thrust backend vectors, if available
+
+#include <thrust/version.h>
+
+#if THRUST_VERSION >= 100600
+
+// specialization for thrust cpp vector
+#include <thrust/system/cpp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::cpp::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+
+// specialization for thrust omp vector
+#ifdef _OPENMP
+#include <thrust/system/omp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::omp::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+#endif // _OPENMP
+
+// specialization for thrust tbb vector
+#ifdef TBB_VERSION_MAJOR
+#include <thrust/system/tbb/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::tbb::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+#endif // TBB_VERSION_MAJOR
+
+// specialization for thrust cuda vector
+#ifdef __CUDACC__
+#include <thrust/system/cuda/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::cuda::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+#endif // __CUDACC__
+
+#endif // THRUST_VERSION >= 100600
+
 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_ALGEBRA_DISPATCHER_HPP_DEFINED
 

include/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp

@@ -6,8 +6,8 @@
   operations_dispatcher specialization for thrust
   [end_description]
 
-  Copyright 2013 Karsten Ahnert
-  Copyright 2013 Mario Mulansky
+  Copyright 2013-2014 Karsten Ahnert
+  Copyright 2013-2014 Mario Mulansky
 
   Distributed under the Boost Software License, Version 1.0.
   (See accompanying file LICENSE_1_0.txt or
@@ -24,6 +24,7 @@
 #include <boost/numeric/odeint/external/thrust/thrust_operations.hpp>
 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
 
+// support for the standard thrust containers
 
 namespace boost {
 namespace numeric {
@@ -47,6 +48,60 @@ struct operations_dispatcher< thrust::device_vector< T , A > >
 } // namespace numeric
 } // namespace boost
 
+// add support for thrust backend vectors, if available
 
-#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_ALGEBRA_DISPATCHER_HPP_DEFINED
+#include <thrust/version.h>
+
+#if THRUST_VERSION >= 100600
+
+// specialization for thrust cpp vector
+#include <thrust/system/cpp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::cpp::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+
+// specialization for thrust omp vector
+#ifdef _OPENMP
+#include <thrust/system/omp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::omp::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+#endif // _OPENMP
+
+// specialization for thrust tbb vector
+#ifdef TBB_VERSION_MAJOR
+#include <thrust/system/tbb/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::tbb::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+#endif // TBB_VERSION_MAJOR
+
+// specialization for thrust cuda vector
+#ifdef __CUDACC__
+#include <thrust/system/cuda/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::cuda::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+#endif // __CUDACC__
+
+#endif // THRUST_VERSION >= 100600
+
+
+#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_OPERATIONS_DISPATCHER_HPP_DEFINED
 

include/boost/numeric/odeint/external/thrust/thrust_resize.hpp

@@ -6,7 +6,7 @@
  Enable resizing for thrusts device and host_vector.
  [end_description]
 
- Copyright 2010-2012 Mario Mulansky
+ Copyright 2010-2014 Mario Mulansky
  Copyright 2010-2011 Karsten Ahnert
 
  Distributed under the Boost Software License, Version 1.0.
@@ -18,102 +18,168 @@
 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED
 #define BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED
 
+#include <boost/range.hpp>
 
 #include <thrust/device_vector.h>
 #include <thrust/host_vector.h>
+#include <thrust/distance.h>
 
+#include <boost/numeric/odeint/util/resize.hpp>
+#include <boost/numeric/odeint/util/same_size.hpp>
 #include <boost/numeric/odeint/util/copy.hpp>
 
 namespace boost {
 namespace numeric {
 namespace odeint {
 
-template< class T >
-struct is_resizeable< thrust::device_vector< T > >
-{
-    struct type : public boost::true_type { };
-    const static bool value = type::value;
-};
+// some macros that define the necessary utilities
 
-template< class T >
-struct same_size_impl< thrust::device_vector< T > , thrust::device_vector< T > >
-{
-    static bool same_size( const thrust::device_vector< T > &x , const thrust::device_vector< T > &y )
-    {
-        return x.size() == y.size();
-    }
-};
+#define ODEINT_THRUST_VECTOR_IS_RESIZEABLE( THRUST_VECTOR ) \
+template< class T , class A >                               \
+struct is_resizeable< THRUST_VECTOR<T,A> >                  \
+{                                                           \
+    struct type : public boost::true_type { };              \
+    const static bool value = type::value;                  \
+};                                                          \
 
-template< class T >
-struct resize_impl< thrust::device_vector< T > , thrust::device_vector< T > >
-{
-    static void resize( thrust::device_vector< T > &x , const thrust::device_vector< T > &y )
-    {
-        x.resize( y.size() );
-    }
-};
+#define ODEINT_TRHUST_VECTOR_RESIZE_IMPL( THRUST_VECTOR )     \
+template< class T, class A >                                  \
+struct resize_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> > \
+{                                                             \
+    static void resize( THRUST_VECTOR<T,A> &x ,               \
+                        const THRUST_VECTOR<T,A> &y )         \
+    {                                                         \
+        x.resize( y.size() );                                 \
+    }                                                         \
+};                                                            \
+template< class T, class A, typename Range >                  \
+struct resize_impl< THRUST_VECTOR<T,A> , Range >              \
+{                                                             \
+    static void resize( THRUST_VECTOR<T,A> &x ,               \
+                        const Range &y )                      \
+    {                                                         \
+        x.resize( thrust::distance(boost::begin(y),           \
+                                   boost::end(y)));           \
+    }                                                         \
+};                                                            \
 
 
-template< class T >
-struct is_resizeable< thrust::host_vector< T > >
-{
-    struct type : public boost::true_type { };
-    const static bool value = type::value;
-};
-
-template< class T >
-struct same_size_impl< thrust::host_vector< T > , thrust::host_vector< T > >
-{
-    static bool same_size( const thrust::host_vector< T > &x , const thrust::host_vector< T > &y )
-    {
-        return x.size() == y.size();
-    }
-};
-
-template< class T >
-struct resize_impl< thrust::host_vector< T > , thrust::host_vector< T > >
-{
-    static void resize( thrust::host_vector< T > &x , const thrust::host_vector< T > &y )
-    {
-        x.resize( y.size() );
-    }
-};
+#define ODEINT_THRUST_SAME_SIZE_IMPL( THRUST_VECTOR )            \
+template< class T , class A >                                    \
+struct same_size_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> > \
+{                                                                \
+    static bool same_size( const THRUST_VECTOR<T,A> &x ,         \
+                           const THRUST_VECTOR<T,A> &y )         \
+    {                                                            \
+        return x.size() == y.size();                             \
+    }                                                            \
+};                                                               \
+template< class T , class A, typename Range >                    \
+struct same_size_impl< THRUST_VECTOR<T,A> , Range >              \
+{                                                                \
+    static bool same_size( const THRUST_VECTOR<T,A> &x ,         \
+                           const Range &y )                      \
+    {                                                            \
+        return x.size() == thrust::distance(boost::begin(y),     \
+                                            boost::end(y));      \
+    }                                                            \
+};                                                               \
 
 
+#define ODEINT_THRUST_COPY_IMPL( THRUST_VECTOR )                        \
+template< class Container1 , class T , class A >                        \
+struct copy_impl< Container1 , THRUST_VECTOR<T,A> >                     \
+{                                                                       \
+    static void copy( const Container1 &from , THRUST_VECTOR<T,A> &to ) \
+    {                                                                   \
+        thrust::copy( boost::begin( from ) , boost::end( from ) ,       \
+                      boost::begin( to ) );                             \
+    }                                                                   \
+};                                                                      \
+                                                                        \
+template< class T , class A , class Container2 >                        \
+struct copy_impl< THRUST_VECTOR<T,A> , Container2 >                     \
+{                                                                       \
+    static void copy( const THRUST_VECTOR<T,A> &from , Container2 &to ) \
+    {                                                                   \
+        thrust::copy( boost::begin( from ) , boost::end( from ) ,       \
+                      boost::begin( to ) );                             \
+    }                                                                   \
+};                                                                      \
+                                                                        \
+template< class T , class A >                                           \
+struct copy_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> >             \
+{                                                                       \
+    static void copy( const THRUST_VECTOR<T,A> &from ,                  \
+                      THRUST_VECTOR<T,A> &to )                          \
+    {                                                                   \
+        thrust::copy( boost::begin( from ) , boost::end( from ) ,       \
+                      boost::begin( to ) );                             \
+    }                                                                   \
+};                                                                      \
 
-template< class Container1, class Value >
-struct copy_impl< Container1 , thrust::device_vector< Value > >
-{
-    static void copy( const Container1 &from , thrust::device_vector< Value > &to )
-    {
-        thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) );
-    }
-};
+// add support for the standard thrust containers
 
-template< class Value , class Container2 >
-struct copy_impl< thrust::device_vector< Value > , Container2 >
-{
-    static void copy( const thrust::device_vector< Value > &from , Container2 &to )
-    {
-        thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) );
-    }
-};
-
-template< class Value >
-struct copy_impl< thrust::device_vector< Value > , thrust::device_vector< Value > >
-{
-    static void copy( const thrust::device_vector< Value > &from , thrust::device_vector< Value > &to )
-    {
-        thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) );
-    }
-};
+ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::device_vector )
+ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::device_vector )
+ODEINT_THRUST_SAME_SIZE_IMPL( thrust::device_vector )
+ODEINT_THRUST_COPY_IMPL( thrust::device_vector )
 
 
+ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::host_vector )
+ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::host_vector )
+ODEINT_THRUST_SAME_SIZE_IMPL( thrust::host_vector )
+ODEINT_THRUST_COPY_IMPL( thrust::host_vector )
 
 
 } // odeint
 } // numeric
 } // boost
 
+// add support for thrust backend vectors, if available
+
+#include <thrust/version.h>
+
+#if THRUST_VERSION >= 100600
+
+#include <thrust/system/cpp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::cpp::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::cpp::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::cpp::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::cpp::vector )
+} } }
+
+#ifdef _OPENMP
+#include <thrust/system/omp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::omp::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::omp::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::omp::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::omp::vector )
+} } }
+#endif // _OPENMP
+
+#ifdef TBB_VERSION_MAJOR
+#include <thrust/system/tbb/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::tbb::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::tbb::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::tbb::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::tbb::vector )
+} } }
+#endif // TBB_VERSION_MAJOR
+
+#ifdef __CUDACC__
+#include <thrust/system/cuda/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::cuda::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::cuda::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::cuda::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::cuda::vector )
+} } }
+#endif // __CUDACC__
+
+#endif // THRUST_VERSION >= 100600
 
 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED

include/boost/numeric/odeint/integrate/integrate.hpp

@@ -48,6 +48,15 @@ integrate( System system , State &start_state , Time start_time , Time end_time
     return integrate_adaptive( stepper_type() , system , start_state , start_time , end_time , dt , observer );
 }
 
+template< class Value , class System , class State , class Time , class Observer >
+size_t 
+integrate( System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
+{
+    typedef controlled_runge_kutta< runge_kutta_dopri5< State , Value , State , Time > > stepper_type;
+    return integrate_adaptive( stepper_type() , system , start_state , start_time , end_time , dt , observer );
+}
+
+
 
 
 /*
@@ -59,6 +68,13 @@ size_t integrate( System system , State &start_state , Time start_time , Time en
     return integrate( system , start_state , start_time , end_time , dt , null_observer() );
 }
 
+template< class Value , class System , class State , class Time >
+size_t integrate( System system , State &start_state , Time start_time , Time end_time , Time dt )
+{
+    return integrate< Value >( system , start_state , start_time , end_time , dt , null_observer() );
+}
+
+
 
 /**
  * \fn integrate( System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
@@ -70,6 +86,9 @@ size_t integrate( System system , State &start_state , Time start_time , Time en
  * integration with step size control, thus dt changes during the integration.
  * This method uses standard error bounds of 1E-6.
  * After each step, the observer is called.
+ * 
+ * \attention A second version of this function template exists which explicitly
+ * expects the value type as template parameter, i.e. integrate< double >( sys , x , t0 , t1 , dt , obs );
  *
  * \param system The system function to solve, hence the r.h.s. of the 
  * ordinary differential equation.
@@ -92,6 +111,9 @@ size_t integrate( System system , State &start_state , Time start_time , Time en
  * integration with step size control, thus dt changes during the integration.
  * This method uses standard error bounds of 1E-6.
  * No observer is called.
+ * 
+ * \attention A second version of this function template exists which explicitly
+ * expects the value type as template parameter, i.e. integrate< double >( sys , x , t0 , t1 , dt );
  *
  * \param system The system function to solve, hence the r.h.s. of the 
  * ordinary differential equation.

include/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp

@@ -205,7 +205,7 @@ namespace odeint {
          * \param sys The system function (ODE) to solve.
          * \param s The initial state.
          */
-        adaptive_iterator_impl( stepper_type stepper , system_type sys , state_type &s )
+        adaptive_iterator_impl( stepper_type stepper , system_type sys , state_type& /* s */ )
             : base_type( stepper , sys ) { }
 
     protected:

include/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp

@@ -77,7 +77,7 @@ namespace odeint {
          * \param sys The system function (ODE) to solve.
          * \param s The initial state. const_step_iterator stores a reference of s and changes its value during the iteration.
          */
-        const_step_iterator_impl( stepper_type stepper , system_type sys , state_type &s )
+        const_step_iterator_impl( stepper_type stepper , system_type sys , state_type& /* s */ )
             : base_type( stepper , sys ) { }
 
     protected:

include/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp

@@ -48,7 +48,7 @@ Time integrate_n_steps(
         Observer observer , stepper_tag )
 {
     // ToDo: is there a better way to extract the final time?
-    Time t;
+    Time t = start_time; // Assignment is only here to avoid warnings.
     boost::for_each( make_n_step_time_range( stepper , system , start_state ,
                                              start_time , dt , num_of_steps ) ,
                      obs_caller_time< Observer , Time >( t , observer ) );
@@ -91,7 +91,7 @@ Time integrate_n_steps(
         Observer observer , dense_output_stepper_tag )
 {
     // ToDo: is there a better way to extract the final time?
-    Time t;
+    Time t = start_time;  // Assignment is only here to avoid warnings.
     boost::for_each( make_n_step_time_range( stepper , system , start_state ,
                                              start_time , dt , num_of_steps ) ,
                      obs_caller_time< Observer , Time >( t , observer ) );

include/boost/numeric/odeint/stepper/adams_bashforth.hpp

@@ -181,7 +181,8 @@ public :
         {
             if( i != 0 ) m_step_storage.rotate();
             sys( x , m_step_storage[0].m_v , t );
-            stepper.do_step( system , x , m_step_storage[0].m_v , t , dt );
+            stepper.do_step_dxdt_impl( system, x, m_step_storage[0].m_v, t,
+                                       dt );
             t += dt;
         }
         m_steps_initialized = steps;
@@ -222,7 +223,8 @@ private:
         {
             if( m_steps_initialized != 0 ) m_step_storage.rotate();
             sys( in , m_step_storage[0].m_v , t );
-            m_initializing_stepper.do_step( system , in , m_step_storage[0].m_v , t , out , dt );
+            m_initializing_stepper.do_step_dxdt_impl(
+                system, in, m_step_storage[0].m_v, t, out, dt );
             ++m_steps_initialized;
         }
         else

include/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp

@@ -49,7 +49,8 @@ class Deriv = State ,
 class Time = Value ,
 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
 class Operations = typename operations_dispatcher< State >::operations_type ,
-class Resizer = initially_resizer
+class Resizer = initially_resizer,
+class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer >
 >
 class adams_bashforth_moulton
 {
@@ -71,12 +72,13 @@ public :
     typedef Operations operations_type;
     typedef Resizer resizer_type;
     typedef stepper_tag stepper_category;
+    typedef InitializingStepper initializing_stepper_type;
 
     static const size_t steps = Steps;
 #ifndef DOXYGEN_SKIP
-    typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_bashforth_type;
+    typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type, initializing_stepper_type > adams_bashforth_type;
     typedef adams_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_moulton_type;
-    typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > stepper_type;
+    typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type , initializing_stepper_type> stepper_type;
 #endif //DOXYGEN_SKIP
     typedef unsigned short order_type;
     static const order_type order_value = steps;
@@ -158,7 +160,7 @@ private:
         {
             m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
             m_adams_bashforth.do_step( system , x , t , m_x.m_v , dt );
-            m_adams_moulton.do_step( system , x , m_x.m_v , t , x , dt , m_adams_bashforth.step_storage() );
+            m_adams_moulton.do_step( system , x , m_x.m_v , t+dt , x , dt , m_adams_bashforth.step_storage() );
         }
         else
         {

include/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp

@@ -151,6 +151,22 @@ public:
     }
 
 
+    /*
+     * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template < class System, class StateInOut, class DerivIn >
+    void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+                            time_type t, time_type dt )
+    {
+        this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
+    }
+
+
+
     /*
      * Version 3 : do_step( sys , in , t , out , dt )
      *
@@ -181,10 +197,21 @@ public:
     {
         this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
     }
-
-
-
-
+    
+    /*
+     * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template < class System, class StateIn, class DerivIn, class StateOut >
+    void do_step_dxdt_impl( System system, const StateIn &in,
+                            const DerivIn &dxdt, time_type t, StateOut &out,
+                            time_type dt )
+    {
+        this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
+    }
 
     /*
      * Version  5 :do_step( sys , x , t , dt , xerr )

include/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp

@@ -150,12 +150,28 @@ public:
     }
 
 
+    /*
+     * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template< class System , class StateInOut , class DerivInOut >
+    void do_step_dxdt_impl( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt )
+    {
+        m_first_call = true;
+        this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt );
+    }
+
     /*
      * version 3 : do_step( sys , in , t , out , dt )
      *
-     * this version does not solve the forwarding problem, boost.range can not be used
+     * this version does not solve the forwarding problem, boost.range can not
+     * be used.
      *
-     * the disable is needed to avoid ambiguous overloads if state_type = time_type
+     * the disable is needed to avoid ambiguous overloads if 
+     * state_type = time_type
      */
     template< class System , class StateIn , class StateOut >
     typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type
@@ -174,12 +190,14 @@ public:
      *
      * this version does not solve the forwarding problem, boost.range can not be used
      */
-    template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
-    void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t ,
-            StateOut &out , DerivOut &dxdt_out , time_type dt )
+    template< class System, class StateIn, class DerivIn, class StateOut,
+               class DerivOut >
+    void do_step( System system, const StateIn &in, const DerivIn &dxdt_in,
+                  time_type t, StateOut &out, DerivOut &dxdt_out, time_type dt )
     {
         m_first_call = true;
-        this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt );
+        this->stepper().do_step_impl( system, in, dxdt_in, t, out, dxdt_out,
+                                      dt );
     }
 
 

include/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp

@@ -138,6 +138,21 @@ public:
     }
 
 
+    /*
+     * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template < class System, class StateInOut, class DerivIn >
+    void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+                            time_type t, time_type dt )
+    {
+        this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
+    }
+
+
     /*
      * Version 3 : do_step( sys , in , t , out , dt )
      *
@@ -164,6 +179,22 @@ public:
         this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
     }
 
+
+    /*
+     * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version. Meant for internal use only.
+     */
+    template < class System, class StateIn, class DerivIn, class StateOut >
+    void do_step_dxdt_impl( System system, const StateIn &in,
+                            const DerivIn &dxdt, time_type t, StateOut &out,
+                            time_type dt )
+    {
+        this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
+    }
+
     template< class StateIn >
     void adjust_size( const StateIn &x )
     {

include/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp

@@ -182,7 +182,7 @@ private:
 
     // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
     template< class System , class StateIn , class StateOut >
-    void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , boost::mpl::true_ )
+    void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , boost::mpl::true_ )
     {
         typedef typename odeint::unwrap_reference< System >::type system_type;
         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;

include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp

@@ -194,7 +194,6 @@ public:
 
         static const value_type val1( 1.0 );
 
-        typename odeint::unwrap_reference< System >::type &sys = system;
         if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
         {
             reset(); // system resized -> reset
@@ -219,12 +218,12 @@ public:
             m_midpoint.set_steps( m_interval_sequence[k] );
             if( k == 0 )
             {
-                m_midpoint.do_step( sys , in , dxdt , t , out , dt );
+                m_midpoint.do_step( system , in , dxdt , t , out , dt );
                 /* the first step, nothing more to do */
             }
             else
             {
-                m_midpoint.do_step( sys , in , dxdt , t , m_table[k-1].m_v , dt );
+                m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt );
                 extrapolate( k , m_table , m_coeff , out );
                 // get error estimate
                 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,

include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp

@@ -155,8 +155,6 @@ public:
         
         static const value_type val1( 1.0 );
 
-        typename odeint::unwrap_reference< System >::type &sys = system;
-
         bool reject( true );
 
         time_vector h_opt( m_k_max+1 );
@@ -261,6 +259,7 @@ public:
         {
 
             //calculate dxdt for next step and dense output
+            typename odeint::unwrap_reference< System >::type &sys = system;
             sys( out , dxdt_new , t+dt );
 
             //prepare dense output

include/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp

@@ -117,7 +117,7 @@ struct adams_bashforth_call_algebra< 7 , Algebra , Operations >
     {
         //BOOST_ASSERT( false ); // not implemented
         typedef typename Coefficients::value_type value_type;
-        Algebra::for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v ,
+        algebra.for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v ,
                             typename Operations::template scale_sum8< value_type , Time , Time , Time , Time , Time , Time >(
                                 1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] ) );
     }
@@ -132,7 +132,7 @@ struct adams_bashforth_call_algebra< 8 , Algebra , Operations >
     {
         //BOOST_ASSERT( false ); // not implemented
         typedef typename Coefficients::value_type value_type;
-        Algebra::for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v ,
+        algebra.for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v ,
                              typename Operations::template scale_sum9< value_type , Time , Time , Time , Time , Time , Time , Time >(
                                  1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] , dt * coef[7] ) );
     }

include/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp

@@ -32,7 +32,7 @@ template< class Algebra , class Operations >
 struct adams_moulton_call_algebra< 1 , Algebra , Operations >
 {
     template< class StateIn , class StateOut , class DerivIn , class StepStorage , class Coefficients , class Time >
-    void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage &steps , const Coefficients &coef , Time dt ) const
+    void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage& /* steps */ , const Coefficients &coef , Time dt ) const
     {
         typedef typename Coefficients::value_type value_type;
         algebra.for_each3( out , in , dxdt , typename Operations::template scale_sum2< value_type , Time >( 1.0 , dt * coef[0] ) );

include/boost/numeric/odeint/stepper/euler.hpp

@@ -76,7 +76,7 @@ public :
     { }
 
     template< class System , class StateIn , class DerivIn , class StateOut >
-    void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
+    void do_step_impl( System /* system */ , const StateIn &in , const DerivIn &dxdt , time_type /* t */ , StateOut &out , time_type dt )
     {
         stepper_base_type::m_algebra.for_each3( out , in , dxdt ,
                 typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt ) );
@@ -84,7 +84,7 @@ public :
     }
 
     template< class StateOut , class StateIn1 , class StateIn2 >
-    void calc_state( StateOut &x , time_type t ,  const StateIn1 &old_state , time_type t_old , const StateIn2 &current_state , time_type t_new ) const
+    void calc_state( StateOut &x , time_type t ,  const StateIn1 &old_state , time_type t_old , const StateIn2 & /*current_state*/ , time_type /* t_new */ ) const
     {
         const time_type delta = t - t_old;
         stepper_base_type::m_algebra.for_each3( x , old_state , stepper_base_type::m_dxdt.m_v ,

include/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp

@@ -121,11 +121,20 @@ public :
 
         // dt * m_dxh = k4
         sys( m_x_tmp.m_v , m_dxh.m_v , t + dt );
+
         //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
         time_type dt6 = dt / static_cast< value_type >( 6 );
         time_type dt3 = dt / static_cast< value_type >( 3 );
         stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_dxt.m_v , m_dxm.m_v , m_dxh.m_v ,
-                typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) );
+                                             typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) );
+        
+        // x += dt/6 * m_dxdt + dt/3 * m_dxt )
+        // stepper_base_type::m_algebra.for_each4( out , in , dxdt , m_dxt.m_v , 
+        //                                         typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt6 , dt3 ) ); 
+        // // x += dt/3 * m_dxm + dt/6 * m_dxh )
+        // stepper_base_type::m_algebra.for_each4( out , out , m_dxm.m_v , m_dxh.m_v , 
+        //                                         typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt3 , dt6 ) ); 
+
     }
 
     template< class StateType >

include/boost/numeric/odeint/stepper/velocity_verlet.hpp

@@ -114,7 +114,7 @@ public:
 
         algebra_stepper_base_type::m_algebra.for_each4(
             qout , qin , pin , ain ,
-            typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one_half * dt , one * dt * dt ) );
+            typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) );
 
         typename odeint::unwrap_reference< System >::type & sys = system;
 

include/boost/numeric/odeint/util/bind.hpp

@@ -41,9 +41,12 @@ using namespace ::std::placeholders;
 
 #else
 
+// unnamed namespace to avoid multiple declarations (#138)
+namespace {
 using ::boost::bind;
 boost::arg<1> _1;
 boost::arg<2> _2;
+}
 // using ::boost::bind;
 // using ::_1;
 // using ::_2;

include/boost/numeric/odeint/util/resize.hpp

@@ -33,8 +33,8 @@
 namespace boost {
 namespace numeric {
 namespace odeint {
-    
-    
+
+
 template< class StateOut , class StateIn , class Enabler = void >
 struct resize_impl_sfinae
 {

include/boost/numeric/odeint/util/same_instance.hpp

@@ -25,7 +25,7 @@ namespace odeint {
 template< class T1 , class T2 , class Enabler=void >
 struct same_instance_impl
 { 
-    static bool same_instance( const T1 &x1 , const T2 &x2 )
+    static bool same_instance( const T1& /* x1 */ , const T2& /* x2 */ )
     {
         return false;
     }

include/boost/numeric/odeint/util/unwrap_reference.hpp

@@ -31,9 +31,9 @@
 namespace boost {
 
 #if BOOST_NUMERIC_ODEINT_CXX11
-template<typename T> class reference_wrapper;
+template<typename T> struct reference_wrapper;
 
-template<typename T> class unwrap_reference;
+template<typename T> struct unwrap_reference;
 #endif
 
 namespace numeric {
@@ -43,24 +43,21 @@ namespace odeint {
 #if BOOST_NUMERIC_ODEINT_CXX11
 
 template<typename T>
-class unwrap_reference
+struct unwrap_reference
 {
-public:
     typedef typename std::remove_reference<T>::type type;
 };
 
 template<typename T>
-class unwrap_reference< std::reference_wrapper<T> >
+struct unwrap_reference< std::reference_wrapper<T> >
 {
-public:
     typedef typename std::remove_reference<T>::type type;
 };
 
 template<typename T>
-class unwrap_reference< boost::reference_wrapper<T> >
+struct unwrap_reference< boost::reference_wrapper<T> >
 {
-public:
-        typedef typename boost::unwrap_reference<T>::type type;
+    typedef typename boost::unwrap_reference<T>::type type;
 };
 
 #else

performance/Jamfile.v2

@@ -11,6 +11,9 @@ project
     : requirements
       <define>BOOST_ALL_NO_LIB=1
       <include>../../../..
+      <cxxflags>-std=c++11
+      <toolset>gcc:<cxxflags>-ffast-math
+      <toolset>intel:<cxxflags>"-fast -inline-forceinline"
     : default-build release
     ;
     
@@ -24,49 +27,6 @@ lib libmkl_intel_thread : : <name>mkl_intel_thread ;
 lib libiomp5 : : <name>iomp5 ;
 lib libpthread : : <name>pthread ;
 
-exe odeint_rk4_lorenz_array
-    : odeint_rk4_lorenz_array.cpp
+exe odeint_rk4_array
+    : odeint_rk4_array.cpp
     ;
-    
-exe odeint_rk4_lorenz_range
-    : odeint_rk4_lorenz_range.cpp
-    ;
-
-# exe odeint_rk4_lorenz_fusion
-#     : odeint_rk4_lorenz_fusion.cpp
-#     ;
-
-exe generic_odeint_rk4_lorenz 
-    : generic_odeint_rk4_lorenz.cpp
-    : <toolset>intel:<cxxflags>-inline-forceinline
-    ;
-
-exe nr_rk4_lorenz
-    : nr_rk4_lorenz.cpp
-    ;
-
-exe rt_generic_rk4_lorenz
-    : rt_generic_rk4_lorenz.cpp
-    ;
-    
-exe gsl_rk4_lorenz
-    : gsl_rk4_lorenz.cpp libgslcblas libgsl
-    ;
-    
-exe odeint_rk4_phase_lattice
-    : odeint_rk4_phase_lattice.cpp
-    ;
-    
-exe odeint_rk4_phase_lattice_mkl
-    : odeint_rk4_phase_lattice_mkl.cpp libpthread libiomp5 libmkl_core libmkl_intel_thread libmkl
-    ;
-    
-exe nr_rk4_phase_lattice
-    : nr_rk4_phase_lattice.cpp
-    ;
-
-exe rt_generic_rk4_phase_lattice
-    : rt_generic_rk4_phase_lattice.cpp
-    ;
-
-

performance/Makefile

@@ -0,0 +1,43 @@
+# Copyright 2011-2014 Mario Mulansky
+# Copyright 2011-2014 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
+
+INCLUDES += -I../../include/ -I$(BOOST_ROOT)
+GCCFLAGS = -O3 -ffast-math -DNDEBUG
+# disabling -ffast-math might give slightly better performance
+ICCFLAGS = -Ofast -xHost -ip -inline-forceinline -DNDEBUG
+# Possible options: -fp-model source -no-fma
+GFORTFLAGS = -Ofast
+
+bin/gcc:
+	mkdir -p bin/gcc
+
+bin/intel:
+	mkdir -p bin/intel
+
+bin/gfort:
+	mkdir -p bin/gfort
+
+bin/gcc/odeint_rk4_array: odeint_rk4_array.cpp bin/gcc
+	g++ ${GCCFLAGS} ${INCLUDES} -o bin/gcc/odeint_rk4_array odeint_rk4_array.cpp
+
+bin/gcc/c_lorenz: c_lorenz.c bin/gcc
+	gcc -std=c99 -Ofast -mtune=corei7-avx c_lorenz.c -o bin/gcc/c_lorenz 
+
+bin/intel/odeint_rk4_array: odeint_rk4_array.cpp bin/intel
+	icpc ${ICCFLAGS} ${INCLUDES} -o bin/intel/odeint_rk4_array odeint_rk4_array.cpp
+
+bin/intel/c_lorenz: c_lorenz.c bin/intel
+	icc -std=c99 -Ofast -xHost -ansi-alias -o bin/intel/c_lorenz c_lorenz.c
+
+bin/gfort/fortran_lorenz: fortran_lorenz.f90 bin/gfort
+	gfortran ${GFORTFLAGS} fortran_lorenz.f90 -o bin/gfort/fortran_lorenz
+
+all: bin/gcc/odeint_rk4_array bin/intel/odeint_rk4_array bin/gcc/c_lorenz bin/intel/c_lorenz bin/gfort/fortran_lorenz

performance/SIMD/Makefile

@@ -0,0 +1,33 @@
+# Copyright 2014 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)
+
+# make sure BOOST_ROOT is pointing to your boost directory
+# otherwise, set it here:
+# BOOST_ROOT = /path/to/boost
+# you also need NT2s SIMD libary available set the include path here:
+# SIMD_INCLUDE = /path/to/simd/include
+
+INCLUDES = -I$(BOOST_ROOT) -I${SIMD_INCLUDE}
+
+# INTEL COMPILER
+# change this if you want to cross-compile
+ARCH = Host
+# ARCH = AVX
+# ARCH = SSE4.2
+
+CXX = icpc
+CC = icpc
+CXXFLAGS = -O3 -x${ARCH} -std=c++0x -fno-alias -inline-forceinline -DNDEBUG ${INCLUDES}
+# -ip
+
+# GCC COMPILER
+# change this if you want to cross-compile
+# ARCH = native
+# # ARCH = core-avx-i
+
+# CXX = g++
+# CC = g++
+# CXXFLAGS = -O3 -ffast-math -mtune=${ARCH} -march=${ARCH} -std=c++0x -DNDEBUG ${INCLUDES}

performance/SIMD/perf_roessler.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+echo "Running on ${HOSTNAME}"
+
+out_dir=perf_${HOSTNAME}
+mkdir -p ${out_dir}
+
+for N in 256 1024 4096 16384 65536 262144 1048576 4194304 16777216 67108864
+do
+    steps=`expr 4 \* 67108864 / ${N}`
+    for exe in "roessler" "roessler_simd"
+    do
+        rm -f ${out_dir}/${exe}_N${N}.times
+        for i in {0..4}
+        do
+            likwid-pin -cS0:0 ./${exe} ${N} ${steps} >> ${out_dir}/${exe}_N${N}.times
+        done
+        for perf_ctr in "FLOPS_DP" "FLOPS_AVX" "L2" "L3" "MEM"
+        do
+            likwid-perfctr -CS0:0 -g ${perf_ctr} ./${exe} ${N} ${steps} > ${out_dir}/${exe}_N${N}_${perf_ctr}.perf
+        done
+    done
+done

performance/SIMD/roessler.cpp

@@ -0,0 +1,125 @@
+/*
+ * Simulation of an ensemble of Roessler attractors
+ *
+ * Copyright 2014 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 <random>
+
+#include <boost/timer.hpp>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+namespace odeint = boost::numeric::odeint;
+
+typedef boost::timer timer_type;
+
+typedef double fp_type;
+//typedef float fp_type;
+
+typedef boost::array<fp_type, 3> state_type;
+typedef std::vector<state_type> state_vec;
+
+//---------------------------------------------------------------------------
+struct roessler_system {
+    const fp_type m_a, m_b, m_c;
+
+    roessler_system(const fp_type a, const fp_type b, const fp_type c)
+        : m_a(a), m_b(b), m_c(c)
+    {}
+
+    void operator()(const state_type &x, state_type &dxdt, const fp_type t) const
+    {
+        dxdt[0] = -x[1] - x[2];
+        dxdt[1] = x[0] + m_a * x[1];
+        dxdt[2] = m_b + x[2] * (x[0] - m_c);
+    }
+};
+
+//---------------------------------------------------------------------------
+int main(int argc, char *argv[]) {
+if(argc<3)
+{
+    std::cerr << "Expected size and steps as parameter" << std::endl;
+    exit(1);
+}
+const size_t n = atoi(argv[1]);
+const size_t steps = atoi(argv[2]);
+//const size_t steps = 50;
+
+const fp_type dt = 0.01;
+
+const fp_type a = 0.2;
+const fp_type b = 1.0;
+const fp_type c = 9.0;
+
+// random initial conditions on the device
+std::vector<fp_type> x(n), y(n), z(n);
+std::default_random_engine generator;
+std::uniform_real_distribution<fp_type> distribution_xy(-8.0, 8.0);
+std::uniform_real_distribution<fp_type> distribution_z(0.0, 20.0);
+auto rand_xy = std::bind(distribution_xy, std::ref(generator));
+auto rand_z = std::bind(distribution_z, std::ref(generator));
+std::generate(x.begin(), x.end(), rand_xy);
+std::generate(y.begin(), y.end(), rand_xy);
+std::generate(z.begin(), z.end(), rand_z);
+
+state_vec state(n);
+for(size_t i=0; i<n; ++i)
+{
+    state[i][0] = x[i];
+    state[i][1] = y[i];
+    state[i][2] = z[i];
+}
+
+std::cout.precision(16);
+
+std::cout << "# n: " << n << std::endl;
+
+std::cout << x[0] << std::endl;
+
+
+// Stepper type - use never_resizer for slight performance improvement
+odeint::runge_kutta4_classic<state_type, fp_type, state_type, fp_type,
+                             odeint::array_algebra,
+                             odeint::default_operations,
+                             odeint::never_resizer> stepper;
+
+roessler_system sys(a, b, c);
+
+timer_type timer;
+
+fp_type t = 0.0;
+
+for (int step = 0; step < steps; step++)
+{
+    for(size_t i=0; i<n; ++i)
+    {
+        stepper.do_step(sys, state[i], t, dt);
+    }
+    t += dt;
+}
+
+std::cout << "Integration finished, runtime for " << steps << " steps: ";
+std::cout << timer.elapsed() << " s" << std::endl;
+
+// compute some accumulation to make sure all results have been computed
+fp_type s = 0.0;
+for(size_t i = 0; i < n; ++i)
+{
+    s += state[i][0];
+}
+
+std::cout << state[0][0] << std::endl;
+std::cout << s/n << std::endl;
+
+}

performance/SIMD/roessler_simd.cpp

@@ -0,0 +1,149 @@
+/*
+ * Simulation of an ensemble of Roessler attractors using NT2 SIMD library
+ * This requires the SIMD library headers.
+ *
+ * Copyright 2014 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 <random>
+
+#include <boost/timer.hpp>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+#include <boost/simd/sdk/simd/pack.hpp>
+#include <boost/simd/sdk/simd/io.hpp>
+#include <boost/simd/memory/allocator.hpp>
+#include <boost/simd/include/functions/splat.hpp>
+#include <boost/simd/include/functions/plus.hpp>
+#include <boost/simd/include/functions/multiplies.hpp>
+
+
+namespace odeint = boost::numeric::odeint;
+namespace simd = boost::simd;
+
+typedef boost::timer timer_type;
+
+static const size_t dim = 3;  // roessler is 3D
+
+typedef double fp_type;
+//typedef float fp_type;
+ 
+typedef simd::pack<fp_type> simd_pack;
+typedef boost::array<simd_pack, dim> state_type;
+// use the simd allocator to get properly aligned memory
+typedef std::vector< state_type, simd::allocator< state_type > > state_vec;
+
+static const size_t pack_size = simd_pack::static_size;
+
+//---------------------------------------------------------------------------
+struct roessler_system {
+    const fp_type m_a, m_b, m_c;
+
+    roessler_system(const fp_type a, const fp_type b, const fp_type c)
+        : m_a(a), m_b(b), m_c(c)
+    {}
+
+    void operator()(const state_type &x, state_type &dxdt, const fp_type t) const
+    {
+        dxdt[0] = -1.0*x[1] - x[2];
+        dxdt[1] = x[0] + m_a * x[1];
+        dxdt[2] = m_b + x[2] * (x[0] - m_c);
+    }
+};
+
+//---------------------------------------------------------------------------
+int main(int argc, char *argv[]) {
+if(argc<3)
+{
+    std::cerr << "Expected size and steps as parameter" << std::endl;
+    exit(1);
+}
+const size_t n = atoi(argv[1]);
+const size_t steps = atoi(argv[2]);
+
+const fp_type dt = 0.01;
+
+const fp_type a = 0.2;
+const fp_type b = 1.0;
+const fp_type c = 9.0;
+
+// random initial conditions on the device
+std::vector<fp_type> x(n), y(n), z(n);
+std::default_random_engine generator;
+std::uniform_real_distribution<fp_type> distribution_xy(-8.0, 8.0);
+std::uniform_real_distribution<fp_type> distribution_z(0.0, 20.0);
+auto rand_xy = std::bind(distribution_xy, std::ref(generator));
+auto rand_z = std::bind(distribution_z, std::ref(generator));
+std::generate(x.begin(), x.end(), rand_xy);
+std::generate(y.begin(), y.end(), rand_xy);
+std::generate(z.begin(), z.end(), rand_z);
+
+state_vec state(n/pack_size);
+for(size_t i=0; i<n/pack_size; ++i)
+{
+    for(size_t p=0; p<pack_size; ++p)
+    {
+        state[i][0][p] = x[i*pack_size+p];
+        state[i][1][p] = y[i*pack_size+p];
+        state[i][2][p] = z[i*pack_size+p];
+    }
+}
+
+std::cout << "Systems: " << n << std::endl;
+std::cout << "Steps: " << steps << std::endl;
+std::cout << "SIMD pack size: " << pack_size << std::endl;
+
+std::cout << state[0][0] << std::endl;
+
+// Stepper type
+odeint::runge_kutta4_classic<state_type, fp_type, state_type, fp_type,
+                             odeint::array_algebra, odeint::default_operations,
+                             odeint::never_resizer> stepper;
+
+roessler_system sys(a, b, c);
+
+timer_type timer;
+
+fp_type t = 0.0;
+
+for(int step = 0; step < steps; step++)
+{
+    for(size_t i = 0; i < n/pack_size; ++i)
+    {
+        stepper.do_step(sys, state[i], t, dt);
+    }
+    t += dt;
+}
+
+std::cout.precision(16);
+
+std::cout << "Integration finished, runtime for " << steps << " steps: ";
+std::cout << timer.elapsed() << " s" << std::endl;
+
+// compute some accumulation to make sure all results have been computed
+simd_pack s_pack = 0.0;
+for(size_t i = 0; i < n/pack_size; ++i)
+{
+    s_pack += state[i][0];
+}
+
+fp_type s = 0.0;
+for(size_t p=0; p<pack_size; ++p)
+{
+    s += s_pack[p];
+}
+
+
+std::cout << state[0][0] << std::endl;
+std::cout << s/n << std::endl;
+
+}

performance/c_lorenz.c

@@ -0,0 +1,57 @@
+#include <stdio.h>
+#include <time.h>
+#include <math.h>
+
+void lorenz(const double *x, double *restrict y) {
+    y[0] = 10.0 * (x[1] - x[0]);
+    y[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
+    y[2] = x[0] * x[1] - (8.0 / 3.0) * x[2];
+}
+
+int main(int argc, const char *argv[])
+{
+    const int nb_steps = 20000000;
+    const double h = 1.0e-10;
+    const double h2 = 0.5 * h;
+    const double nb_loops = 21;
+    double x[3];
+    double y[3];
+    double f1[3];
+    double f2[3];
+    double f3[3];
+    double f4[3];
+    double min_time = 1E6;
+    clock_t begin, end;
+    double time_spent;
+    
+    for (int j = 0; j < nb_loops; j++) {
+        x[0] = 8.5;
+        x[1] = 3.1;
+        x[2] = 1.2;
+        begin = clock();
+        for (int k = 0; k < nb_steps; k++) {
+            lorenz(x, f1);
+            for (int i = 0; i < 3; i++) {
+                y[i] = x[i] + h2 * f1[i];
+            }
+            lorenz(y, f2);
+            for (int i = 0; i < 3; i++) {
+                y[i] = x[i] + h2 * f2[i];
+            }
+            lorenz(y, f3);
+            for (int i = 0; i < 3; i++) {
+                y[i] = x[i] + h * f3[i];
+            }
+            lorenz(y, f4);
+            for (int i = 0; i < 3; i++) {
+                x[i] = x[i] + h * (f1[i] + 2 * (f2[i] + f3[i]) + f4[i]) / 6.0;
+            }
+        }
+        end = clock();
+        min_time = fmin(min_time, (double)(end-begin)/CLOCKS_PER_SEC);
+        printf("Result: %f\t runtime: %f\n", x[0], (double)(end-begin)/CLOCKS_PER_SEC);
+    }
+    printf("Minimal Runtime: %f\n", min_time);
+    
+    return 0;
+}

performance/fortran_lorenz.f90

@@ -0,0 +1,60 @@
+program main
+    implicit none
+
+    integer, parameter :: dp = 8
+    real(dp), dimension(1:3) :: x
+    integer, parameter :: nstep = 20000000
+    real(dp) :: t = 0.0_dp
+    real(dp) :: h = 1.0e-10_dp
+    integer, parameter :: nb_loops = 21
+    integer, parameter :: n = 3
+    integer :: k
+    integer :: time_begin
+    integer :: time_end
+    integer :: count_rate
+    real(dp) :: time
+    real(dp) :: min_time = 100.0
+
+    do k = 1, nb_loops
+        x = [ 8.5_dp, 3.1_dp, 1.2_dp ]
+        call system_clock(time_begin, count_rate)
+        call rk4sys(n, t, x, h, nstep)
+        call system_clock(time_end, count_rate)
+        time = real(time_end - time_begin, dp) / real(count_rate, dp)
+        min_time = min(time, min_time)
+        write (*,*) time, x(1)
+    end do
+    write (*,*) "Minimal Runtime:", min_time
+contains
+    subroutine xpsys(x,f)
+        real(dp), dimension(1:3), intent(in) :: x
+        real(dp), dimension(1:3), intent(out) :: f
+        f(1) = 10.0_dp * ( x(2) - x(1) )
+        f(2) = 28.0_dp * x(1) - x(2) - x(1) * x(3)
+        f(3) = x(1) * x(2) - (8.0_dp / 3.0_dp) * x(3)
+    end subroutine xpsys
+
+    subroutine rk4sys(n, t, x, h, nstep)
+        integer, intent(in) :: n
+        real(dp), intent(in) :: t
+        real(dp), dimension(1:n), intent(inout) :: x
+        real(dp), intent(in) :: h
+        integer, intent(in) :: nstep
+        ! Local variables
+        real(dp) :: h2
+        real(dp), dimension(1:n) :: y, f1, f2, f3, f4
+        integer :: i, k
+
+        h2 = 0.5_dp * h
+        do k = 1, nstep
+            call xpsys(x, f1)
+            y = x + h2 * f1
+            call xpsys(y, f2)
+            y = x + h2 * f2
+            call xpsys(y, f3)
+            y = x + h * f3
+            call xpsys(y, f4)
+            x = x + h * (f1 + 2.0_dp * (f2 + f3) + f4) / 6.0_dp
+        end do
+    end subroutine rk4sys
+end program main

performance/fusion_algebra.hpp

@@ -1,202 +0,0 @@
-/*
- * fusion_algebra.hpp
- *
- * 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)
- */
-
-
-#ifndef FUSION_ALGEBRA_HPP_
-#define FUSION_ALGEBRA_HPP_
-
-#include <boost/array.hpp>
-
-#include <iostream>
-
-
-template< size_t n >
-struct fusion_algebra
-{
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
-            const boost::array< double , n > &a ,
-            const boost::array< T , dim > k_vector[n] , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i];// + a[0]*dt*k_vector[0][i];
-            for( size_t j = 0 ; j<n ; ++j )
-                x_tmp[i] += a[j]*dt*k_vector[j][i];
-        }
-    }
-
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp ,
-                const boost::array< double , n > &a ,
-                const boost::array< T , dim > k_vector[n] , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = a[0]*dt*k_vector[0][i];
-            for( size_t j = 1 ; j<n ; ++j )
-                x_tmp[i] += a[j]*dt*k_vector[j][i];
-         }
-    }
-
-};
-
-
-
-
-/** hand-wise implementation for performance improvement for n = 1..4 **/
-
-/* !!!!!!!   Actually, this is factor 3 slower with intel compiler, so we don'y use it !!!!!
- * Update: It increases performance on msvc 9.0 by about 30%, so it is activated for MSVC
- */
-
-//#ifdef BOOST_MSVC
-
-template<>
-struct fusion_algebra< 1 >
-{
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
-            const boost::array< double , 1 > &a ,
-            const boost::array< T , dim > *k_vector , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i]
-                + a[0]*dt*k_vector[0][i];
-        }
-    }
-
-};
-
-
-template<>
-struct fusion_algebra< 2 >
-{
-
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
-            const boost::array< double , 2 > &a ,
-            const boost::array< T , dim > *k_vector , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i]
-             + a[0]*dt*k_vector[0][i]
-             + a[1]*dt*k_vector[1][i];
-        }
-    }
-
-};
-
-
-template<>
-struct fusion_algebra< 3 >
-{
-
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
-            const boost::array< double , 3 > &a ,
-            const boost::array< T , dim > *k_vector , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i]
-             + a[0]*dt*k_vector[0][i]
-             + a[1]*dt*k_vector[1][i]
-             + a[2]*dt*k_vector[2][i];
-        }
-    }
-
-};
-
-template<>
-struct fusion_algebra< 4 >
-{
-
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
-            const boost::array< double , 4 > &a ,
-            const boost::array< T , dim > *k_vector , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i]
-             + a[0]*dt*k_vector[0][i]
-             + a[1]*dt*k_vector[1][i]
-             + a[2]*dt*k_vector[2][i]
-             + a[3]*dt*k_vector[3][i];
-        }
-    }
-
-};
-
-template<>
-struct fusion_algebra< 5 >
-{
-
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
-            const boost::array< double , 5 > &a ,
-            const boost::array< T , dim > *k_vector , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i]
-             + a[0]*dt*k_vector[0][i]
-             + a[1]*dt*k_vector[1][i]
-             + a[2]*dt*k_vector[2][i]
-             + a[3]*dt*k_vector[3][i]
-             + a[4]*dt*k_vector[4][i];
-        }
-    }
-
-};
-
-template<>
-struct fusion_algebra< 6 >
-{
-
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > &x_tmp , const boost::array< T , dim > &x ,
-            const boost::array< double , 6 > &a ,
-            const boost::array< T , dim > *k_vector , const double dt )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i]
-             + a[0]*dt*k_vector[0][i]
-             + a[1]*dt*k_vector[1][i]
-             + a[2]*dt*k_vector[2][i]
-             + a[3]*dt*k_vector[3][i]
-             + a[4]*dt*k_vector[4][i]
-             + a[5]*dt*k_vector[5][i];
-        }
-    }
-
-    template< typename T , size_t dim >
-    inline static void foreach(boost::array<T , dim> &x_tmp ,
-            const boost::array<double , 6> &a ,
-            const boost::array<T , dim> *k_vector , const double dt)
-    {
-        for (size_t i = 0 ; i < dim ; ++i)
-        {
-            x_tmp[i] = a[0] * dt * k_vector[0][i] + a[1] * dt * k_vector[1][i]
-                    + a[2] * dt * k_vector[2][i] + a[3] * dt * k_vector[3][i]
-                    + a[4] * dt * k_vector[4][i] + a[5] * dt * k_vector[5][i];
-        }
-    }
-
-};
-
-//#endif /* BOOST_MSVC */
-
-#endif /* FUSION_ALGEBRA_HPP_ */

performance/fusion_explicit_error_rk.hpp

@@ -1,63 +0,0 @@
-/*
- * fusion_explicit_error_rk.hpp
- *
- * 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)
- */
-
-
-#ifndef FUSION_EXPLICIT_ERROR_RK_HPP_
-#define FUSION_EXPLICIT_ERROR_RK_HPP_
-
-#include "fusion_explicit_rk_new.hpp"
-#include "fusion_algebra.hpp"
-
-namespace mpl = boost::mpl;
-namespace fusion = boost::fusion;
-
-using namespace std;
-
-template< class StateType , size_t stage_count >
-class explicit_error_rk : public explicit_rk< StateType , stage_count >
-{
-
-public:
-
-    typedef explicit_rk< StateType , stage_count > base;
-
-    typedef StateType state_type;
-
-    typedef typename base::stage_indices stage_indices;
-
-    typedef typename base::coef_a_type coef_a_type;
-
-    typedef typename base::coef_b_type coef_b_type;
-    typedef typename base::coef_c_type coef_c_type;
-
- public:
-
-    explicit_error_rk( const coef_a_type &a ,
-                          const coef_b_type &b ,
-                          const coef_b_type &b2 ,
-                          const coef_c_type &c )
-    : base( a , b , c ) , m_b2( b2 )
-    { }
-
-    template< class System >
-    void inline do_step( System system , state_type &x , const double t , const double dt , state_type &x_err )
-    {
-        base::do_step( system , x , t , dt );
-        // compute error estimate
-        fusion_algebra< stage_count >::foreach( x_err , m_b2 , base::m_F , dt );
-    }
-
-private:
-
-    const coef_b_type m_b2;
-};
-
-#endif /* FUSION_EXPLICIT_ERROR_RK_HPP_ */

performance/fusion_explicit_rk_new.hpp

@@ -1,217 +0,0 @@
-/*
- * fusion_runge_kutta.hpp
- *
- * Copyright 2010-2011 Mario Mulansky
- * Copyright 2010-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)
- */
-
-
-#ifndef FUSION_EXPLICIT_RK_HPP_
-#define FUSION_EXPLICIT_RK_HPP_
-
-#include <boost/mpl/vector.hpp>
-#include <boost/mpl/push_back.hpp>
-#include <boost/mpl/for_each.hpp>
-#include <boost/mpl/range_c.hpp>
-#include <boost/mpl/copy.hpp>
-#include <boost/mpl/size_t.hpp>
-
-#include <boost/fusion/container.hpp>
-#include <boost/fusion/algorithm/iteration.hpp>
-
-#include <boost/array.hpp>
-
-#include "fusion_algebra.hpp"
-//#include "fusion_foreach_performance.hpp"
-
-namespace mpl = boost::mpl;
-namespace fusion = boost::fusion;
-
-using namespace std;
-
-struct intermediate_stage {};
-struct last_stage {};
-
-
-
-template< class T , class Constant >
-struct array_wrapper
-{
-    typedef const typename boost::array< T , Constant::value > type;
-};
-
-template< class T , size_t i , class StageCategory >
-struct stage
-{
-    T c;
-    boost::array< T , i > a;
-    typedef StageCategory category;
-};
-
-template< class T , size_t i>
-struct stage< T , i , last_stage >
-{
-    T c;
-    boost::array< T , i > b;
-    typedef last_stage category;
-};
-
-
-
-template< class T , class Constant , class StageCategory >
-struct stage_wrapper
-{
-    typedef stage< T , Constant::value , StageCategory > type;
-};
-
-
-template< class StateType , size_t stage_count >
-class explicit_rk
-{
-
-public:
-
-    typedef StateType state_type;
-
-    typedef mpl::range_c< size_t , 1 , stage_count > stage_indices;
-
-    typedef typename fusion::result_of::as_vector
-    <
-        typename mpl::copy
-        <
-            stage_indices ,
-            mpl::inserter
-            <
-                mpl::vector0< > ,
-                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
-            >
-        >::type
-    >::type coef_a_type;
-
-    typedef boost::array< double , stage_count > coef_b_type;
-    typedef boost::array< double , stage_count > coef_c_type;
-
-    typedef typename fusion::result_of::as_vector
-    <
-        typename mpl::push_back
-        <
-            typename mpl::copy
-            <
-                stage_indices,
-                mpl::inserter
-                <
-                    mpl::vector0<> ,
-                    mpl::push_back< mpl::_1 , stage_wrapper< double , mpl::_2 , intermediate_stage > >
-                >
-            >::type ,
-            stage< double , stage_count , last_stage >
-        >::type
-    >::type stage_vector_base;
-
-
-    struct stage_vector : public stage_vector_base
-    {
-        struct do_insertion
-        {
-            stage_vector_base &m_base;
-            const coef_a_type &m_a;
-            const coef_c_type &m_c;
-
-            do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
-            : m_base( base ) , m_a( a ) , m_c( c ) { }
-
-            template< class Index >
-            void operator()( Index ) const
-            {
-                //fusion::at< Index >( m_base ) = stage< double , Index::value+1 , intermediate_stage >( m_c[ Index::value ] , fusion::at< Index >( m_a ) );
-                fusion::at< Index >( m_base ).c  = m_c[ Index::value ];
-                fusion::at< Index >( m_base ).a = fusion::at< Index >( m_a );
-            }
-        };
-
-        stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
-        {
-            typedef mpl::range_c< size_t , 0 , stage_count - 1 > indices;
-            mpl::for_each< indices >( do_insertion( *this , a , c ) );
-            //fusion::at_c< 0 >( fusion::at_c< stage_count - 1 >( *this ) ) = stage_count - 1 ;
-            fusion::at_c< stage_count - 1 >( *this ).c = c[ stage_count - 1 ];
-            fusion::at_c< stage_count - 1 >( *this ).b = b;
-        }
-    };
-
-
-
-    template< class System >
-    struct calculate_stage
-    {
-        System &system;
-        state_type &x , &x_tmp;
-        state_type *F;
-        const double t;
-        const double dt;
-
-        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_F ,
-                            const double _t , const double _dt )
-        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , F( _F ) , t( _t ) , dt( _dt )
-        {}
-
-
-        template< typename T , size_t stage_number >
-        void inline operator()( stage< T , stage_number , intermediate_stage > const &stage ) const
-        //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
-        {
-            if( stage_number == 1 )
-                system( x , F[stage_number-1] , t + stage.c * dt );
-            else
-                system( x_tmp , F[stage_number-1] , t + stage.c * dt );
-
-            fusion_algebra<stage_number>::foreach( x_tmp , x , stage.a , F , dt);
-        }
-
-
-        template< typename T , size_t stage_number >
-        void inline operator()( stage< T , stage_number , last_stage > const &stage ) const
-        //void operator()( typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , last_stage >::type const &stage ) const
-        {
-            if( stage_number == 1 )
-                system( x , F[stage_number-1] , t + stage.c * dt );
-            else
-                system( x_tmp , F[stage_number-1] , t + stage.c * dt );
-
-            fusion_algebra<stage_number>::foreach( x , x , stage.b , F , dt);
-        }
-
-
-    };
-
-public:
-
-    explicit_rk( const coef_a_type &a ,
-                  const coef_b_type &b ,
-                  const coef_c_type &c )
-    : m_stages( a , b , c )
-
-    { }
-
-
-    template< class System >
-    void inline do_step( System system , state_type &x , const double t , const double dt )
-    {
-        fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_F , t , dt ) );
-    }
-
-private:
-
-    stage_vector m_stages;
-    state_type m_x_tmp;
-
-protected:
-    state_type m_F[stage_count];
-
-};
-
-#endif /* FUSION_EXPLICIT_RK_HPP_ */

performance/generic_odeint_rk4_lorenz.cpp

@@ -1,84 +0,0 @@
-/*
- * 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)
- */
-
-
-#include <boost/array.hpp>
-
-//#include <boost/numeric/odeint/stepper/explicit_generic_rk.hpp>
-#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
-#include <boost/numeric/odeint/algebra/array_algebra.hpp>
-
-#include "rk_performance_test_case.hpp"
-
-#include "lorenz.hpp"
-
-using namespace boost::numeric::odeint;
-
-typedef boost::array< double , 3 > state_type;
-
-/*
-typedef explicit_generic_rk< 4 , 4 , state_type , double , state_type , double , array_algebra > rk4_type;
-
-typedef rk4_type::coef_a_type coef_a_type;
-typedef rk4_type::coef_b_type coef_b_type;
-typedef rk4_type::coef_c_type coef_c_type;
-
-const boost::array< double , 1 > a1 = {{ 0.5 }};
-const boost::array< double , 2 > a2 = {{ 0.0 , 0.5 }};
-const boost::array< double , 3 > a3 = {{ 0.0 , 0.0 , 1.0 }};
-
-const coef_a_type a = fusion::make_vector( a1 , a2 , a3 );
-const coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
-const coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
-*/
-
-typedef runge_kutta4< state_type , double , state_type , double , array_algebra > rk4_type;
-
-
-class rk4_wrapper
-{
-
-public:
-
-    rk4_wrapper()
-    //        : m_stepper( a , b , c ) 
-    {}
-
-    void reset_init_cond()
-    {
-        m_x[0] = 10.0 * rand() / RAND_MAX;
-        m_x[1] = 10.0 * rand() / RAND_MAX;
-        m_x[2] = 10.0 * rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        m_stepper.do_step( lorenz(), m_x , m_t , dt );
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-    rk4_type m_stepper;
-};
-
-
-
-int main()
-{
-    srand( 12312354 );
-
-    rk4_wrapper stepper;
-
-    run( stepper );
-}

performance/gsl_rk4_lorenz.cpp

@@ -1,71 +0,0 @@
-/*
- * gsl_rk4_lorenz.cpp
- *
- * 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)
- */
-
-
-#include <gsl/gsl_odeiv.h>
-
-#include "rk_performance_test_case.hpp"
-
-#include "lorenz_gsl.hpp"
-
-const size_t dim = 3;
-
-class gsl_wrapper
-{
-public:
-
-    gsl_wrapper()
-    {
-        m_s = gsl_odeiv_step_alloc( gsl_odeiv_step_rk4 , dim);
-        m_sys.function = lorenz_gsl;
-        m_sys.jacobian = 0;
-        m_sys.dimension = dim;
-        m_sys.params = 0;
-    }
-
-    void reset_init_cond()
-    {
-        m_x[0] = 10.0 * rand() / RAND_MAX;
-        m_x[1] = 10.0 * rand() / RAND_MAX;
-        m_x[2] = 10.0 * rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        gsl_odeiv_step_apply ( m_s , m_t , dt , m_x , m_x_err , 0 , 0 , &m_sys );
-        //m_t += dt;
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-    ~gsl_wrapper()
-    {
-        gsl_odeiv_step_free( m_s );
-    }
-
-private:
-    double m_x[dim];
-    double m_x_err[dim];
-    double m_t;
-    gsl_odeiv_step *m_s;
-    gsl_odeiv_system m_sys;
-};
-
-
-
-int main()
-{
-    gsl_wrapper stepper;
-
-    run( stepper , 20000000 / 3 , 1E-10 * 3);
-}

performance/lorenz.hpp

@@ -30,17 +30,4 @@ struct lorenz
 };
 
 
-typedef boost::array< double , 3 > state_type;
-
-
-inline void lorenz_func( 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] = x[0]*x[1] - b * x[2];
-}
-
 #endif /* LORENZ_HPP_ */

performance/lorenz_gsl.hpp

@@ -1,30 +0,0 @@
-/*
- * lorenz_gsl.hpp
- *
- * 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)
- */
-
-
-#ifndef LORENZ_GSL_HPP_
-#define LORENZ_GSL_HPP_
-
-#include <gsl/gsl_matrix.h>
-
-int lorenz_gsl( const double t , const double y[] , double f[] , void *params)
-{
-    const double sigma = 10.0;
-    const double R = 28.0;
-    const double b = 8.0 / 3.0;
-
-    f[0] = sigma * ( y[1] - y[0] );
-    f[1] = R * y[0] - y[1] - y[0] * y[2];
-    f[2] = y[0]*y[1] - b * y[2];
-    return GSL_SUCCESS;
-}
-
-#endif

performance/mpi/Jamfile.v2

@@ -1,20 +0,0 @@
-# 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)
-
-import mpi ;
-
-project
-    : requirements
-      <include>../../../../..
-      <include>..
-      <define>BOOST_ALL_NO_LIB=1
-      <library>/boost//timer
-      <library>/boost//mpi
-      <library>/boost//program_options
-    ;
-
-exe osc_chain_1d : osc_chain_1d.cpp ;

performance/mpi/osc_chain_1d.cpp

@@ -1,112 +0,0 @@
-/* Boost libs/numeric/odeint/performance/openmp/osc_chain_1d.cpp
-
- Copyright 2013 Karsten Ahnert
- Copyright 2013 Mario Mulansky
- Copyright 2013 Pascal Germroth
-
- stronlgy nonlinear hamiltonian lattice in 2d
-
- 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>
-#include <boost/numeric/odeint/external/mpi/mpi.hpp>
-
-#include <boost/program_options.hpp>
-#include <boost/random.hpp>
-#include <boost/timer/timer.hpp>
-#include <boost/foreach.hpp>
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics/stats.hpp>
-#include <boost/accumulators/statistics/mean.hpp>
-#include <boost/accumulators/statistics/median.hpp>
-#include "osc_chain_1d_system.hpp"
-
-using namespace std;
-using namespace boost::numeric::odeint;
-using namespace boost::accumulators;
-using namespace boost::program_options;
-
-using boost::timer::cpu_timer;
-
-const double p_kappa = 3.3;
-const double p_lambda = 4.7;
-
-int main( int argc , char* argv[] )
-{
-    boost::mpi::environment env(argc, argv);
-    boost::mpi::communicator world;
-
-    size_t N, steps, repeat;
-    bool dump;
-    options_description desc("Options");
-    desc.add_options()
-        ("help,h", "show this help")
-        ("length", value(&N)->default_value(1024), "length of chain")
-        ("steps", value(&steps)->default_value(100), "simulation steps")
-        ("repeat", value(&repeat)->default_value(25), "repeat runs")
-        ("dump", bool_switch(&dump), "dump final state to stderr (on node 0)")
-        ;
-    variables_map vm;
-    store(command_line_parser(argc, argv).options(desc).run(), vm);
-    notify(vm);
-    if(vm.count("help"))
-    {
-        if(world.rank() == 0)
-            cerr << desc << endl;
-        return EXIT_FAILURE;
-    }
-    cout << "length\tsteps\tthreads\ttime" << endl;
-
-    accumulator_set< double, stats<tag::mean, tag::median> > acc_time;
-
-    vector<double> p( N ), q( N, 0 );
-    if(world.rank() == 0) {
-        boost::random::uniform_real_distribution<double> distribution;
-        boost::random::mt19937 engine( 0 );
-        generate( p.begin() , p.end() , boost::bind( distribution , engine ) );
-    }
-
-    typedef vector<double> inner_state_type;
-    typedef mpi_state< inner_state_type > state_type;
-    typedef symplectic_rkn_sb3a_mclachlan<
-              state_type , state_type , double
-            > stepper_type;
-    state_type p_split( world ), q_split( world );
-    split(p, p_split);
-    split(q, q_split);
-
-    for(size_t n_run = 0 ; n_run != repeat ; n_run++) {
-        cpu_timer timer;
-        world.barrier();
-        integrate_n_steps( stepper_type() , osc_chain( p_kappa , p_lambda ) ,
-                           make_pair( boost::ref(q_split) , boost::ref(p_split) ) ,
-                           0.0 , 0.01 , steps );
-        world.barrier();
-        if(world.rank() == 0) {
-            double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
-            acc_time(run_time);
-            cout << N << '\t' << steps << '\t' << world.size() << '\t' << run_time << endl;
-        }
-    }
-
-    if(dump) {
-        unsplit(p_split, p);
-        if(world.rank() == 0) {
-            copy(p.begin(), p.end(), ostream_iterator<double>(cerr, "\t"));
-            cerr << endl;
-        }
-    }
-
-    if(world.rank() == 0)
-        cout << "# mean=" << mean(acc_time)
-             << " median=" << median(acc_time) << endl;
-
-    return 0;
-}
-

performance/mpi/osc_chain_1d_system.hpp

@@ -1,88 +0,0 @@
-/* Boost libs/numeric/odeint/performance/openmp/osc_chain_1d_system.hpp
-
- Copyright 2013 Karsten Ahnert
- Copyright 2013 Mario Mulansky
- Copyright 2013 Pascal Germroth
-
- stronlgy nonlinear hamiltonian lattice
-
- 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 SYSTEM_HPP
-#define SYSTEM_HPP
-
-#include <vector>
-#include <cmath>
-#include <iostream>
-
-#include <boost/math/special_functions/sign.hpp>
-#include <boost/numeric/odeint/external/mpi/mpi.hpp>
-
-namespace checked_math {
-    inline double pow( double x , double y )
-    {
-        if( x==0.0 )
-            // 0**y = 0, don't care for y = 0 or NaN
-            return 0.0;
-        using std::pow;
-        using std::abs;
-        return pow( abs(x) , y );
-    }
-}
-
-double signed_pow( double x , double k )
-{
-    using boost::math::sign;
-    return checked_math::pow( x , k ) * sign(x);
-}
-
-struct osc_chain {
-    const double m_kap, m_lam;
-    osc_chain( const double kap , const double lam )
-        : m_kap( kap ) , m_lam( lam ) { }
-
-    void operator()( const boost::numeric::odeint::mpi_state< std::vector<double> > &q ,
-                           boost::numeric::odeint::mpi_state< std::vector<double> > &dpdt ) const
-    {
-        const bool have_left = q.world.rank() > 0;
-        const bool have_right = q.world.rank() + 1 < q.world.size();
-        double q_left = 0, q_right = 0;
-        boost::mpi::request r_left, r_right;
-        if(have_left)
-        {
-            q.world.isend(q.world.rank() - 1, 0, q().front());
-            r_left = q.world.irecv(q.world.rank() - 1, 0, q_left);
-        }
-        if(have_right)
-        {
-            q.world.isend(q.world.rank() + 1, 0, q().back());
-            r_right = q.world.irecv(q.world.rank() + 1, 0, q_right);
-        }
-
-        double coupling_lr = 0;
-        if(have_left)
-        {
-            r_left.wait();
-            coupling_lr = signed_pow( q_left - q()[0] , m_lam-1 );
-        }
-        const size_t N = q().size();
-        for(size_t i = 0 ; i < N-1 ; ++i)
-        {
-            dpdt()[i] = -signed_pow( q()[i] , m_kap-1 ) + coupling_lr;
-            coupling_lr = signed_pow( q()[i] - q()[i+1] , m_lam-1 );
-            dpdt()[i] -= coupling_lr;
-        }
-        dpdt()[N-1] = -signed_pow( q()[N-1] , m_kap-1 ) + coupling_lr;
-        if(have_right)
-        {
-            r_right.wait();
-            dpdt()[N-1] -= signed_pow( q()[N-1] - q_right , m_lam-1 );
-        }
-    }
-    //]
-};
-
-#endif

performance/mpi/osc_chain_speedup.gnu

@@ -1,38 +0,0 @@
-#!/usr/bin/env gnuplot
-
-set terminal pngcairo size 1000,1000
-set output "osc_chain_speedup.png"
-
-set multiplot layout 2,2
-
-set key left
-
-set xrange [1:16]
-set x2range [1:16]
-set x2tics 8 format ""
-set grid x2tics
-set yrange [0:8]
-
-set title "short: speedup"
-plot \
-    "osc_chain_speedup-short.dat" i 0 u "block":"mul" w lp t "MPI" , \
-    (x < 4 ? x : 4) lc 0 lt 0 t "target"
-
-unset key
-
-set title "long: speedup"
-plot \
-    "osc_chain_speedup-long.dat" i 0 u "block":"mul" w lp, \
-    (x < 4 ? x : 4) lc 0 lt 0
-
-set yrange [0:*]
-
-set title "short: time[s]"
-plot \
-    "osc_chain_speedup-short.dat" i 0 u "block":"med" w lp
-
-set title "long: time[s]"
-plot \
-    "osc_chain_speedup-long.dat" i 0 u "block":"med" w lp
-
-unset multiplot

performance/mpi/osc_chain_speedup.sh

@@ -1,24 +0,0 @@
-#!/bin/bash
-
-bench="bin/clang-linux-3.2/release/threading-multi/osc_chain_1d"
-
-repeat=5
-maxnodes=16
-
-function run {
-    n=$1
-    steps=$2
-    for ((nodes=1 ; nodes < $maxnodes ; nodes++)) ; do
-        # swap stderr & stdout
-        mpirun -np $nodes $bench $n $steps $repeat 3>&1 1>&2 2>&3
-    done
-}
-
-function run_all {
-    printf "n\tsteps\tnodes\ttime\n"
-    run 256 1024
-    run 4096 1024
-    run 4194304 1
-}
-
-run_all | tee osc_chain_speedup.dat

performance/nr_rk4_lorenz.cpp

@@ -1,84 +0,0 @@
-/*
- * nr_rk4_lorenz.cpp
- *
- * 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)
- */
-
-
-#include <boost/array.hpp>
-
-#include "rk_performance_test_case.hpp"
-
-#include "lorenz.hpp"
-
-const size_t dim = 3;
-
-typedef boost::array< double , dim > state_type;
-
-
-template< class System , typename T , size_t dim >
-void rk4_step( const System sys , boost::array< T , dim > &x , const double t , const double dt )
-{   // fast rk4 implementation adapted from the book 'Numerical Recipes'
-    size_t i;
-    const double hh = dt*0.5;
-    const double h6 = dt/6.0;
-    const double th = t+hh;
-    boost::array< T , dim > dydx , dym , dyt , yt;
-
-    sys( x , dydx , t );
-
-    for( i=0 ; i<dim ; i++ )
-        yt[i] = x[i] + hh*dydx[i];
-
-    sys( yt , dyt , th );
-    for( i=0 ; i<dim ; i++ )
-        yt[i] = x[i] + hh*dyt[i];
-
-    sys( yt , dym , th );
-    for( i=0 ; i<dim ; i++ ) {
-        yt[i] = x[i] + dt*dym[i];
-        dym[i] += dyt[i];
-    }
-    sys( yt , dyt , t+dt );
-    for( i=0 ; i<dim ; i++ )
-        x[i] += h6*( dydx[i] + dyt[i] + 2.0*dym[i] );
-}
-
-
-class nr_wrapper
-{
-public:
-    void reset_init_cond()
-    {
-        m_x[0] = 10.0 * rand() / RAND_MAX;
-        m_x[1] = 10.0 * rand() / RAND_MAX;
-        m_x[2] = 10.0 * rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        rk4_step( lorenz() , m_x , m_t , dt );
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-};
-
-
-
-int main()
-{
-    nr_wrapper stepper;
-
-    run( stepper );
-}

performance/nr_rk4_phase_lattice.cpp

@@ -1,85 +0,0 @@
-/*
- * nr_rk4_phase_lattice.cpp
- *
- * 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)
- */
-
-
-#include <boost/array.hpp>
-
-#include "rk_performance_test_case.hpp"
-
-#include "phase_lattice.hpp"
-
-const size_t dim = 1024;
-
-typedef boost::array< double , dim > state_type;
-
-
-template< class System , typename T , size_t dim >
-void rk4_step( const System sys , boost::array< T , dim > &x , const double t , const double dt )
-{   // fast rk4 implementation adapted from the book 'Numerical Recipes'
-    size_t i;
-    const double hh = dt*0.5;
-    const double h6 = dt/6.0;
-    const double th = t+hh;
-    boost::array< T , dim > dydx , dym , dyt , yt;
-
-    sys( x , dydx , t );
-
-    for( i=0 ; i<dim ; i++ )
-        yt[i] = x[i] + hh*dydx[i];
-
-    sys( yt , dyt , th );
-    for( i=0 ; i<dim ; i++ )
-        yt[i] = x[i] + hh*dyt[i];
-
-    sys( yt , dym , th );
-    for( i=0 ; i<dim ; i++ ) {
-        yt[i] = x[i] + dt*dym[i];
-        dym[i] += dyt[i];
-    }
-    sys( yt , dyt , t+dt );
-    for( i=0 ; i<dim ; i++ )
-        x[i] += h6*( dydx[i] + dyt[i] + 2.0*dym[i] );
-}
-
-
-class nr_wrapper
-{
-public:
-    void reset_init_cond()
-    {
-        for( size_t i = 0 ; i<dim ; ++i )
-            m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        rk4_step( phase_lattice<dim>() , m_x , m_t , dt );
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-};
-
-
-
-int main()
-{
-    srand( 12312354 );
-
-    nr_wrapper stepper;
-
-    run( stepper , 10000 , 1E-6 );
-}

performance/odeint_rk4_array.cpp

@@ -0,0 +1,63 @@
+/*
+ * odeint_rk4_array
+ *
+ * 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)
+ */
+
+#include <iostream>
+
+#include <boost/timer.hpp>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/algebra/array_algebra.hpp>
+
+#include "lorenz.hpp"
+
+typedef boost::timer timer_type;
+
+typedef boost::array< double , 3 > state_type;
+
+using namespace boost::numeric::odeint;
+
+//typedef boost::numeric::odeint::runge_kutta4_classic< state_type > rk4_odeint_type;
+
+// use the never resizer explicitely for optimal performance with gcc, 
+// for the intel compiler this doesnt matter and the above definition
+// gives the same performance
+typedef runge_kutta4_classic< state_type , double , state_type , double ,
+                              array_algebra, default_operations, never_resizer > rk4_odeint_type;
+
+
+const int loops = 21;
+const int num_of_steps = 20000000;
+const double dt = 1E-10;
+
+
+int main()
+{
+    double min_time = 1E6; // something big
+    rk4_odeint_type stepper;
+    std::clog.precision(16);
+    std::cout.precision(16);
+    for( int n=0; n<loops; n++ )
+    {
+        state_type x = {{ 8.5, 3.1, 1.2 }};
+        double t = 0.0;
+        timer_type timer;
+        for( size_t i = 0 ; i < num_of_steps ; ++i )
+        {
+            stepper.do_step( lorenz(), x, t, dt );
+            t += dt;
+        }
+        min_time = std::min( timer.elapsed() , min_time );
+        std::clog << timer.elapsed() << '\t' << x[0] << std::endl;
+    }
+    std::cout << "Minimal Runtime: " << min_time << std::endl;
+}

performance/odeint_rk4_lorenz_array.cpp

@@ -1,71 +0,0 @@
-/*
- * odeint_rk4_lorenz.cpp
- *
- * Copyright 2011-2012 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 <boost/array.hpp>
-
-#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
-#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
-#include <boost/numeric/odeint/algebra/array_algebra.hpp>
-
-#include "rk_performance_test_case.hpp"
-
-#include "lorenz.hpp"
-
-typedef boost::array< double , 3 > state_type;
-typedef boost::numeric::odeint::runge_kutta4< state_type , double , state_type , double ,
-                                              boost::numeric::odeint::array_algebra > rk4_odeint_type;
-
-
-class odeint_wrapper
-{
-public:
-    void reset_init_cond()
-    {
-        /* random */
-        /*
-        m_x[0] = 10.0 * rand() / RAND_MAX;
-        m_x[1] = 10.0 * rand() / RAND_MAX;
-        m_x[2] = 10.0 * rand() / RAND_MAX;
-        */
-        /* hand chosen random (cf fortran) */
-        m_x[0] = 8.5;
-        m_x[1] = 3.1;
-        m_x[2] = 1.2;
-
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        m_stepper.do_step( lorenz() , m_x , m_t , dt );
-        //m_t += dt;
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-    rk4_odeint_type m_stepper;
-};
-
-
-
-int main()
-{
-    srand( 12312354 );
-
-    odeint_wrapper stepper;
-
-    run( stepper );
-}

performance/odeint_rk4_lorenz_range.cpp

@@ -1,59 +0,0 @@
-/*
- * odeint_rk4_lorenz_def_alg.cpp
- *
- * 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)
- */
-
-
-#include <boost/array.hpp>
-
-#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
-#include <boost/numeric/odeint/algebra/array_algebra.hpp>
-
-#include "rk_performance_test_case.hpp"
-
-#include "lorenz.hpp"
-
-typedef boost::array< double , 3 > state_type;
-typedef boost::numeric::odeint::runge_kutta4_classic< state_type > rk4_odeint_type;
-
-
-class odeint_wrapper
-{
-public:
-    void reset_init_cond()
-    {
-        m_x[0] = 10.0 * rand() / RAND_MAX;
-        m_x[1] = 10.0 * rand() / RAND_MAX;
-        m_x[2] = 10.0 * rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        m_stepper.do_step( lorenz() , m_x , m_t , dt );
-        //m_t += dt;
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-    rk4_odeint_type m_stepper;
-};
-
-
-
-int main()
-{
-    odeint_wrapper stepper;
-
-    run( stepper );
-}

performance/odeint_rk4_phase_lattice.cpp

@@ -1,64 +0,0 @@
-/*
- * odeint_rk4_phase_lattice.cpp
- *
- * 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)
- */
-
-
-#include <cmath>
-
-#include <boost/array.hpp>
-
-#include <boost/numeric/odeint/stepper/runge_kutta4_classic.hpp>
-#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
-#include <boost/numeric/odeint/algebra/array_algebra.hpp>
-
-#include "rk_performance_test_case.hpp"
-
-#include "phase_lattice.hpp"
-
-const size_t N = 1024;
-
-typedef boost::array< double , N > state_type;
-typedef boost::numeric::odeint::runge_kutta4_classic< state_type , double , state_type , double , boost::numeric::odeint::array_algebra> rk4_odeint_type;
-
-class odeint_wrapper
-{
-public:
-    void reset_init_cond()
-    {
-        for( size_t i = 0 ; i<N ; ++i )
-            m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        m_stepper.do_step( phase_lattice<N>() , m_x , m_t , dt );
-        //m_t += dt;
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-    
-private:
-    state_type m_x;
-    double m_t;
-    rk4_odeint_type m_stepper;
-};
-
-
-
-int main()
-{
-    srand( 12312354 );
-
-    odeint_wrapper stepper;
-
-    run( stepper , 10000 , 1E-6 );
-}

performance/odeint_rk4_phase_lattice_mkl.cpp

@@ -1,69 +0,0 @@
-/*
- * odeint_rk4_phase_lattice_mkl.cpp
- *
- * 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)
- */
-
-
-#include <boost/array.hpp>
-
-#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
-#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
-#include <boost/numeric/odeint/external/mkl/mkl_operations.hpp>
-
-
-#include "rk_performance_test_case.hpp"
-
-#include "phase_lattice.hpp"
-#include "phase_lattice_mkl.hpp"
-
-const size_t N = 1024;
-
-typedef boost::array< double , N > state_type;
-typedef boost::numeric::odeint::runge_kutta4< state_type 
-                                              , double , state_type , double ,
-                                              boost::numeric::odeint::vector_space_algebra , 
-                                              boost::numeric::odeint::mkl_operations 
-                                              > rk4_odeint_type;
-
-
-class odeint_wrapper
-{
-public:
-    void reset_init_cond()
-    {
-        for( size_t i = 0 ; i<N ; ++i )
-            m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        m_stepper.do_step( phase_lattice_mkl<N>() , m_x , m_t , dt );
-        //m_t += dt;
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-    rk4_odeint_type m_stepper;
-};
-
-
-
-int main()
-{
-    srand( 12312354 );
-
-    odeint_wrapper stepper;
-
-    run( stepper , 10000 , 1E-6 );
-}

performance/openmp/Jamfile.v2

@@ -1,21 +0,0 @@
-# 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>../../../../..
-      <include>..
-      <define>BOOST_ALL_NO_LIB=1
-      <library>/boost//timer
-      <library>/boost//program_options
-      [ openmp ]
-    ;
-
-exe osc_chain_1d : osc_chain_1d.cpp ;

performance/openmp/osc_chain_1d.cpp

@@ -1,127 +0,0 @@
-/* Boost libs/numeric/odeint/performance/openmp/osc_chain_1d.cpp
-
- Copyright 2013 Karsten Ahnert
- Copyright 2013 Pascal Germroth
- Copyright 2013 Mario Mulansky
-
- stronlgy nonlinear hamiltonian lattice in 2d
-
- 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 <omp.h>
-#include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/external/openmp/openmp.hpp>
-
-#include <boost/program_options.hpp>
-#include <boost/random.hpp>
-#include <boost/timer/timer.hpp>
-#include <boost/foreach.hpp>
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics/stats.hpp>
-#include <boost/accumulators/statistics/mean.hpp>
-#include <boost/accumulators/statistics/median.hpp>
-#include "osc_chain_1d_system.hpp"
-
-using namespace std;
-using namespace boost::numeric::odeint;
-using namespace boost::accumulators;
-using namespace boost::program_options;
-
-using boost::timer::cpu_timer;
-
-const double p_kappa = 3.3;
-const double p_lambda = 4.7;
-
-int main( int argc , char* argv[] )
-{
-    size_t N, blocks, steps, repeat;
-    bool split_range, dump;
-    options_description desc("Options");
-    desc.add_options()
-        ("help,h", "show this help")
-        ("length", value(&N)->default_value(1024), "length of chain")
-        ("steps", value(&steps)->default_value(100), "simulation steps")
-        ("blocks", value(&blocks)->default_value(omp_get_max_threads()), "number of blocks (split) or threads (non-split)")
-        ("split", bool_switch(&split_range), "split range")
-        ("repeat", value(&repeat)->default_value(25), "repeat runs")
-        ("dump", bool_switch(&dump), "dump final state to stderr")
-        ;
-    variables_map vm;
-    store(command_line_parser(argc, argv).options(desc).run(), vm);
-    notify(vm);
-    if(vm.count("help"))
-    {
-        cerr << desc << endl;
-        return EXIT_FAILURE;
-    }
-    cout << "length\tsteps\tthreads\ttime" << endl;
-
-    accumulator_set< double, stats<tag::mean, tag::median> > acc_time;
-
-    vector<double> p( N ), q( N, 0 );
-    boost::random::uniform_real_distribution<double> distribution;
-    boost::random::mt19937 engine( 0 );
-    generate( p.begin() , p.end() , boost::bind( distribution , engine ) );
-
-    if(split_range) {
-        typedef openmp_state<double> state_type;
-        typedef symplectic_rkn_sb3a_mclachlan<
-                  state_type , state_type , double
-                > stepper_type;
-        state_type p_split(blocks), q_split(blocks);
-        split(p, p_split);
-        split(q, q_split);
-
-        for(size_t n_run = 0 ; n_run != repeat ; n_run++) {
-            cpu_timer timer;
-            integrate_n_steps( stepper_type() , osc_chain( p_kappa , p_lambda ) ,
-                               make_pair( boost::ref(q_split) , boost::ref(p_split) ) ,
-                               0.0 , 0.01 , steps );
-            double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
-            acc_time(run_time);
-            cout << N << '\t' << steps << '\t' << blocks << '\t' << run_time << endl;
-        }
-
-        if(dump) {
-            unsplit(p_split, p);
-            copy(p.begin(), p.end(), ostream_iterator<double>(cerr, "\t"));
-            cerr << endl;
-        }
-
-    } else {
-        typedef vector<double> state_type;
-        typedef symplectic_rkn_sb3a_mclachlan<
-                  state_type , state_type , double ,
-                  state_type , state_type , double ,
-                  openmp_range_algebra
-                > stepper_type;
-        omp_set_num_threads(blocks);
-
-        for(size_t n_run = 0 ; n_run != repeat ; n_run++) {
-            cpu_timer timer;
-            integrate_n_steps( stepper_type() , osc_chain( p_kappa , p_lambda ) ,
-                               make_pair( boost::ref(q) , boost::ref(p) ) ,
-                               0.0 , 0.01 , steps );
-            double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
-            acc_time(run_time);
-            cout << N << '\t' << steps << '\t' << blocks << '\t' << run_time << endl;
-        }
-
-        if(dump) {
-            copy(p.begin(), p.end(), ostream_iterator<double>(cerr, "\t"));
-            cerr << endl;
-        }
-
-    }
-
-    cout << "# mean=" << mean(acc_time)
-         << " median=" << median(acc_time) << endl;
-
-    return 0;
-}

performance/openmp/osc_chain_1d_system.hpp

@@ -1,98 +0,0 @@
-/* Boost libs/numeric/odeint/performance/openmp/osc_chain_1d_system.hpp
-
- Copyright 2013 Karsten Ahnert
- Copyright 2013 Mario Mulansky
- Copyright 2013 Pascal Germroth
-
- stronlgy nonlinear hamiltonian lattice
-
- 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 SYSTEM_HPP
-#define SYSTEM_HPP
-
-#include <vector>
-#include <cmath>
-#include <iostream>
-
-#include <omp.h>
-
-#include <boost/math/special_functions/sign.hpp>
-#include <boost/numeric/odeint/external/openmp/openmp.hpp>
-
-namespace checked_math {
-    inline double pow( double x , double y )
-    {
-        if( x==0.0 )
-            // 0**y = 0, don't care for y = 0 or NaN
-            return 0.0;
-        using std::pow;
-        using std::abs;
-        return pow( abs(x) , y );
-    }
-}
-
-double signed_pow( double x , double k )
-{
-    using boost::math::sign;
-    return checked_math::pow( x , k ) * sign(x);
-}
-
-struct osc_chain {
-
-    const double m_kap, m_lam;
-
-    osc_chain( const double kap , const double lam )
-        : m_kap( kap ) , m_lam( lam )
-    { }
-
-    // Simple case with openmp_range_algebra
-    void operator()( const std::vector<double> &q ,
-                           std::vector<double> &dpdt ) const
-    {
-        const size_t N = q.size();
-        double coupling_lr = 0;
-        size_t last_i = N;
-        #pragma omp parallel for firstprivate(coupling_lr, last_i) lastprivate(coupling_lr) schedule(runtime)
-        for(size_t i = 0 ; i < N - 1 ; ++i)
-        {
-            if(i > 0 && i != last_i + 1)
-                coupling_lr = signed_pow( q[i-1]-q[i] , m_lam-1 );
-            dpdt[i] = -signed_pow( q[i] , m_kap-1 ) + coupling_lr;
-            coupling_lr = signed_pow( q[i] - q[i+1] , m_lam-1 );
-            dpdt[i] -= coupling_lr;
-            last_i = i;
-        }
-        dpdt[N-1] = -signed_pow( q[N-1] , m_kap-1 ) + coupling_lr;
-    }
-
-    // Split case with openmp_algebra
-    void operator()( const boost::numeric::odeint::openmp_state<double> &q ,
-                           boost::numeric::odeint::openmp_state<double> &dpdt ) const
-    {
-        const size_t M = q.size();
-        #pragma omp parallel for schedule(runtime)
-        for(size_t i = 0 ; i < M ; ++i)
-        {
-            const std::vector<double> &_q = q[i];
-            std::vector<double> &_dpdt = dpdt[i];
-            const size_t N = q[i].size();
-            double coupling_lr = 0;
-            if(i > 0) coupling_lr = signed_pow( q[i-1].back() - _q[0] , m_lam-1 );
-            for(size_t j = 0 ; j < N-1 ; ++j)
-            {
-                _dpdt[j] = -signed_pow( _q[j] , m_kap-1 ) + coupling_lr;
-                coupling_lr = signed_pow( _q[j] - _q[j+1] , m_lam-1 );
-                _dpdt[j] -= coupling_lr;
-            }
-            _dpdt[N-1] = -signed_pow( _q[N-1] , m_kap-1 ) + coupling_lr;
-            if(i + 1 < M) _dpdt[N-1] -= signed_pow( _q[N-1] - q[i+1].front() , m_lam-1 );
-        }
-    }
-
-};
-
-#endif

performance/openmp/osc_chain_speedup.gnu

@@ -1,50 +0,0 @@
-#!/usr/bin/env gnuplot
-
-set terminal pngcairo size 1000,1000
-set output "osc_chain_speedup.png"
-
-set multiplot layout 2,2
-
-set key left
-
-set xrange [1:64]
-set x2range [1:64]
-set x2tics 8 format ""
-set grid x2tics
-set yrange [0:8]
-
-set title "short: speedup"
-plot \
-    "osc_chain_speedup-short.dat" i 0 u "block":"gcc-s-mul" w lp t "gcc (split)" , \
-    "osc_chain_speedup-short.dat" i 0 u "block":"gcc-t-mul" w lp t "gcc (simple)", \
-    "osc_chain_speedup-short.dat" i 0 u "block":"icc-s-mul" w lp t "icc (split)" , \
-    "osc_chain_speedup-short.dat" i 0 u "block":"icc-t-mul" w lp t "icc (simple)", \
-    (x < 4 ? x : 4) lc 0 lt 0 t "target"
-
-unset key
-
-set title "long: speedup"
-plot \
-    "osc_chain_speedup-long.dat" i 0 u "block":"gcc-s-mul" w lp, \
-    "osc_chain_speedup-long.dat" i 0 u "block":"gcc-t-mul" w lp, \
-    "osc_chain_speedup-long.dat" i 0 u "block":"icc-s-mul" w lp, \
-    "osc_chain_speedup-long.dat" i 0 u "block":"icc-t-mul" w lp, \
-    (x < 4 ? x : 4) lc 0 lt 0
-
-set yrange [0:*]
-
-set title "short: time[s]"
-plot \
-    "osc_chain_speedup-short.dat" i 0 u "block":"gcc-s-med" w lp, \
-    "osc_chain_speedup-short.dat" i 0 u "block":"gcc-t-med" w lp, \
-    "osc_chain_speedup-short.dat" i 0 u "block":"icc-s-med" w lp, \
-    "osc_chain_speedup-short.dat" i 0 u "block":"icc-t-med" w lp
-
-set title "long: time[s]"
-plot \
-    "osc_chain_speedup-long.dat" i 0 u "block":"gcc-s-med" w lp, \
-    "osc_chain_speedup-long.dat" i 0 u "block":"gcc-t-med" w lp, \
-    "osc_chain_speedup-long.dat" i 0 u "block":"icc-s-med" w lp, \
-    "osc_chain_speedup-long.dat" i 0 u "block":"icc-t-med" w lp
-
-unset multiplot

performance/openmp/osc_chain_speedup.sh

@@ -1,38 +0,0 @@
-#!/bin/zsh
-
-export LC_NUMERIC=en_US.UTF-8
-declare -A times
-
-export OMP_SCHEDULE=static
-export OMP_PROC_BIND=true
-repeat=2
-
-function run {
-    n=$1
-    steps=$2
-    printf "# n=$n steps=$steps repeat=$repeat\n"
-    printf '"block"'
-    for b in gcc icc ; do
-        for s in s t ; do
-            for t in med mul ; do
-                printf "\t\"$b-$s-$t\""
-            done
-        done
-    done
-    for block in 1 2 4 8 16 32 64; do
-        printf '\n%d' $block
-        for build in gcc-4.7 intel-linux ; do
-            bench="bin/$build/release/osc_chain_1d"
-            for split in 1 0 ; do
-                med=$($bench $n $block $steps $repeat $split | tail -1 | awk '{print $4}')
-                times[$build-$split-$block]=$med
-                speedup=$((${times[$build-$split-1]}/$med))
-                printf '\t%f\t%f' $med $speedup
-            done
-        done
-    done
-    printf '\n\n\n'
-}
-
-run 4096 1024 | tee osc_chain_speedup-short.dat
-run 524288 10 | tee osc_chain_speedup-long.dat

performance/performance.py

@@ -1,70 +0,0 @@
-"""
- 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)
-"""
-
-
-from os import popen
-from os import system
-from os.path import isfile
-from numpy import *
-#from pylab import *
-
-#toolset = "gcc-4.5"
-#toolset = "intel-11.1"
-toolset = "msvc"
-#toolset = "msvc-10.0"
-
-#bin_path = "bin/gcc-4.5/release/"
-#bin_path = "bin/intel-linux-11.1/release/"
-bin_path = "bin\\msvc-10.0\\release\\threading-multi\\"
-extension = ".exe"
-#extension = ""
-
-bins = [ "odeint_rk4_lorenz_array" , "odeint_rk4_lorenz_range" , "generic_odeint_rk4_lorenz" , "nr_rk4_lorenz" ,  "rt_generic_rk4_lorenz" , "gsl_rk4_lorenz" ]
-
-results = []
-
-print "Performance tests for " , bin_path
-print
-
-for bin in bins:
-        #system( "bjam toolset=" + toolset + " -a " + bin );
-        if isfile( bin_path + bin + extension):
-                print "Running" , bin
-                res = popen( bin_path+bin+extension ).read()
-                print bin , res
-                results.append( res )
-        else:
-                print "no executable found:" , bin_path + bin + extension
-                results.append( 0 )
-
-print "Results from" , bin_path
-print
-
-for i in range(len(bins)):
-        print bins[i] , results[i]
-
-res = array( results , dtype='float' )
-savetxt( bin_path + "rk4_lorenz.dat" , res )
-
-res = 100*res[0]/res
-
-bar_width = 0.6
-
-"""
-figure(1)
-title("Runge-Kutta 4 with " + toolset , fontsize=20)
-bar( arange(6) , res , bar_width , color='blue' , linewidth=4 , edgecolor='blue' , ecolor='red') #, elinewidth=2, ecolor='red' )
-xlim( -0.5 , 5.5+bar_width )
-xticks( arange(6)+bar_width/2 , ('array' , 'range' , 'generic' , 'NR' , 'rt gen' , 'gsl' ) )
-ylabel('Performance in %' , fontsize=20)
-
-savefig( bin_path + "rk4_lorenz.png" )
-
-show()
-"""

performance/phase_lattice.hpp

@@ -1,43 +0,0 @@
-/*
- * 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)
- */
-
-#include <cmath>
-
-#include <boost/array.hpp>
-
-template< size_t N >
-struct phase_lattice
-{
-    typedef double value_type;
-    typedef boost::array< value_type , N > state_type;
-
-    value_type m_epsilon; 
-    state_type m_omega;
-
-    phase_lattice() : m_epsilon( 6.0/(N*N) ) // should be < 8/N^2 to see phase locking
-    {
-        for( size_t i=1 ; i<N-1 ; ++i )
-            m_omega[i] = m_epsilon*(N-i);
-    }
-
-    void inline operator()( const state_type &x , state_type &dxdt , const double t ) const
-    {
-        double c = 0.0;
-
-        for( size_t i=0 ; i<N-1 ; ++i )
-        {
-            dxdt[i] = m_omega[i] + c;
-            c = ( x[i+1] - x[i] );
-            dxdt[i] += c;
-        }
-
-        //dxdt[N-1] = m_omega[N-1] + sin( x[N-1] - x[N-2] );
-    }
-
-};

performance/phase_lattice_mkl.hpp

@@ -1,57 +0,0 @@
-/*
- * phase_lattice_mkl.hpp
- *
- * 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)
- */
-
-
-
-
-#ifndef PHASE_LATTICE_MKL_HPP_
-#define PHASE_LATTICE_MKL_HPP_
-
-#include <cmath>
-
-#include <mkl_blas.h>
-#include <mkl_vml_functions.h>
-#include <boost/array.hpp>
-
-template< size_t N >
-struct phase_lattice_mkl
-{
-    typedef double value_type;
-    typedef boost::array< value_type , N > state_type;
-
-    value_type m_epsilon;
-    state_type m_omega;
-    state_type m_tmp;
-
-    phase_lattice_mkl() : m_epsilon( 6.0/(N*N) ) // should be < 8/N^2 to see phase locking
-    {
-        for( size_t i=1 ; i<N-1 ; ++i )
-            m_omega[i] = m_epsilon*(N-i);
-    }
-
-    void inline operator()( const state_type &x , state_type &dxdt , const double t )
-    {
-        const int n = x.size();
-
-        dxdt[0] = m_omega[0] + sin( x[1] - x[0] );
-
-        vdSub( n-1 , &(x[1]) , &(x[0]) , &(m_tmp[0]) );
-        vdSin( n-1 , &(m_tmp[0]) , &(m_tmp[0]) );
-        vdAdd( n-2 , &(m_tmp[0]) , &(m_tmp[1]) , &(dxdt[1]) );
-        vdAdd( n-2 , &(dxdt[1]) , &(m_omega[1]) , &(dxdt[1]) );
-
-        dxdt[N-1] = m_omega[N-1] + sin( x[N-1] - x[N-2] );
-    }
-
-};
-
-
-#endif /* PHASE_LATTICE_MKL_HPP_ */

performance/plot_result.py

@@ -1,40 +1,64 @@
 """
- Copyright 2011 Mario Mulansky
- Copyright 2012 Karsten Ahnert
- 
+ Copyright 2011-2014 Mario Mulansky
+ Copyright 2011-2014 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)
 """
 
+import numpy as np
+from matplotlib import pyplot as plt
 
-from pylab import *
+plt.rc("font", size=16)
 
-#toolset = "gcc-4.5"
-toolset = "intel-11.1"
-#toolset = "msvc"
-#toolset = "msvc-10.0"
 
-#bin_path = "bin/gcc-4.5/release/"
-bin_path = "bin/intel-linux-11.1/release/"
-#bin_path = "bin\\msvc-10.0\\release\\" #threading-multi\\"
-#extension = ".exe"
-extension = ""
+def get_runtime_from_file(filename):
+    gcc_perf_file = open(filename, 'r')
+    for line in gcc_perf_file:
+        if "Minimal Runtime:" in line:
+            return float(line.split(":")[-1])
 
-res = loadtxt( bin_path + "rk4_lorenz.dat" )
 
-res = 100*res[0]/res
+t_gcc = [get_runtime_from_file("perf_workbook/odeint_rk4_array_gcc.perf"),
+         get_runtime_from_file("perf_ariel/odeint_rk4_array_gcc.perf"),
+         get_runtime_from_file("perf_lyra/odeint_rk4_array_gcc.perf")]
 
-bar_width = 0.6
+t_intel = [get_runtime_from_file("perf_workbook/odeint_rk4_array_intel.perf"),
+           get_runtime_from_file("perf_ariel/odeint_rk4_array_intel.perf"),
+           get_runtime_from_file("perf_lyra/odeint_rk4_array_intel.perf")]
 
-figure(1)
-title("Runge-Kutta 4 with " + toolset , fontsize=20)
-bar( arange(6) , res , bar_width , color='blue' , linewidth=4 , edgecolor='blue' , ecolor='red') #, elinewidth=2, ecolor='red' )
-xlim( -0.5 , 5.5+bar_width )
-ylim( 0 , max( res ) + 10 )
-xticks( arange(6)+bar_width/2 , ('array' , 'range' , 'generic' , 'NR' , 'rt gen' , 'gsl' ) )
-ylabel('Performance in %' , fontsize=20)
+t_gfort = [get_runtime_from_file("perf_workbook/rk4_gfort.perf"),
+           get_runtime_from_file("perf_ariel/rk4_gfort.perf"),
+           get_runtime_from_file("perf_lyra/rk4_gfort.perf")]
 
-savefig( bin_path + "rk4_lorenz.png" )
+t_c_intel = [get_runtime_from_file("perf_workbook/rk4_c_intel.perf"),
+             get_runtime_from_file("perf_ariel/rk4_c_intel.perf"),
+             get_runtime_from_file("perf_lyra/rk4_c_intel.perf")]
 
-show()
+print t_c_intel
+
+
+ind = np.arange(3)  # the x locations for the groups
+width = 0.15         # the width of the bars
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+rects1 = ax.bar(ind, t_gcc, width, color='b', label="odeint gcc")
+rects2 = ax.bar(ind+width, t_intel, width, color='g', label="odeint intel")
+rects3 = ax.bar(ind+2*width, t_c_intel, width, color='y', label="C intel")
+rects4 = ax.bar(ind+3*width, t_gfort, width, color='c', label="gfort")
+
+ax.axis([-width, 2.0+5*width, 0.0, 0.85])
+ax.set_ylabel('Runtime (s)')
+ax.set_title('Performance for integrating the Lorenz system')
+ax.set_xticks(ind + 1.5*width)
+ax.set_xticklabels(('Core i5-3210M\n3.1 GHz',
+                    'Xeon E5-2690\n3.8 GHz',
+                    'Opteron 8431\n 2.4 GHz'))
+ax.legend(loc='upper left', prop={'size': 16})
+
+plt.savefig("perf.pdf")
+plt.savefig("perf.png", dpi=50)
+
+plt.show()

performance/rk4_lorenz.f

@@ -1,53 +0,0 @@
-C
-C  NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, (c) 1985
-C
-C FILE: rk4sys.f
-C
-C RUNGE-KUTTA METHOD OF ORDER 4 FOR A SYSTEM OF ODE'S (RK4SYS,XPSYS)
-C
-      DOUBLE PRECISION X,T,H
-      DIMENSION  X(3) 
-      DATA  N/3/, T/0.0/, X/8.5,3.1,1.2/
-      DATA  H/1E-10/, NSTEP/20000000/ 
-      CALL RK4SYS(N,T,X,H,NSTEP)
-      STOP
-      END 
-  
-      SUBROUTINE XPSYS(X,F)
-      DOUBLE PRECISION X,F
-      DIMENSION  X(3),F(3)  
-      F(1) = 10.0 * ( X(2) - X(1) )
-      F(2) = 28.0 * X(1) - X(2) - X(1) * X(3)
-      F(3) = X(1)*X(2) - (8.0/3.0) * X(3)
-      RETURN
-      END 
-
-      SUBROUTINE RK4SYS(N,T,X,H,NSTEP)
-      DOUBLE PRECISION X,Y,T,H,F1,F2,F3,F4
-      DIMENSION  X(N),Y(N),F1(N),F2(N),F3(N),F4(N)  
-      PRINT 7,T,(X(I),I=1,N)
-      H2 = 0.5*H  
-      START = T   
-      DO 6 K = 1,NSTEP      
-        CALL XPSYS(X,F1)    
-        DO 2 I = 1,N
-          Y(I) = X(I) + H2*F1(I)      
-   2    CONTINUE
-        CALL XPSYS(Y,F2)    
-        DO 3 I = 1,N
-          Y(I) = X(I) + H2*F2(I)      
-   3    CONTINUE
-        CALL XPSYS(Y,F3)    
-        DO 4 I = 1,N
-          Y(I) = X(I) + H*F3(I)       
-   4    CONTINUE
-        CALL XPSYS(Y,F4)    
-        DO 5 I = 1,N
-          X(I) = X(I) + H*(F1(I) + 2.0*(F2(I) + F3(I)) + F4(I))/6.0 
-   5    CONTINUE
-   6  CONTINUE
-      PRINT 7,T,(X(I),I = 1,N)
-   7  FORMAT(2X,'T,X:',E10.3,5(2X,E22.14))
-      RETURN      
-      END 
-

performance/rk_performance_test_case.hpp

@@ -1,68 +0,0 @@
-/*
- * rk_performance_test_case.hpp
- *
- * Copyright 2011-2012 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/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics.hpp>
-#include <boost/timer.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace boost::accumulators;
-
-typedef accumulator_set<
-    double , stats< tag::mean , tag::variance >
-    > accumulator_type;
-
-ostream& operator<<( ostream& out , accumulator_type &acc )
-{
-    out << boost::accumulators::mean( acc ) << tab;
-//    out << boost::accumulators::variance( acc ) << tab;
-    return out;
-}
-
-typedef boost::timer timer_type;
-
-
-template< class Stepper >
-void run( Stepper &stepper , const size_t num_of_steps = 20000000 , const double dt = 1E-10 )
-{
-    const size_t loops = 20;
-
-    accumulator_type acc;
-    timer_type timer;
-
-    srand( 12312354 );
-
-    // transient
-    //stepper.reset_init_cond( );
-    //for( size_t i = 0 ; i < num_of_steps ; ++i )
-    //    stepper.do_step( dt );
-
-    for( size_t n=0 ; n<loops+1 ; ++n )
-    {
-        stepper.reset_init_cond( );
-
-        timer.restart();
-        for( size_t i = 0 ; i < num_of_steps ; ++i )
-            stepper.do_step( dt );
-        if( n>0 )
-        {   // take first run as transient
-            acc(timer.elapsed());
-            clog.precision(8);
-            clog.width(10);
-            clog << acc << " " << stepper.state(0) << endl;
-        }
-    }
-    cout << acc << endl;
-}

performance/rt_algebra.hpp

@@ -1,32 +0,0 @@
-/*
- * 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)
- */
-
-
-#include <vector>
-
-using namespace std;
-
-struct rt_algebra
-{
-    template< typename T , size_t dim >
-    inline static void foreach( boost::array< T , dim > & x_tmp , 
-                                const boost::array< T , dim > &x ,
-                                //const vector< double > &a ,
-                                const double* a ,
-                                const boost::array< T , dim > *k_vector , 
-                                const double dt , const size_t s )
-    {
-        for( size_t i=0 ; i<dim ; ++i )
-        {
-            x_tmp[i] = x[i];
-            for( size_t j = 0 ; j<s ; ++j )
-                x_tmp[i] += a[j]*dt*k_vector[j][i];
-        }
-    }
-};

performance/rt_explicit_rk.hpp

@@ -1,87 +0,0 @@
-/*
- * rt_explicit_rk.hpp
- *
- * 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)
- */
-
-
-#ifndef RT_EXPLICIT_RK_HPP_
-#define RT_EXPLICIT_RK_HPP_
-
-#include <vector>
-
-#include "rt_algebra.hpp"
-
-using namespace std;
-
-template< class StateType >
-class rt_explicit_rk
-{
-public:
-    typedef StateType state_type;
-    typedef double* const * coeff_a_type;
-    typedef vector< double > coeff_b_type;
-    typedef vector< double > coeff_c_type;
-
-    rt_explicit_rk( size_t stage_count ) : m_s( stage_count )
-    {
-        m_F = new state_type[ m_s ];
-    }
-
-    rt_explicit_rk( const size_t stage_count , 
-                    const coeff_a_type a , 
-                    const coeff_b_type &b , const coeff_c_type &c )
-        : m_s( stage_count ) , m_a( a ) , m_b( b ) , m_c( c )
-    {
-        m_F = new state_type[ m_s ];
-    }
-
-    ~rt_explicit_rk()
-    {
-        delete[] m_F;
-    }
-
-    /* void set_params( coeff_a_type &a , coeff_b_type &b , coeff_c_type &c )
-    {
-        m_a = a;
-        m_b = b;
-        m_c = c;
-        }*/
-
-    template< class System >
-    void do_step( System sys , state_type &x , const double t , const double dt )
-    {
-        // first stage separately
-        sys( x , m_F[0] , t + m_c[0]*t );
-        if( m_s == 1 )
-            rt_algebra::foreach( x , x , &m_b[0] , m_F , dt , 1 );
-        else
-            rt_algebra::foreach( m_x_tmp , x , m_a[0] , m_F , dt , 1 );
-
-        for( size_t stage = 2 ; stage <= m_s ; ++stage )
-        {
-            sys( m_x_tmp , m_F[stage-1] , t + m_c[stage-1]*dt );
-            if( stage == m_s )
-                rt_algebra::foreach( x , x , &m_b[0] , m_F , dt , stage-1 );
-            else
-                rt_algebra::foreach( m_x_tmp , x , m_a[stage-1] , m_F , dt , stage-1 );
-        }
-    }
-
-
-private:
-    const size_t m_s;
-    const coeff_a_type m_a;
-    const coeff_b_type m_b;
-    const coeff_c_type m_c;
-
-    state_type m_x_tmp;
-    state_type *m_F;
-};
-
-#endif /* RT_EXPLICIT_RK_HPP_ */

performance/rt_generic_rk4_lorenz.cpp

@@ -1,81 +0,0 @@
-/*
- * rt_generic_rk4_lorenz.cpp
- *
- * 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)
- */
-
-
-#include <boost/array.hpp>
-
-#include "rt_explicit_rk.hpp"
-
-#include "rk_performance_test_case.hpp"
-
-#include "lorenz.hpp"
-
-typedef boost::array< double , 3 > state_type;
-
-typedef rt_explicit_rk< state_type > rk_stepper_type;
-
-const size_t stage_count = 4;
-
-
-class rt_generic_wrapper
-{
-public:
-
-    rt_generic_wrapper( const double * const * a , 
-                        const rk_stepper_type::coeff_b_type &b ,
-                        const rk_stepper_type::coeff_c_type &c ) 
-        : m_stepper( stage_count , 
-                     (rk_stepper_type::coeff_a_type) a , b , c )
-    { }
-
-    void reset_init_cond()
-    {
-        m_x[0] = 10.0 * rand() / RAND_MAX;
-        m_x[1] = 10.0 * rand() / RAND_MAX;
-        m_x[2] = 10.0 * rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        m_stepper.do_step( lorenz() , m_x , m_t , dt );
-        //m_t += dt;
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-    rk_stepper_type m_stepper;
-};
-
-
-
-int main()
-{
-
-    const double a_tmp[3*4/2] = { 0.5 ,
-                                  0.0 , 1.0 ,
-                                  0.0 , 0.0 , 1.0 };
-    const double* const a[3] = { a_tmp , a_tmp+1 , a_tmp+3 };
-
-    rk_stepper_type::coeff_b_type b( stage_count );
-    b[0] = 1.0/6; b[1] = 1.0/3; b[2] = 1.0/3; b[3] = 1.0/6;
-    
-    rk_stepper_type::coeff_c_type c( stage_count );
-    c[0] = 0.0; c[1] = 0.5; c[2] = 0.5; c[3] = 1.0;
-
-    rt_generic_wrapper stepper( a , b , c );
-
-    run( stepper );
-}

performance/rt_generic_rk4_phase_lattice.cpp

@@ -1,83 +0,0 @@
-/*
- * rt_generic_rk4_lorenz.cpp
- *
- * 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)
- */
-
-
-#include <boost/array.hpp>
-
-#include "rt_explicit_rk.hpp"
-
-#include "rk_performance_test_case.hpp"
-
-#include "phase_lattice.hpp"
-
-const size_t N = 1024;
-
-typedef boost::array< double , N > state_type;
-
-typedef rt_explicit_rk< state_type > rk_stepper_type;
-
-const size_t stage_count = 4;
-
-
-class rt_generic_wrapper
-{
-public:
-
-    rt_generic_wrapper( const double * const * a , 
-                        const rk_stepper_type::coeff_b_type &b ,
-                        const rk_stepper_type::coeff_c_type &c ) 
-        : m_stepper( stage_count , 
-                     (rk_stepper_type::coeff_a_type) a , b , c )
-    { }
-
-    void reset_init_cond()
-    {
-        for( size_t i = 0 ; i<N ; ++i )
-            m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
-        m_t = 0.0;
-    }
-
-    inline void do_step( const double dt )
-    {
-        m_stepper.do_step( phase_lattice<N>() , m_x , m_t , dt );
-        //m_t += dt;
-    }
-
-    double state( const size_t i ) const
-    { return m_x[i]; }
-
-private:
-    state_type m_x;
-    double m_t;
-    rk_stepper_type m_stepper;
-};
-
-
-
-int main()
-{
-    srand( 12312354 );
-
-    const double a_tmp[3*4/2] = { 0.5 ,
-                                  0.0 , 1.0 ,
-                                  0.0 , 0.0 , 1.0 };
-    const double* const a[3] = { a_tmp , a_tmp+1 , a_tmp+3 };
-
-    rk_stepper_type::coeff_b_type b( stage_count );
-    b[0] = 1.0/6; b[1] = 1.0/3; b[2] = 1.0/3; b[3] = 1.0/6;
-    
-    rk_stepper_type::coeff_c_type c( stage_count );
-    c[0] = 0.0; c[1] = 0.5; c[2] = 0.5; c[3] = 1.0;
-
-    rt_generic_wrapper stepper( a , b , c );
-
-    run( stepper , 10000 , 1E-6 );
-}

test/Jamfile.v2

@@ -6,18 +6,25 @@
 
 # bring in rules for testing
 
-
 import testing ;
 
+# make sure you are using a new version of boost.build, otherwise the local
+# odeint will not be included properly
+# you can fix older boost.build versions by applying the patch provided in
+# odeint's root, e.g.:
+# cd ~/odeint-v2
+# sudo patch /usr/share/boost-build/build/toolset.jam toolset.jam.patch
+
 use-project boost : $(BOOST_ROOT) ;
 
 project
     : requirements
-      <library>/boost/test//boost_unit_test_framework
       <define>BOOST_ALL_NO_LIB=1
+      # use test library
+      <library>/boost//unit_test_framework
       <link>static
       <toolset>clang:<cxxflags>-Wno-unused-variable
-# <cxxflags>-D_SCL_SECURE_NO_WARNINGS
+      # <cxxflags>-D_SCL_SECURE_NO_WARNINGS
     ;
 
 test-suite "odeint"
@@ -74,12 +81,16 @@ test-suite "odeint"
     [ compile unwrap_boost_reference.cpp ]
     [ compile unwrap_reference.cpp : <cxxflags>-std=c++0x : unwrap_reference_C++11 ]
     [ compile-fail unwrap_reference.cpp : <cxxflags>-std=c++98 : unwrap_reference_C++98 ]
-    : <testing.launcher>valgrind
+    [ compile std_array.cpp : <cxxflags>-std=c++0x ]
+    : 
+        <testing.launcher>valgrind
     ;
     
 # also run numeric tests
 build-project numeric ;
 
+build-project regression ;
+
 # test-suite "odeint-iterator_integrate"
 #    :
 #    [ run integrate.cpp : : : : integrate_iterator ]

test/adams_bashforth.cpp

@@ -69,14 +69,14 @@ public:
     size_t do_count;
 
     template< class System , class StateIn , class DerivIn , class StateOut >
-    void do_step( System system , const StateIn &in , const DerivIn &dxdt , value_type t , StateOut &out , value_type dt )
+    void do_step_dxdt_impl( System system , const StateIn &in , const DerivIn &dxdt , value_type t , StateOut &out , value_type dt )
     {
         m_stepper.do_step( system , in , dxdt , t , out , dt );
         ++do_count;
     }
 
     template< class System , class StateInOut , class DerivIn >
-    void do_step( System system , StateInOut &x , const DerivIn &dxdt , value_type t , value_type dt )
+    void do_step_dxdt_impl( System system , StateInOut &x , const DerivIn &dxdt , value_type t , value_type dt )
     {
         m_stepper.do_step( system , x , dxdt , t , dt );
         ++do_count;

test/default_operations.cpp

@@ -229,12 +229,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( scale_sum2_units_test , T , test_types )
     typedef unit_fixture< T > fix_type;
     typedef typename fix_type::value_type value_type;
     typedef typename fix_type::time_type time_type;
-    typedef typename fix_type::time_2_type time_2_type;
-    typedef typename fix_type::time_3_type time_3_type;
-    typedef typename fix_type::time_4_type time_4_type;
-    typedef typename fix_type::time_5_type time_5_type;
-    typedef typename fix_type::time_6_type time_6_type;
-    typedef typename fix_type::time_7_type time_7_type;
 
     fix_type f;
     typedef default_operations::scale_sum2< value_type , time_type > Op;
@@ -249,11 +243,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( scale_sum3_units_test , T , test_types )
     typedef typename fix_type::value_type value_type;
     typedef typename fix_type::time_type time_type;
     typedef typename fix_type::time_2_type time_2_type;
-    typedef typename fix_type::time_3_type time_3_type;
-    typedef typename fix_type::time_4_type time_4_type;
-    typedef typename fix_type::time_5_type time_5_type;
-    typedef typename fix_type::time_6_type time_6_type;
-    typedef typename fix_type::time_7_type time_7_type;
 
     fix_type f;
     typedef default_operations::scale_sum3< value_type , time_type , time_2_type > Op;