diff --git a/test_external/eigen/Jamfile.v2 b/test_external/eigen/Jamfile.v2
new file mode 100644
index 0000000..489f67c
--- /dev/null
+++ b/test_external/eigen/Jamfile.v2
@@ -0,0 +1,36 @@
+# Copyright 2012-2013 Karsten Ahnert
+# Copyright 2012-2013 Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# bring in rules for testing
+
+
+import testing ;
+import os ;
+
+use-project boost : $(BOOST_ROOT) ;
+
+local EIGEN_ROOT = [ os.environ EIGEN_ROOT ] ;
+
+project
+    : requirements
+      <library>/boost/test//boost_unit_test_framework
+      <define>BOOST_ALL_NO_LIB=1
+      <include>$(EIGEN_ROOT)
+      <include>../../test
+      <link>static
+      # <cxxflags>-D_SCL_SECURE_NO_WARNINGS
+    ;
+
+test-suite "odeint"
+    :
+    [ compile is_resizeable.cpp ]
+    [ run same_size.cpp ]
+    [ run resize.cpp ]
+    [ run runge_kutta4.cpp ]
+    [ run runge_kutta_dopri5.cpp ]
+    [ run integrate.cpp ]
+    [ compile-fail fail_integrate.cpp ]
+    : <testing.launcher>valgrind
+    ;
diff --git a/test_external/eigen/fail_integrate.cpp b/test_external/eigen/fail_integrate.cpp
new file mode 100644
index 0000000..4c71b01
--- /dev/null
+++ b/test_external/eigen/fail_integrate.cpp
@@ -0,0 +1,47 @@
+/*
+ [auto_generated]
+ fail_integrate.cpp
+
+ [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)
+ */
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_dummy
+
+#include <boost/numeric/odeint/integrate/integrate.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include "dummy_odes.hpp"
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+
+BOOST_AUTO_TEST_SUITE( eigen_fail_integrate )
+
+BOOST_AUTO_TEST_CASE( test )
+{
+    typedef Eigen::Matrix< double , 1 , 1 > state_type;
+    state_type x;
+    x[0] = 10.0;
+    double t_start = 0.0 , t_end = 1.0 , dt = 0.1;
+    integrate( constant_system_functor_standard() , x , t_start , t_end , dt );
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
diff --git a/test_external/eigen/integrate.cpp b/test_external/eigen/integrate.cpp
new file mode 100644
index 0000000..c64496a
--- /dev/null
+++ b/test_external/eigen/integrate.cpp
@@ -0,0 +1,69 @@
+/*
+ [auto_generated]
+ integrate.cpp
+
+ [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)
+ */
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_eigen_integrate
+
+#include <boost/numeric/odeint/integrate/integrate.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+#include "dummy_odes.hpp"
+
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+
+BOOST_AUTO_TEST_SUITE( eigen_integrate )
+
+BOOST_AUTO_TEST_CASE( test_const_sys )
+{
+    typedef Eigen::Matrix< double , 1 , 1 > state_type;
+    state_type x;
+    x[0] = 10.0;
+    double t_start = 0.0 , t_end = 1.0 , dt = 0.1;
+    integrate< double >( constant_system_functor_standard() , x , t_start , t_end , dt );
+    BOOST_CHECK_CLOSE( x[0] , 11.0 , 1.0e-13 );
+}
+
+BOOST_AUTO_TEST_CASE( test_lorenz )
+{
+    typedef Eigen::Matrix< double , 3 , 1 > state_type;
+    state_type x;
+    x[0] = 10.0;
+    x[1] = 10.0;
+    x[2] = 10.0;
+    double t_start = 0.0 , t_end = 1000.0 , dt = 0.1;
+    integrate< double >( lorenz() , x , t_start , t_end , dt );
+    
+    std::vector< double > x2( 3 );
+    x2[0] = 10.0;
+    x2[1] = 10.0;
+    x2[2] = 10.0;
+    integrate( lorenz() , x2 , t_start , t_end , dt );
+    
+    BOOST_CHECK_CLOSE( x[0] , x2[0] , 1.0e-13 );
+    BOOST_CHECK_CLOSE( x[1] , x2[1] , 1.0e-13 );
+    BOOST_CHECK_CLOSE( x[2] , x2[2] , 1.0e-13 );
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/eigen/is_resizeable.cpp b/test_external/eigen/is_resizeable.cpp
new file mode 100644
index 0000000..5f29ee5
--- /dev/null
+++ b/test_external/eigen/is_resizeable.cpp
@@ -0,0 +1,47 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/test_external/eigen/is_resizeable.cpp
+
+ [begin_description]
+ tba.
+ [end_description]
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_eigen_is_resizeable
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/external/eigen/eigen_resize.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+
+BOOST_AUTO_TEST_SUITE( is_resizeable )
+
+BOOST_AUTO_TEST_CASE( test_compile_time_matrix )
+{
+    typedef Eigen::Matrix< double , 1 , 1 > matrix_type;
+    BOOST_STATIC_ASSERT(( boost::numeric::odeint::is_resizeable< matrix_type >::value ));
+}
+
+BOOST_AUTO_TEST_CASE( test_compile_time_array )
+{
+    typedef Eigen::Array< double , 1 , 1 > array_type;
+    BOOST_STATIC_ASSERT(( boost::numeric::odeint::is_resizeable< array_type >::value ));
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/eigen/resize.cpp b/test_external/eigen/resize.cpp
new file mode 100644
index 0000000..4efcdd4
--- /dev/null
+++ b/test_external/eigen/resize.cpp
@@ -0,0 +1,145 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/test_external/eigen/resize.cpp
+
+ [begin_description]
+ tba.
+ [end_description]
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_eigen_resize
+
+#include <boost/test/unit_test.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen_resize.hpp>
+
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+
+BOOST_AUTO_TEST_SUITE( eigen_resize )
+
+BOOST_AUTO_TEST_CASE( test_compile_time_matrix )
+{
+    typedef Eigen::Matrix< double , 1 , 1 > matrix_type;
+    matrix_type a , b;
+    boost::numeric::odeint::resize( a , b );
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+    BOOST_CHECK_EQUAL( a.rows() , 1 );
+    BOOST_CHECK_EQUAL( a.cols() , 1 );
+}
+
+BOOST_AUTO_TEST_CASE( test_rumtime_matrix )
+{
+    typedef Eigen::Matrix< double , Eigen::Dynamic , Eigen::Dynamic > matrix_type;
+    matrix_type a( 5 , 2 ) , b;
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 0 );
+    BOOST_CHECK_EQUAL( b.cols() , 0 );
+    BOOST_CHECK( !boost::numeric::odeint::same_size( a , b ) );
+
+    boost::numeric::odeint::resize( b , a );
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 5 );
+    BOOST_CHECK_EQUAL( b.cols() , 2 );
+
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+BOOST_AUTO_TEST_CASE( test_rumtime_matrix2 )
+{
+    typedef Eigen::Matrix< double , Eigen::Dynamic , Eigen::Dynamic > matrix_type;
+    matrix_type a( 5 , 2 ) , b( 2 , 3 );
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 2 );
+    BOOST_CHECK_EQUAL( b.cols() , 3 );
+    BOOST_CHECK( !boost::numeric::odeint::same_size( a , b ) );
+
+    boost::numeric::odeint::resize( b , a );
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 5 );
+    BOOST_CHECK_EQUAL( b.cols() , 2 );
+
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE( test_compile_time_array )
+{
+    typedef Eigen::Array< double , 1 , 1 > array_type;
+    array_type a , b;
+    boost::numeric::odeint::resize( a , b );
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+    BOOST_CHECK_EQUAL( a.rows() , 1 );
+    BOOST_CHECK_EQUAL( a.cols() , 1 );
+}
+
+BOOST_AUTO_TEST_CASE( test_rumtime_array )
+{
+    typedef Eigen::Array< double , Eigen::Dynamic , Eigen::Dynamic > array_type;
+    array_type a( 5 , 2 ) , b;
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 0 );
+    BOOST_CHECK_EQUAL( b.cols() , 0 );
+    BOOST_CHECK( !boost::numeric::odeint::same_size( a , b ) );
+
+    boost::numeric::odeint::resize( b , a );
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 5 );
+    BOOST_CHECK_EQUAL( b.cols() , 2 );
+
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+BOOST_AUTO_TEST_CASE( test_rumtime_array2 )
+{
+    typedef Eigen::Array< double , Eigen::Dynamic , Eigen::Dynamic > array_type;
+    array_type a( 5 , 2 ) , b( 2 , 3 );
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 2 );
+    BOOST_CHECK_EQUAL( b.cols() , 3 );
+    BOOST_CHECK( !boost::numeric::odeint::same_size( a , b ) );
+
+    boost::numeric::odeint::resize( b , a );
+
+    BOOST_CHECK_EQUAL( a.rows() , 5 );
+    BOOST_CHECK_EQUAL( a.cols() , 2 );
+    BOOST_CHECK_EQUAL( b.rows() , 5 );
+    BOOST_CHECK_EQUAL( b.cols() , 2 );
+
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/eigen/runge_kutta4.cpp b/test_external/eigen/runge_kutta4.cpp
new file mode 100644
index 0000000..c4054fd
--- /dev/null
+++ b/test_external/eigen/runge_kutta4.cpp
@@ -0,0 +1,94 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/test_external/eigen/runge_kutta4.cpp
+
+ [begin_description]
+ tba.
+ [end_description]
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_eigen_runge_kutta4
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen_resize.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+struct sys
+{
+    template< class State , class Deriv >
+    void operator()( const State &x , Deriv &dxdt , double t ) const
+    {
+        dxdt[0] = 1.0;
+    }
+};
+
+
+BOOST_AUTO_TEST_SUITE( eigen_runge_kutta4 )
+
+BOOST_AUTO_TEST_CASE( compile_time_matrix )
+{
+    typedef Eigen::Matrix< double , 1 , 1 > state_type;
+    state_type x;
+    x[0] = 10.0;
+    runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    rk4.do_step( sys() , x , 0.0 , 0.1 );
+    BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+    
+}
+
+BOOST_AUTO_TEST_CASE( runtime_matrix )
+{
+    typedef Eigen::Matrix< double , Eigen::Dynamic , 1 > state_type;
+    state_type x( 1 );
+    x[0] = 10.0;
+    runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    rk4.do_step( sys() , x , 0.0 , 0.1 );
+    BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+}
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE( compile_time_array )
+{
+    typedef Eigen::Array< double , 1 , 1 > state_type;
+    state_type x;
+    x[0] = 10.0;
+    runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    rk4.do_step( sys() , x , 0.0 , 0.1 );
+    BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+    
+}
+
+BOOST_AUTO_TEST_CASE( runtime_array )
+{
+    typedef Eigen::Array< double , Eigen::Dynamic , 1 > state_type;
+    state_type x( 1 );
+    x[0] = 10.0;
+    runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    rk4.do_step( sys() , x , 0.0 , 0.1 );
+    BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+}
+
+
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/eigen/runge_kutta_dopri5.cpp b/test_external/eigen/runge_kutta_dopri5.cpp
new file mode 100644
index 0000000..516f04f
--- /dev/null
+++ b/test_external/eigen/runge_kutta_dopri5.cpp
@@ -0,0 +1,132 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/test_external/eigen/runge_kutta_dopri5.cpp
+
+ [begin_description]
+ tba.
+ [end_description]
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_eigen_runge_kutta4
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
+#include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
+
+// #include <boost/numeric/odeint/external/eigen/eigen_resize.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>
+
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+struct sys
+{
+    template< class State , class Deriv >
+    void operator()( const State &x , Deriv &dxdt , double t ) const
+    {
+        dxdt[0] = 1.0;
+    }
+};
+
+template< class State >
+struct stepper
+{
+    typedef runge_kutta_dopri5< State , double , State , double , vector_space_algebra > type;
+};
+
+template< class State >
+struct controlled_stepper
+{
+    typedef controlled_runge_kutta< typename stepper< State >::type > type;
+};
+
+template< class State >
+struct dense_output_stepper
+{
+    typedef dense_output_runge_kutta< typename controlled_stepper< State >::type > type;
+};
+
+
+
+
+BOOST_AUTO_TEST_SUITE( eigen_runge_kutta_dopri5 )
+
+BOOST_AUTO_TEST_CASE( compile_time_matrix )
+{
+    typedef Eigen::Matrix< double , 1 , 1 > state_type;
+    state_type x;
+    x[0] = 10.0;
+    double t = 0.0 , dt = 0.1 ;
+    
+    // dense_output_stepper< state_type >::type s;
+    // s.initialize( x , t , dt );
+
+    // controlled_stepper< state_type >::type s;
+    // s.try_step( sys() , x , t , dt );
+
+    stepper< state_type >::type s;
+    s.do_step( sys() , x , t , dt );
+
+    // runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    // rk4.do_step( sys() , x , 0.0 , 0.1 );
+    // BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+    
+}
+
+BOOST_AUTO_TEST_CASE( runtime_matrix )
+{
+    typedef Eigen::Matrix< double , Eigen::Dynamic , 1 > state_type;
+    state_type x( 1 );
+    x[0] = 10.0;
+    
+    // runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    // rk4.do_step( sys() , x , 0.0 , 0.1 );
+    // BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+}
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE( compile_time_array )
+{
+    typedef Eigen::Array< double , 1 , 1 > state_type;
+    state_type x;
+    x[0] = 10.0;
+    // runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    // rk4.do_step( sys() , x , 0.0 , 0.1 );
+    // BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+    
+}
+
+BOOST_AUTO_TEST_CASE( runtime_array )
+{
+    typedef Eigen::Array< double , Eigen::Dynamic , 1 > state_type;
+    state_type x( 1 );
+    x[0] = 10.0;
+    // runge_kutta4< state_type , double , state_type , double , vector_space_algebra > rk4;
+    // rk4.do_step( sys() , x , 0.0 , 0.1 );
+    // BOOST_CHECK_CLOSE( x[0] , 10.1 , 1.0e-13 );
+}
+
+
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/eigen/same_size.cpp b/test_external/eigen/same_size.cpp
new file mode 100644
index 0000000..ea771f9
--- /dev/null
+++ b/test_external/eigen/same_size.cpp
@@ -0,0 +1,83 @@
+/*
+ [auto_generated]
+ libs/numeric/odeint/test_external/eigen/same_size.cpp
+
+ [begin_description]
+ tba.
+ [end_description]
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_eigen_same_size
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/external/eigen/eigen_resize.hpp>
+
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+
+BOOST_AUTO_TEST_SUITE( eigen_same_size )
+
+BOOST_AUTO_TEST_CASE( compile_time_matrix )
+{
+    typedef Eigen::Matrix< double , 1 , 1 > matrix_type;
+    matrix_type a , b;
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+BOOST_AUTO_TEST_CASE( runtime_matrix )
+{
+    typedef Eigen::Matrix< double , Eigen::Dynamic , Eigen::Dynamic > matrix_type;
+    matrix_type a( 10 , 2 ) , b( 10 , 2 );
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+BOOST_AUTO_TEST_CASE( fail_runtime_matrix )
+{
+    typedef Eigen::Matrix< double , Eigen::Dynamic , Eigen::Dynamic > matrix_type;
+    matrix_type a( 11 , 2 ) , b( 10 , 2 );
+    BOOST_CHECK( !boost::numeric::odeint::same_size( a , b ) );
+
+}
+
+
+
+BOOST_AUTO_TEST_CASE( compile_time_array )
+{
+    typedef Eigen::Array< double , 1 , 1 > array_type;
+    array_type a , b;
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+BOOST_AUTO_TEST_CASE( runtime_array )
+{
+    typedef Eigen::Array< double , Eigen::Dynamic , Eigen::Dynamic > array_type;
+    array_type a( 10 , 2 ) , b( 10 , 2 );
+    BOOST_CHECK( boost::numeric::odeint::same_size( a , b ) );
+}
+
+BOOST_AUTO_TEST_CASE( fail_runtime_array )
+{
+    typedef Eigen::Array< double , Eigen::Dynamic , Eigen::Dynamic > array_type;
+    array_type a( 11 , 2 ) , b( 10 , 2 );
+    BOOST_CHECK( !boost::numeric::odeint::same_size( a , b ) );
+
+}
+
+
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/gmp/Jamfile.v2 b/test_external/gmp/Jamfile.v2
new file mode 100644
index 0000000..f0e6af0
--- /dev/null
+++ b/test_external/gmp/Jamfile.v2
@@ -0,0 +1,26 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2012 Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# bring in rules for testing
+
+import testing ;
+use-project boost : $(BOOST_ROOT) ;
+
+project gmp
+    : requirements
+      <library>/boost/test//boost_unit_test_framework
+    ;
+    
+    
+lib libgmp : : <name>gmp <link>shared ;
+lib libgmpxx : : <name>gmpxx <link>shared ;
+
+test-suite "gmp"
+    : 
+      [ run check_gmp.cpp libgmpxx libgmp : : : <link>shared:<define>BOOST_TEST_DYN_LINK=1 ] 
+      [ run gmp_integrate.cpp libgmpxx libgmp : : : <link>shared:<define>BOOST_TEST_DYN_LINK=1 ] 
+    ;
+
+
diff --git a/test_external/gmp/check_gmp.cpp b/test_external/gmp/check_gmp.cpp
new file mode 100644
index 0000000..04b3e10
--- /dev/null
+++ b/test_external/gmp/check_gmp.cpp
@@ -0,0 +1,165 @@
+/* Boost check_gmp.cpp test file
+
+ Copyright 2010-2012 Mario Mulansky
+ Copyright 2011-2012 Karsten Ahnert
+
+ This file tests the odeint library with the gmp arbitrary precision types
+
+ 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)
+*/
+
+#define BOOST_TEST_MODULE odeint_gmp
+
+#include <iostream>
+
+#include <gmpxx.h>
+
+#include <boost/test/unit_test.hpp>
+#include <boost/array.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+//#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+namespace mpl = boost::mpl;
+
+const int precision = 1024;
+
+typedef mpf_class value_type;
+typedef mpf_class state_type;
+
+//provide min, max and pow functions for mpf types - required for controlled steppers
+value_type min( const value_type a , const value_type b )
+{
+    if( a<b ) return a;
+    else return b;
+}
+value_type max( const value_type a , const value_type b )
+{
+    if( a>b ) return a;
+    else return b;
+}
+value_type pow( const value_type a , const value_type b )
+{
+    // do calculation in double precision
+    return value_type( std::pow( a.get_d() , b.get_d() ) );
+}
+
+
+//provide vector_space reduce:
+
+namespace boost { namespace numeric { namespace odeint {
+
+template<>
+struct vector_space_reduce< state_type >
+{
+  template< class Op >
+  state_type operator()( state_type x , Op op , state_type init ) const
+  {
+      init = op( init , x );
+      return init;
+  }
+};
+
+} } }
+
+
+void constant_system( const state_type &x , state_type &dxdt , value_type t )
+{
+    dxdt = value_type( 1.0 , precision );
+}
+
+
+/* check runge kutta stepers */
+typedef mpl::vector<
+    euler< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    modified_midpoint< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta4< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta4_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_cash_karp54< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra >
+    > stepper_types;
+
+
+template< class Stepper >
+struct perform_runge_kutta_test {
+
+    void operator()( void )
+    {
+        /* We have to specify the desired precision in advance! */
+        mpf_set_default_prec( precision );
+
+        mpf_t eps_ , unity;
+        mpf_init( eps_ ); mpf_init( unity );
+        mpf_set_d( unity , 1.0 );
+        mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
+        value_type eps( eps_ );
+
+        Stepper stepper;
+        state_type x;
+        x = 0.0;
+
+        stepper.do_step( constant_system , x , 0.0 , 0.1 );
+
+        BOOST_MESSAGE( eps );
+        BOOST_CHECK_MESSAGE( abs( x - value_type( 0.1 , precision ) ) < eps , x - 0.1 );
+    }
+};
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_stepper_test , Stepper , stepper_types )
+{
+    perform_runge_kutta_test< Stepper > tester;
+    tester();
+}
+
+
+/* check controlled steppers */
+typedef mpl::vector<
+    controlled_runge_kutta< runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
+    controlled_runge_kutta< runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > > , 
+    controlled_runge_kutta< runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
+    bulirsch_stoer< state_type , value_type , state_type , value_type , vector_space_algebra >
+    > controlled_stepper_types;
+
+
+template< class Stepper >
+struct perform_controlled_test {
+
+    void operator()( void )
+    {
+        mpf_set_default_prec( precision );
+
+        mpf_t eps_ , unity;
+        mpf_init( eps_ ); mpf_init( unity );
+        mpf_set_d( unity , 1.0 );
+        mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
+        value_type eps( eps_ );
+
+        Stepper stepper;
+        state_type x;
+        x = 0.0;
+
+        value_type t(0.0);
+        value_type dt(0.1);
+
+        stepper.try_step( constant_system , x , t , dt );
+
+        BOOST_MESSAGE( eps );
+        BOOST_CHECK_MESSAGE( abs( x - value_type( 0.1 , precision ) ) < eps , x - 0.1 );
+    }
+};
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( controlled_stepper_test , Stepper , controlled_stepper_types )
+{
+    perform_controlled_test< Stepper > tester;
+    tester();
+}
diff --git a/test_external/gmp/gmp_integrate.cpp b/test_external/gmp/gmp_integrate.cpp
new file mode 100644
index 0000000..0a6be30
--- /dev/null
+++ b/test_external/gmp/gmp_integrate.cpp
@@ -0,0 +1,169 @@
+/* Boost check_gmp.cpp test file
+
+ Copyright 2010-2012 Mario Mulansky
+ Copyright 2011-2012 Karsten Ahnert
+
+ This file tests the odeint library with the gmp arbitrary precision types
+
+ 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)
+*/
+
+#define BOOST_TEST_MODULE odeint_gmp
+
+#include <gmpxx.h>
+
+#include <boost/test/unit_test.hpp>
+#include <boost/array.hpp>
+
+#include <boost/mpl/vector.hpp>
+
+#include <boost/numeric/odeint.hpp>
+//#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+namespace mpl = boost::mpl;
+
+const int precision = 1024;
+
+typedef mpf_class value_type;
+typedef mpf_class state_type;
+
+//provide min, max and pow functions for mpf types - required for controlled steppers
+value_type min( const value_type a , const value_type b )
+{
+    if( a<b ) return a;
+    else return b;
+}
+value_type max( const value_type a , const value_type b )
+{
+    if( a>b ) return a;
+    else return b;
+}
+value_type pow( const value_type a , const value_type b )
+{
+    // do the calculation in double precision...
+    return value_type( std::pow( a.get_d() , b.get_d() ) );
+}
+
+
+//provide vector_space reduce:
+
+namespace boost { namespace numeric { namespace odeint {
+
+template<>
+struct vector_space_reduce< state_type >
+{
+  template< class Op >
+  state_type operator()( state_type x , Op op , state_type init ) const
+  {
+      init = op( init , x );
+      return init;
+  }
+};
+
+} } }
+
+
+void constant_system( const state_type &x , state_type &dxdt , value_type t )
+{
+    dxdt = value_type( 1.0 , precision );
+}
+
+/* check runge kutta stepers */
+typedef mpl::vector<
+    euler< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    modified_midpoint< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta4< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta4_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_cash_karp54< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > ,
+    runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra >
+    > stepper_types;
+
+
+template< class Stepper >
+struct perform_integrate_const_test {
+
+    void operator()( void )
+    {
+        /* We have to specify the desired precision in advance! */
+        mpf_set_default_prec( precision );
+
+        mpf_t eps_ , unity;
+        mpf_init( eps_ ); mpf_init( unity );
+        mpf_set_d( unity , 1.0 );
+        mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
+        value_type eps( eps_ );
+
+        Stepper stepper;
+        state_type x;
+        x = 0.0;
+        value_type t0( 0.0 );
+        value_type tend( 1.0 );
+        value_type dt(0.1);
+
+        integrate_const( stepper , constant_system , x , t0 , tend , dt );
+
+        x = 0.0;
+        t0 = 0.0;
+        dt = 0.1;
+        size_t steps = 10;
+        
+        integrate_n_steps( stepper , constant_system , x , t0 , dt , steps );
+
+        BOOST_CHECK_MESSAGE( abs( x - 10*dt ) < eps , x );
+    }
+};
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_const_test , Stepper , stepper_types )
+{
+    perform_integrate_const_test< Stepper > tester;
+    tester();
+}
+
+
+typedef mpl::vector<
+    controlled_runge_kutta< runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
+    controlled_runge_kutta< runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > > , 
+    controlled_runge_kutta< runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
+    bulirsch_stoer< state_type , value_type , state_type , value_type , vector_space_algebra >
+    > controlled_stepper_types;
+
+
+template< class Stepper >
+struct perform_integrate_adaptive_test {
+
+    void operator()( void )
+    {
+        mpf_set_default_prec( precision );
+
+        mpf_t eps_ , unity;
+        mpf_init( eps_ ); mpf_init( unity );
+        mpf_set_d( unity , 1.0 );
+        mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
+        value_type eps( eps_ );
+
+        Stepper stepper;
+        state_type x;
+        x = 0.0;
+        value_type t0( 0.0 );
+        value_type tend( 1.0 );
+        value_type dt(0.1);
+
+        integrate_adaptive( stepper , constant_system , x , t0 , tend , dt );
+
+        BOOST_CHECK_MESSAGE( abs( x - tend ) < eps , x - 0.1 );
+    }
+};
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_adaptive__test , Stepper , controlled_stepper_types )
+{
+    perform_integrate_adaptive_test< Stepper > tester;
+    tester();
+}
diff --git a/test_external/gsl/Jamfile.v2 b/test_external/gsl/Jamfile.v2
new file mode 100644
index 0000000..086edf6
--- /dev/null
+++ b/test_external/gsl/Jamfile.v2
@@ -0,0 +1,28 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2012 Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# bring in rules for testing
+
+import testing ;
+use-project boost : $(BOOST_ROOT) ;
+
+project
+    : requirements
+      <library>/boost/test//boost_unit_test_framework
+    ;
+    
+    
+lib libgsl : : <name>gsl <link>shared ;
+lib libgslcblas : : <name>gslcblas <link>shared ;
+
+test-suite "gsl"
+    : 
+      [ run check_gsl.cpp libgslcblas libgsl
+      :
+      :
+      : <link>shared:<define>BOOST_TEST_DYN_LINK=1
+      ] 
+    ;
+
diff --git a/test_external/gsl/check_gsl.cpp b/test_external/gsl/check_gsl.cpp
new file mode 100644
index 0000000..3e7856d
--- /dev/null
+++ b/test_external/gsl/check_gsl.cpp
@@ -0,0 +1,58 @@
+/* Boost check_gmp.cpp test file
+
+ Copyright 2010-2011 Karsten Ahnert
+ Copyright 2011 Mario Mulansky
+
+ This file tests the odeint library with the gmp arbitrary precision types
+
+ 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)
+*/
+
+#define BOOST_TEST_MODULE odeint_gsl
+
+#include <gsl/gsl_vector.h>
+#include <boost/numeric/odeint/stepper/euler.hpp>
+#include <boost/numeric/odeint/external/gsl/gsl_wrapper.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef gsl_vector *state_type;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type x , state_type dxdt , double t )
+{
+    gsl_vector_set( dxdt , 0 , sigma * ( gsl_vector_get(x , 1 ) - gsl_vector_get( x , 0 ) ) );
+    gsl_vector_set( dxdt , 1 , R * gsl_vector_get( x , 0 ) - gsl_vector_get( x , 1 ) - gsl_vector_get( x , 0 ) * gsl_vector_get( x , 2) );
+    gsl_vector_set( dxdt , 2 , gsl_vector_get( x , 0 ) * gsl_vector_get( x , 1 ) - b * gsl_vector_get( x , 2) );
+}
+
+BOOST_AUTO_TEST_CASE( gsl )
+{
+    euler< state_type > euler;
+
+    state_type x = gsl_vector_alloc( 3 );
+
+    // check resizing
+    state_type y = 0;
+    boost::numeric::odeint::resize( y , x );
+    BOOST_CHECK( 0 != y );
+
+    gsl_vector_set( x , 0 , 1.0);
+    gsl_vector_set( x , 1 , 1.0);
+    gsl_vector_set( x , 2 , 2.0);
+
+    euler.do_step( lorenz , x , 0.0 , 0.1 );
+
+    //cout << gsl_vector_get( x , 0 ) << "  " << gsl_vector_get( x , 1 ) << "  " << gsl_vector_get( x , 2 ) << endl;
+
+    gsl_vector_free( x );
+
+}
diff --git a/test_external/mkl/Jamfile.v2 b/test_external/mkl/Jamfile.v2
new file mode 100644
index 0000000..f585fa8
--- /dev/null
+++ b/test_external/mkl/Jamfile.v2
@@ -0,0 +1,30 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2012 Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# bring in rules for testing
+
+import testing ;
+use-project boost : $(BOOST_ROOT) ;
+
+project
+    : requirements
+      <library>/boost/test//boost_unit_test_framework
+    ;
+    
+    
+lib libmkl : : <name>mkl_intel_lp64 <link>shared ;
+lib libmkl_core : : <name>mkl_core <link>shared ;
+lib libmkl_intel_thread : : <name>mkl_intel_thread ;
+lib libiomp5 : : <name>iomp5 ;
+lib libpthread : : <name>pthread ;
+
+test-suite "mkl"
+    : 
+      [ run check_mkl.cpp  libpthread libiomp5 libmkl_core libmkl_intel_thread libmkl
+      :
+      :
+      : <link>shared:<define>BOOST_TEST_DYN_LINK=1
+      ] 
+    ;
diff --git a/test_external/mkl/check_mkl.cpp b/test_external/mkl/check_mkl.cpp
new file mode 100644
index 0000000..bedbd66
--- /dev/null
+++ b/test_external/mkl/check_mkl.cpp
@@ -0,0 +1,51 @@
+/* Boost check_mkl.cpp test file
+
+ Copyright 2010-2011 Mario Mulansky
+ Copyright 2011 Karsten Ahnert
+
+ This file tests the odeint library with the intel mkl blas1 routines
+
+ 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)
+*/
+
+#define BOOST_TEST_MODULE test_mkl
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/stepper/euler.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>
+
+using namespace boost::numeric::odeint;
+
+typedef double value_type;
+typedef boost::array< value_type , 1 > state_type;
+
+
+void constant_system( state_type &x , state_type &dxdt , value_type t )
+{
+    dxdt[0] = 1.0;
+}
+
+const double eps = 1E-14;
+
+
+BOOST_AUTO_TEST_CASE( test_mkl )
+{
+
+    //to use mkl routines we have to use the vector_space_algebra and the mkl_operations
+    runge_kutta4< state_type , value_type , state_type , value_type , vector_space_algebra , mkl_operations > stepper;
+    state_type x;
+    x[0] = 0.0;
+
+    stepper.do_step( constant_system , x , 0.0 , 0.1 );
+
+    using std::abs;
+
+    std::cout << x[0] << " ?= " << 0.1 << std::endl;
+    BOOST_CHECK_SMALL( abs( x[0] - 0.1 ) , eps );
+
+}
diff --git a/test_external/mpi/Jamfile.v2 b/test_external/mpi/Jamfile.v2
new file mode 100644
index 0000000..0a02bcc
--- /dev/null
+++ b/test_external/mpi/Jamfile.v2
@@ -0,0 +1,27 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2012-2013 Mario Mulansky
+# Copyright 2013 Pascal Germroth
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+import testing ;
+import mpi : mpi-test ;
+
+use-project boost : $(BOOST_ROOT) ;
+
+project
+    : requirements
+      <library>/boost/test//boost_unit_test_framework
+      <library>/boost//mpi
+      <link>static
+      <define>BOOST_ALL_NO_LIB=1
+    ;
+
+# mpi-test name : source : req : np=1 2 3 4 7 8 13 17
+test-suite "odeint-mpi"
+    :
+    [ mpi-test split_test ]
+    [ mpi-test state_test ]
+    [ mpi-test norm_test ]
+    ;
+
diff --git a/test_external/mpi/norm_test.cpp b/test_external/mpi/norm_test.cpp
new file mode 100644
index 0000000..69f1d40
--- /dev/null
+++ b/test_external/mpi/norm_test.cpp
@@ -0,0 +1,62 @@
+/*
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <iostream>
+#include <sstream>
+#include <cstdlib>
+
+#define BOOST_TEST_MODULE odeint_mpi
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/external/mpi/mpi.hpp>
+
+using namespace boost::numeric::odeint;
+
+boost::mpi::environment env;
+
+BOOST_AUTO_TEST_SUITE( norm_test_suite )
+
+BOOST_AUTO_TEST_CASE( norm_test )
+{
+    boost::mpi::communicator world;
+
+    int ref_value = 0;
+    std::vector<int> in_data;
+    mpi_state< std::vector<int> > state(world);
+
+    // generate data and reference value on master
+    if(world.rank() == 0) {
+        for(size_t i = 0 ; i < 400 ; i++)
+            in_data.push_back( rand() % 10000 );
+        ref_value = *std::max_element(in_data.begin(), in_data.end());
+    }
+    boost::mpi::broadcast(world, ref_value, 0);
+
+    // copy to nodes
+    split( in_data, state );
+
+    int value = mpi_nested_algebra< range_algebra >::norm_inf( state );
+
+    {
+        std::ostringstream ss;
+        ss << "state[" << world.rank() << "]"
+           << " local:" << range_algebra::norm_inf( state() )
+           << " global:" << value
+           << " ref:" << ref_value << "\n";
+        std::clog << ss.str() << std::flush;
+    }
+
+    BOOST_REQUIRE_EQUAL( value, ref_value );
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
diff --git a/test_external/mpi/split_test.cpp b/test_external/mpi/split_test.cpp
new file mode 100644
index 0000000..f78a601
--- /dev/null
+++ b/test_external/mpi/split_test.cpp
@@ -0,0 +1,61 @@
+/*
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <iostream>
+#include <sstream>
+
+#define BOOST_TEST_MODULE odeint_mpi
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/external/mpi/mpi.hpp>
+
+using namespace boost::numeric::odeint;
+
+boost::mpi::environment env;
+
+BOOST_AUTO_TEST_SUITE( split_test_suite )
+
+BOOST_AUTO_TEST_CASE( split_test )
+{
+    boost::mpi::communicator world;
+
+    const size_t total_size = 31;
+
+    std::vector<size_t> in_data, out_data;
+    mpi_state< std::vector<size_t> > state(world);
+
+    // generate data on master
+    if(world.rank() == 0)
+        for(size_t i = 0 ; i < total_size ; i++) in_data.push_back(i);
+
+    // copy to nodes
+    split( in_data, state );
+
+    BOOST_REQUIRE((state().size() == total_size / world.size())
+               || (state().size() == total_size / world.size() + 1));
+
+    {
+        std::ostringstream ss;
+        ss << "state[" << world.rank() << "].data = {";
+        std::copy(state().begin(), state().end(), std::ostream_iterator<size_t>(ss, ", "));
+        ss << "}\n";
+        std::clog << ss.str() << std::flush;
+    }
+
+    // copy back to master
+    if(world.rank() == 0) out_data.resize(in_data.size());
+    unsplit( state, out_data );
+
+    if(world.rank() == 0)
+        BOOST_REQUIRE_EQUAL_COLLECTIONS(in_data.begin(), in_data.end(), out_data.begin(), out_data.end());
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/mpi/state_test.cpp b/test_external/mpi/state_test.cpp
new file mode 100644
index 0000000..490b46e
--- /dev/null
+++ b/test_external/mpi/state_test.cpp
@@ -0,0 +1,78 @@
+/*
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include <iostream>
+#include <sstream>
+
+#define BOOST_TEST_MODULE odeint_mpi
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/external/mpi/mpi.hpp>
+
+using namespace boost::numeric::odeint;
+
+boost::mpi::environment env;
+
+BOOST_AUTO_TEST_SUITE( state_test_suite )
+
+BOOST_AUTO_TEST_CASE( state_test )
+{
+    boost::mpi::communicator world;
+
+    std::vector<size_t> in_data1, in_data2;
+    mpi_state< std::vector<size_t> > state1(world), state2(world);
+
+    // generate data on master
+    if(world.rank() == 0) {
+        in_data1.resize(31);
+        in_data2.resize(33);
+        for(size_t i = 0 ; i < in_data2.size() ; i++)
+            in_data2[i] = i;
+    }
+
+    // copy to nodes
+    split( in_data1, state1 );
+    split( in_data2, state2 );
+
+    {
+        std::ostringstream ss;
+        ss << "state[" << world.rank() << "] {"
+           << state1().size() << ", "
+           << state2().size() << "}\n";
+        std::clog << ss.str() << std::flush;
+    }
+
+    // compare size
+    BOOST_REQUIRE( !same_size( state1, state2 ) );
+
+    // resize state1 to match state2.
+    resize( state1, state2 );
+
+    {
+        std::ostringstream ss;
+        ss << "state[" << world.rank() << "] 1:"
+           << state1().size() << " 2:"
+           << state2().size() << "\n";
+        std::clog << ss.str() << std::flush;
+    }
+
+    // compare size
+    BOOST_REQUIRE( same_size( state1, state2 ) );
+
+    // copy state2 to state1
+    copy( state2, state1 );
+
+    BOOST_REQUIRE_EQUAL_COLLECTIONS(state1().begin(), state1().end(),
+        state2().begin(), state2().end());
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
+
diff --git a/test_external/mtl4/Jamfile.v2 b/test_external/mtl4/Jamfile.v2
new file mode 100644
index 0000000..de9d87f
--- /dev/null
+++ b/test_external/mtl4/Jamfile.v2
@@ -0,0 +1,30 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2013 Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0. (See
+# accompanying file LICENSE_1_0.txt or copy at
+# http://www.boost.org/LICENSE_1_0.txt)
+
+import testing ;
+import boost ;
+
+# boost.use-project ;
+use-project boost : $(BOOST_ROOT) ;
+
+# set your MTL4 directory here
+MTL4_INCLUDE = /home/mario/MTL4/usr/include ;
+
+project
+    : requirements 
+      <library>/boost/test//boost_unit_test_framework
+      <include>$(MTL4_INCLUDE)
+      <define>BOOST_ALL_NO_LIB=1
+      <link>static
+    : 
+    : default-build release
+    ;
+    
+test-suite "odeint-mtl4"
+    :
+    [ run mtl4_resize.cpp ]
+    : <testing.launcher>valgrind
+    ;
diff --git a/test_external/mtl4/mtl4_resize.cpp b/test_external/mtl4/mtl4_resize.cpp
new file mode 100644
index 0000000..7ef985d
--- /dev/null
+++ b/test_external/mtl4/mtl4_resize.cpp
@@ -0,0 +1,89 @@
+/* Boost mtl4_resize.cpp test file
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 Mario Mulansky
+
+ This file tests the odeint library with the mtl4 routines.
+
+ 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)
+*/
+
+#define BOOST_TEST_MODULE test_mtl4_resize
+
+#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/external/mtl4/mtl4_resize.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+
+namespace odeint = boost::numeric::odeint;
+
+
+BOOST_AUTO_TEST_CASE( test_dense_vector_resizeability )
+{
+    BOOST_CHECK( odeint::is_resizeable< mtl::dense_vector< double > >::value );
+}
+
+BOOST_AUTO_TEST_CASE( test_dense2D_resizeability )
+{
+    BOOST_CHECK( odeint::is_resizeable< mtl::dense2D< double > >::value );
+}
+
+BOOST_AUTO_TEST_CASE( test_compressed2D_resizeability )
+{
+    BOOST_CHECK( odeint::is_resizeable< mtl::compressed2D< double > >::value );
+}
+
+
+
+BOOST_AUTO_TEST_CASE( test_dense_vector_vector_same_size )
+{
+    mtl::dense_vector< double > v1( 10 ) , v2( 10 );
+    BOOST_CHECK( odeint::same_size( v2 , v1 ) );
+}
+
+BOOST_AUTO_TEST_CASE( test_dense_vector_dense2D_same_size )
+{
+    mtl::dense_vector< double > v( 10 );
+    mtl::dense2D< double > m( 10 , 10 );
+    BOOST_CHECK( odeint::same_size( m , v ) );
+}
+
+BOOST_AUTO_TEST_CASE( test_dense_vector_compressed2D_same_size )
+{
+    mtl::dense_vector< double > v( 10 );
+    mtl::compressed2D< double > m( 10 , 10 );
+    BOOST_CHECK( odeint::same_size( m , v ) );
+}
+
+
+
+
+BOOST_AUTO_TEST_CASE( test_dense_vector_vector_resize )
+{
+    mtl::dense_vector< double > v1( 10 );
+    mtl::dense_vector< double > v2;
+    odeint::resize( v2 , v1 );
+    BOOST_CHECK( mtl::size( v2 ) == mtl::size( v1 ) );
+}
+
+BOOST_AUTO_TEST_CASE( test_dense_vector_dense2D_resize )
+{
+    mtl::dense_vector< double > v( 10 );
+    mtl::dense2D< double > m;
+
+    odeint::resize( m , v );
+    BOOST_CHECK( m.num_cols() == mtl::size( v ) );
+    BOOST_CHECK( m.num_rows() == mtl::size( v ) );
+}
+
+BOOST_AUTO_TEST_CASE( test_dense_vector_compressed2D_resize )
+{
+    mtl::dense_vector< double > v( 10 );
+    mtl::compressed2D< double > m;
+
+    odeint::resize( m , v );
+    BOOST_CHECK( m.num_cols() == mtl::size( v ) );
+    BOOST_CHECK( m.num_rows() == mtl::size( v ) );
+}
diff --git a/test_external/nt2/Jamfile.v2 b/test_external/nt2/Jamfile.v2
new file mode 100644
index 0000000..26763da
--- /dev/null
+++ b/test_external/nt2/Jamfile.v2
@@ -0,0 +1,44 @@
+#==============================================================================
+#         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+#         Copyright 2014          NumScale SAS
+#
+#          Distributed under the Boost Software License, Version 1.0.
+#                 See accompanying file LICENSE.txt or copy at
+#                     http://www.boost.org/LICENSE_1_0.txt
+#==============================================================================
+
+import testing ;
+import os ;
+
+# This must be built using an NT2 installation.
+# NT2_ROOT_PATH should point to the build directory.
+# Currently, cxxflags needs to be set to the required architecture
+# if using avx/avx2, set the environemnt variable NT2_SIMD_FLAGS to the
+# required value for your compiler (i.e. -mavx2 on g++)
+# If using sse2/3/4 in 64 bits, this is set automatically.
+
+local NT2_ROOT_PATH = [ os.environ NT2_ROOT_PATH ] ;
+local NT2_SIMD_FLAGS = [ os.environ NT2_SIMD_FLAGS ] ;
+
+use-project boost : $(BOOST_ROOT) ;
+
+project
+  : requirements
+    <library>$(BOOST_ROOT)/boost/test/included/unit_test_framework.hpp
+    <define>BOOST_ALL_NO_LIB=1
+    <include>$(NT2_ROOT_PATH)/include/
+    <link>static
+    <toolset>gcc:<cxxflags>-DBOOST_SIMD_NO_STRICT_ALIASING
+    <toolset>gcc:<cxxflags>-fno-strict-aliasing
+    <cxxflags>$(NT2_SIMD_FLAGS)
+  ;
+
+test-suite "odeint"
+  :
+    [ run copy.cpp ]
+    [ run norm_inf.cpp ]
+    [ run resize.cpp ]
+    [ run is_resizeable.cpp ]
+    [ run algebra_dispatcher.cpp ]
+    : <testing.launcher>valgrind
+  ;
diff --git a/test_external/nt2/algebra_dispatcher.cpp b/test_external/nt2/algebra_dispatcher.cpp
new file mode 100644
index 0000000..bea6349
--- /dev/null
+++ b/test_external/nt2/algebra_dispatcher.cpp
@@ -0,0 +1,55 @@
+//==============================================================================
+//         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+//         Copyright 2014          NumScale SAS
+//
+//          Distributed under the Boost Software License, Version 1.0.
+//                 See accompanying file LICENSE.txt or copy at
+//                     http://www.boost.org/LICENSE_1_0.txt
+//==============================================================================
+#define BOOST_TEST_MODULE odeint_nt2_algebra_dispatcher
+
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/numeric/odeint/external/nt2/nt2_algebra_dispatcher.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+#include <boost/mpl/list.hpp>
+
+#include <boost/preprocessor/repetition.hpp>
+#include <boost/preprocessor/arithmetic/mul.hpp>
+
+#include <nt2/table.hpp>
+#include <nt2/sdk/meta/as.hpp>
+#include <nt2/include/functions/ones.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef boost::mpl::list< float , double > fp_types;
+
+#define TABLE(z,n,text) nt2::table<T> y ## n =                                   \
+                        nt2::ones(1,2,nt2::meta::as_<T>() )*T(BOOST_PP_ADD(n,1));
+
+#define PARAMS(z,n,text) T(BOOST_PP_ADD(n,1)),
+
+#define SUM(z,n,text) +BOOST_PP_MUL(BOOST_PP_ADD(n,3),BOOST_PP_ADD(n,2))
+
+#define TEST(z,n,text) BOOST_CHECK_SMALL( y0(BOOST_PP_ADD(n,1))                  \
+                       -T( 2 BOOST_PP_REPEAT(text, SUM, text) ), T(1e-10) );
+
+#define TEST_CASE(z,n,text) BOOST_AUTO_TEST_CASE_TEMPLATE (                      \
+                            BOOST_PP_CAT(odeint_foreach, n), T, fp_types )       \
+{                                                                                \
+  vector_space_algebra algebra;                                                  \
+  BOOST_PP_REPEAT(BOOST_PP_ADD(n,2),TABLE,tt)                                    \
+  BOOST_PP_CAT(algebra.for_each,BOOST_PP_ADD(n,2))(                              \
+               BOOST_PP_ENUM_PARAMS(BOOST_PP_ADD(n,2),y), default_operations::   \
+  BOOST_PP_CAT(scale_sum,BOOST_PP_ADD(n,1)) <T>(                                 \
+               BOOST_PP_REPEAT(n, PARAMS, text ) T(BOOST_PP_ADD(n,1))));         \
+  BOOST_PP_REPEAT(2,TEST,n)                                                      \
+}
+
+BOOST_AUTO_TEST_SUITE( nt2_algebra )
+
+BOOST_PP_REPEAT(7,TEST_CASE,dummy)
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/nt2/copy.cpp b/test_external/nt2/copy.cpp
new file mode 100644
index 0000000..841728e
--- /dev/null
+++ b/test_external/nt2/copy.cpp
@@ -0,0 +1,44 @@
+//==============================================================================
+//         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+//         Copyright 2014          NumScale SAS
+//
+//          Distributed under the Boost Software License, Version 1.0.
+//                 See accompanying file LICENSE.txt or copy at
+//                     http://www.boost.org/LICENSE_1_0.txt
+//==============================================================================
+#include <boost/numeric/odeint.hpp>
+#include <nt2/table.hpp>
+#include <nt2/include/functions/linspace.hpp>
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_nt2_copy
+
+#include <boost/test/included/unit_test.hpp>
+#include <boost/numeric/odeint/external/nt2/nt2_copy.hpp>
+
+#include <boost/mpl/list.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef boost::mpl::list< float , double > fp_types;
+
+BOOST_AUTO_TEST_SUITE( nt2_copy )
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( test_copy, T, fp_types )
+{
+  nt2::table<T> x = nt2::linspace(T(1),T(0),7);
+
+  nt2::table<T> y;
+
+  copy(y,x);
+
+  for (std::size_t ii=1; ii<=x.size();ii++)
+    BOOST_CHECK_EQUAL(x(ii),y(ii));
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/nt2/is_resizeable.cpp b/test_external/nt2/is_resizeable.cpp
new file mode 100644
index 0000000..10ecdbb
--- /dev/null
+++ b/test_external/nt2/is_resizeable.cpp
@@ -0,0 +1,36 @@
+//==============================================================================
+//         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+//         Copyright 2014          NumScale SAS
+//
+//          Distributed under the Boost Software License, Version 1.0.
+//                 See accompanying file LICENSE.txt or copy at
+//                     http://www.boost.org/LICENSE_1_0.txt
+//==============================================================================
+#include <boost/numeric/odeint.hpp>
+#include <nt2/table.hpp>
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_nt2_resize
+
+#include <boost/test/included/unit_test.hpp>
+#include <boost/numeric/odeint/external/nt2/nt2_resize.hpp>
+
+#include <boost/mpl/list.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef boost::mpl::list< float , double > fp_types;
+
+BOOST_AUTO_TEST_SUITE( nt2_is_resizeable )
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( is_resizeable, T, fp_types )
+{
+  BOOST_STATIC_ASSERT(( boost::numeric::odeint::is_resizeable< nt2::table<T> >::value ));
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/nt2/norm_inf.cpp b/test_external/nt2/norm_inf.cpp
new file mode 100644
index 0000000..1e12061
--- /dev/null
+++ b/test_external/nt2/norm_inf.cpp
@@ -0,0 +1,46 @@
+//==============================================================================
+//         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+//         Copyright 2014          NumScale SAS
+//
+//          Distributed under the Boost Software License, Version 1.0.
+//                 See accompanying file LICENSE.txt or copy at
+//                     http://www.boost.org/LICENSE_1_0.txt
+//==============================================================================
+#include <boost/numeric/odeint.hpp>
+#include <nt2/table.hpp>
+#include <nt2/include/functions/zeros.hpp>
+#include <nt2/include/functions/ones.hpp>
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_nt2_copy
+
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/numeric/odeint/external/nt2/nt2_norm_inf.hpp>
+
+#include <boost/mpl/list.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef boost::mpl::list< float , double > fp_types;
+
+BOOST_AUTO_TEST_SUITE( nt2_norm_inf )
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( test_norm_inf, T, fp_types )
+{
+  nt2::table<T> x = nt2::ones(10,1, nt2::meta::as_<T>() );
+  x(4) = 55;
+
+  nt2::table<T> y = nt2::zeros(8,8, nt2::meta::as_<T>() );
+  y(6,4) = -42;
+
+  BOOST_CHECK_SMALL(vector_space_norm_inf<nt2::table<T> >()(x) - T(55), T(1e-10));
+  BOOST_CHECK_SMALL(vector_space_norm_inf<nt2::table<T> >()(y) - T(42), T(1e-10));
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/nt2/resize.cpp b/test_external/nt2/resize.cpp
new file mode 100644
index 0000000..56d7ce4
--- /dev/null
+++ b/test_external/nt2/resize.cpp
@@ -0,0 +1,45 @@
+//==============================================================================
+//         Copyright 2014          LRI    UMR 8623 CNRS/Univ Paris Sud XI
+//         Copyright 2014          NumScale SAS
+//
+//          Distributed under the Boost Software License, Version 1.0.
+//                 See accompanying file LICENSE.txt or copy at
+//                     http://www.boost.org/LICENSE_1_0.txt
+//==============================================================================
+#include <boost/numeric/odeint.hpp>
+#include <nt2/table.hpp>
+
+#include <boost/config.hpp>
+#ifdef BOOST_MSVC
+    #pragma warning(disable:4996)
+#endif
+
+#define BOOST_TEST_MODULE odeint_nt2_resize
+
+#include <boost/test/included/unit_test.hpp>
+#include <boost/numeric/odeint/external/nt2/nt2_resize.hpp>
+
+#include <boost/mpl/list.hpp>
+
+using namespace boost::unit_test;
+using namespace boost::numeric::odeint;
+
+typedef boost::mpl::list< float , double > fp_types;
+
+BOOST_AUTO_TEST_SUITE( nt2_resize )
+
+BOOST_AUTO_TEST_CASE_TEMPLATE( test_resize, T, fp_types )
+{
+  nt2::table<T> x;
+  x.resize(nt2::of_size(10,10));
+
+  nt2::table<T> y;
+
+  BOOST_CHECK_EQUAL(same_size(x,y),false);
+
+  resize(y,x);
+
+  BOOST_CHECK_EQUAL(same_size(x,y),true);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test_external/thrust/Makefile b/test_external/thrust/Makefile
new file mode 100644
index 0000000..49d9cd7
--- /dev/null
+++ b/test_external/thrust/Makefile
@@ -0,0 +1,34 @@
+# Copyright 2010-2014 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)
+
+# make sure BOOST_ROOT is pointing to your boost directory
+# otherwise, set it here:
+# BOOST_ROOT = /path/to/boost
+
+# path to the cuda installation
+CUDA_ROOT = /usr/local/cuda
+# target architecture
+ARCH = sm_13
+
+NVCC = $(CUDA_ROOT)/bin/nvcc
+
+INCLUDES += -I../../include/ -I$(BOOST_ROOT)
+
+NVCCFLAGS = -O3 $(INCLUDES) -arch $(ARCH)
+
+%.o : %.cu
+	$(NVCC) $(NVCCFLAGS) -c $< -o $@
+
+% : %.o
+	$(NVCC) $(NVCCFLAGS) -o $@ $<
+
+
+all : check_thrust
+
+
+clean :
+	-rm *~ *.o check_thrust
diff --git a/test_external/thrust/check_thrust.cu b/test_external/thrust/check_thrust.cu
new file mode 100644
index 0000000..3623c37
--- /dev/null
+++ b/test_external/thrust/check_thrust.cu
@@ -0,0 +1,72 @@
+/* Boost check_thrust.cu test file
+ 
+ Copyright 2010-2013 Mario Mulansky
+ Copyright 2010-2011 Karsten Ahnert
+ 
+ This file tests the use of the euler stepper
+  
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+//#include <boost/test/unit_test.hpp>
+
+#include <boost/numeric/odeint/stepper/euler.hpp>
+#include <boost/numeric/odeint/external/thrust/thrust.hpp>
+
+#include <thrust/device_vector.h>
+#include <thrust/fill.h>
+
+using namespace boost::numeric::odeint;
+
+typedef float base_type;
+// typedef thrust::device_vector< base_type > state_type;
+typedef thrust::host_vector< base_type > state_type;
+
+void constant_system( const state_type &x , state_type &dxdt , base_type t )
+{
+	thrust::fill( dxdt.begin() , dxdt.end() , static_cast<base_type>(1.0) );
+}
+
+const base_type eps = 1.0e-7;
+
+
+template< class Stepper , class System >
+void check_stepper_concept( Stepper &stepper , System system , typename Stepper::state_type &x )
+{
+    typedef Stepper stepper_type;
+    typedef typename stepper_type::state_type container_type;
+    typedef typename stepper_type::order_type order_type;
+    typedef typename stepper_type::time_type time_type;
+
+    stepper.do_step( system , x , 0.0 , 0.1 );
+    base_type xval = *boost::begin( x );
+    if( fabs( xval - 0.1 ) < eps )
+    	std::clog << "TEST PASSED" << std::endl;
+    else
+    	std::clog << "TEST FAILED" << std::endl;
+}
+
+void test_euler_with_thrust( void )
+{
+	state_type x(1);
+	thrust::fill( x.begin() , x.end() , static_cast<base_type>(0.0) );
+	euler< state_type , base_type , state_type , base_type > euler;
+	check_stepper_concept( euler , constant_system , x );
+
+
+}
+
+/*test_suite* init_unit_test_suite( int argc, char* argv[] )
+{
+    test_suite *test = BOOST_TEST_SUITE("check stepper with thrust");
+
+    test->add( BOOST_TEST_CASE( &test_euler_with_thrust ) );
+
+    return test;
+}*/
+
+int main() {
+	test_euler_with_thrust();
+}
diff --git a/test_external/vexcl/Jamfile.v2 b/test_external/vexcl/Jamfile.v2
new file mode 100644
index 0000000..9b76b43
--- /dev/null
+++ b/test_external/vexcl/Jamfile.v2
@@ -0,0 +1,35 @@
+# Copyright 2012 Karsten Ahnert
+# Copyright 2013 Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+# bring in rules for testing
+
+
+import testing ;
+
+use-project boost : $(BOOST_ROOT) ;
+VEXCL_INCLUDE = /home/karsten/boost/testing/vexcl ;
+OPENCL_INCLUDE = /usr/local/cuda/include ;
+#OPENCL_INCLUDE = /usr/include ;
+
+project
+    : requirements
+      <library>/boost/test//boost_unit_test_framework
+      <define>BOOST_ALL_NO_LIB=1
+      <include>$(VEXCL_INCLUDE)
+      <include>$(OPENCL_INCLUDE)
+      <cxxflags>-std=c++0x
+      <library>/boost//system/
+    ;
+
+lib OpenCL : : <name>OpenCL <link>shared ;
+
+test-suite "odeint"
+    :
+    [ run lorenz.cpp OpenCL ]
+    [ run norm_inf.cpp OpenCL ]
+    : <testing.launcher>valgrind
+    : 
+    : <link>shared:<define>BOOST_TEST_DYN_LINK=1
+    ;
\ No newline at end of file
diff --git a/test_external/vexcl/lorenz.cpp b/test_external/vexcl/lorenz.cpp
new file mode 100644
index 0000000..6bde91f
--- /dev/null
+++ b/test_external/vexcl/lorenz.cpp
@@ -0,0 +1,147 @@
+/* Boost lorenz.cpp test file
+
+ Copyright 2012 Karsten Ahnert
+ Copyright 2012 Mario Mulansky
+
+ This file tests the odeint library with the vexcl types
+
+ 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)
+*/
+
+#define BOOST_TEST_MODULE odeint_vexcl
+
+#include <vector>
+#include <iostream>
+
+#include <vexcl/vexcl.hpp>
+
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
+#include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
+#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
+#include <boost/numeric/odeint/integrate/integrate_const.hpp>
+#include <boost/numeric/odeint/external/vexcl/vexcl.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+using namespace boost::unit_test;
+namespace odeint = boost::numeric::odeint;
+
+
+typedef vex::vector< double >    vector_type;
+typedef vex::multivector< double, 3 > state_type;
+
+
+
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+struct sys_func
+{
+    const vector_type &R;
+
+    sys_func( const vector_type &_R ) : R( _R ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const
+    {
+        dxdt(0) = -sigma * ( x(0) - x(1) );
+        dxdt(1) = R * x(0) - x(1) - x(0) * x(2);
+        dxdt(2) = - b * x(2) + x(0) * x(1);
+    }
+};
+
+const size_t n = 1024 ;
+const double dt = 0.01;
+const double t_max = 1.0;
+
+
+BOOST_AUTO_TEST_CASE( integrate_const_rk4 )
+{
+    vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
+
+    double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
+    std::vector<double> x( n * 3 ) , r( n );
+    for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
+
+    state_type X( ctx.queue() , n );
+    X(0) = 10.0;
+    X(1) = 10.0;
+    X(2) = 10.0;
+
+    vector_type R( ctx.queue() , r );
+
+    odeint::runge_kutta4<
+        state_type , double , state_type , double ,
+        odeint::vector_space_algebra , odeint::default_operations
+        > stepper;
+
+    odeint::integrate_const( stepper , sys_func( R ) , X , 0.0 , t_max , dt );
+
+    std::vector< double > res( 3 * n );
+    vex::copy( X(0) , res );
+    std::cout << res[0] << std::endl;
+}
+
+BOOST_AUTO_TEST_CASE( integrate_const_controlled_rk54 )
+{
+    vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
+
+    double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
+    std::vector<double> x( n * 3 ) , r( n );
+    for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
+
+    state_type X( ctx.queue() , n );
+    X(0) = 10.0;
+    X(1) = 10.0;
+    X(2) = 10.0;
+
+    vector_type R( ctx.queue() , r );
+
+    typedef odeint::runge_kutta_cash_karp54<
+        state_type , double , state_type , double ,
+        odeint::vector_space_algebra , odeint::default_operations
+        > stepper_type;
+    typedef odeint::controlled_runge_kutta< stepper_type > controlled_stepper_type;
+
+    odeint::integrate_const( controlled_stepper_type() , sys_func( R ) , X , 0.0 , t_max , dt );
+
+    std::vector< double > res( 3 * n );
+    vex::copy( X(0) , res );
+    std::cout << res[0] << std::endl;
+}
+
+BOOST_AUTO_TEST_CASE( integrate_const_dense_output_dopri5 )
+{
+    vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
+
+    double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
+    std::vector<double> x( n * 3 ) , r( n );
+    for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
+
+    state_type X( ctx.queue() , n );
+    X(0) = 10.0;
+    X(1) = 10.0;
+    X(2) = 10.0;
+
+    vector_type R( ctx.queue() , r );
+
+    typedef odeint::runge_kutta_dopri5<
+        state_type , double , state_type , double ,
+        odeint::vector_space_algebra , odeint::default_operations
+        > stepper_type;
+    typedef odeint::controlled_runge_kutta< stepper_type > controlled_stepper_type;
+    typedef odeint::dense_output_runge_kutta< controlled_stepper_type > dense_output_stepper_type;
+
+    odeint::integrate_const( dense_output_stepper_type() , sys_func( R ) , X , 0.0 , t_max , dt );
+
+    std::vector< double > res( 3 * n );
+    vex::copy( X(0) , res );
+    std::cout << res[0] << std::endl;
+}
+
+
+
diff --git a/test_external/vexcl/norm_inf.cpp b/test_external/vexcl/norm_inf.cpp
new file mode 100644
index 0000000..71481d2
--- /dev/null
+++ b/test_external/vexcl/norm_inf.cpp
@@ -0,0 +1,25 @@
+#define BOOST_TEST_MODULE odeint_vexcl_norm_inf
+
+#include <boost/numeric/odeint/external/vexcl/vexcl_norm_inf.hpp>
+#include <boost/test/unit_test.hpp>
+
+template <class T>
+double norm(const T &x) {
+    return boost::numeric::odeint::vector_space_norm_inf<T>()(x);
+}
+
+BOOST_AUTO_TEST_CASE( norm_inf )
+{
+    vex::Context ctx(vex::Filter::Env);
+    std::cout << ctx << std::endl;
+
+    vex::vector<double> x(ctx, 1024);
+    x = 41;
+
+    vex::multivector<double, 2> y(ctx, 1024);
+    y = 42;
+
+    BOOST_CHECK_EQUAL( norm(x), 41 );
+    BOOST_CHECK_EQUAL( norm(y), 42 );
+}
+