diff options
Diffstat (limited to 'libs/numeric/odeint/performance')
42 files changed, 606 insertions, 2372 deletions
diff --git a/libs/numeric/odeint/performance/Jamfile.v2 b/libs/numeric/odeint/performance/Jamfile.v2 index 54b89fa34..e60e4ea12 100644 --- a/libs/numeric/odeint/performance/Jamfile.v2 +++ b/libs/numeric/odeint/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 - ; - - diff --git a/libs/numeric/odeint/performance/Makefile b/libs/numeric/odeint/performance/Makefile new file mode 100644 index 000000000..641cb0313 --- /dev/null +++ b/libs/numeric/odeint/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 diff --git a/libs/numeric/odeint/performance/SIMD/Makefile b/libs/numeric/odeint/performance/SIMD/Makefile new file mode 100644 index 000000000..811acd988 --- /dev/null +++ b/libs/numeric/odeint/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} diff --git a/libs/numeric/odeint/performance/SIMD/perf_roessler.sh b/libs/numeric/odeint/performance/SIMD/perf_roessler.sh new file mode 100755 index 000000000..a1094f63a --- /dev/null +++ b/libs/numeric/odeint/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 diff --git a/libs/numeric/odeint/performance/SIMD/roessler.cpp b/libs/numeric/odeint/performance/SIMD/roessler.cpp new file mode 100644 index 000000000..4e6cc4229 --- /dev/null +++ b/libs/numeric/odeint/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; + +} diff --git a/libs/numeric/odeint/performance/SIMD/roessler_simd.cpp b/libs/numeric/odeint/performance/SIMD/roessler_simd.cpp new file mode 100644 index 000000000..d79af4d8b --- /dev/null +++ b/libs/numeric/odeint/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; + +} diff --git a/libs/numeric/odeint/performance/c_lorenz.c b/libs/numeric/odeint/performance/c_lorenz.c new file mode 100644 index 000000000..85aba7fde --- /dev/null +++ b/libs/numeric/odeint/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; +} diff --git a/libs/numeric/odeint/performance/fortran_lorenz.f90 b/libs/numeric/odeint/performance/fortran_lorenz.f90 new file mode 100644 index 000000000..26869973c --- /dev/null +++ b/libs/numeric/odeint/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 diff --git a/libs/numeric/odeint/performance/fusion_algebra.hpp b/libs/numeric/odeint/performance/fusion_algebra.hpp deleted file mode 100644 index 1233ce006..000000000 --- a/libs/numeric/odeint/performance/fusion_algebra.hpp +++ /dev/null @@ -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_ */ diff --git a/libs/numeric/odeint/performance/fusion_explicit_error_rk.hpp b/libs/numeric/odeint/performance/fusion_explicit_error_rk.hpp deleted file mode 100644 index b0085d641..000000000 --- a/libs/numeric/odeint/performance/fusion_explicit_error_rk.hpp +++ /dev/null @@ -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_ */ diff --git a/libs/numeric/odeint/performance/fusion_explicit_rk_new.hpp b/libs/numeric/odeint/performance/fusion_explicit_rk_new.hpp deleted file mode 100644 index d16a67b5e..000000000 --- a/libs/numeric/odeint/performance/fusion_explicit_rk_new.hpp +++ /dev/null @@ -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_ */ diff --git a/libs/numeric/odeint/performance/generic_odeint_rk4_lorenz.cpp b/libs/numeric/odeint/performance/generic_odeint_rk4_lorenz.cpp deleted file mode 100644 index 0d0975646..000000000 --- a/libs/numeric/odeint/performance/generic_odeint_rk4_lorenz.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/gsl_rk4_lorenz.cpp b/libs/numeric/odeint/performance/gsl_rk4_lorenz.cpp deleted file mode 100644 index f5f13882a..000000000 --- a/libs/numeric/odeint/performance/gsl_rk4_lorenz.cpp +++ /dev/null @@ -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); -} diff --git a/libs/numeric/odeint/performance/lorenz.hpp b/libs/numeric/odeint/performance/lorenz.hpp index 22de81361..c1ea37c9e 100644 --- a/libs/numeric/odeint/performance/lorenz.hpp +++ b/libs/numeric/odeint/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_ */ diff --git a/libs/numeric/odeint/performance/lorenz_gsl.hpp b/libs/numeric/odeint/performance/lorenz_gsl.hpp deleted file mode 100644 index a907aed87..000000000 --- a/libs/numeric/odeint/performance/lorenz_gsl.hpp +++ /dev/null @@ -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 diff --git a/libs/numeric/odeint/performance/mpi/Jamfile.v2 b/libs/numeric/odeint/performance/mpi/Jamfile.v2 deleted file mode 100644 index d80e860d8..000000000 --- a/libs/numeric/odeint/performance/mpi/Jamfile.v2 +++ /dev/null @@ -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 ; diff --git a/libs/numeric/odeint/performance/mpi/osc_chain_1d.cpp b/libs/numeric/odeint/performance/mpi/osc_chain_1d.cpp deleted file mode 100644 index ef2ce7dbd..000000000 --- a/libs/numeric/odeint/performance/mpi/osc_chain_1d.cpp +++ /dev/null @@ -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; -} - diff --git a/libs/numeric/odeint/performance/mpi/osc_chain_1d_system.hpp b/libs/numeric/odeint/performance/mpi/osc_chain_1d_system.hpp deleted file mode 100644 index 697208458..000000000 --- a/libs/numeric/odeint/performance/mpi/osc_chain_1d_system.hpp +++ /dev/null @@ -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 diff --git a/libs/numeric/odeint/performance/mpi/osc_chain_speedup.gnu b/libs/numeric/odeint/performance/mpi/osc_chain_speedup.gnu deleted file mode 100755 index 7f3877257..000000000 --- a/libs/numeric/odeint/performance/mpi/osc_chain_speedup.gnu +++ /dev/null @@ -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 diff --git a/libs/numeric/odeint/performance/mpi/osc_chain_speedup.sh b/libs/numeric/odeint/performance/mpi/osc_chain_speedup.sh deleted file mode 100755 index 7a0ff7891..000000000 --- a/libs/numeric/odeint/performance/mpi/osc_chain_speedup.sh +++ /dev/null @@ -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 diff --git a/libs/numeric/odeint/performance/nr_rk4_lorenz.cpp b/libs/numeric/odeint/performance/nr_rk4_lorenz.cpp deleted file mode 100644 index 36be91df8..000000000 --- a/libs/numeric/odeint/performance/nr_rk4_lorenz.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/nr_rk4_phase_lattice.cpp b/libs/numeric/odeint/performance/nr_rk4_phase_lattice.cpp deleted file mode 100644 index ca0d4419a..000000000 --- a/libs/numeric/odeint/performance/nr_rk4_phase_lattice.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/odeint_rk4_array.cpp b/libs/numeric/odeint/performance/odeint_rk4_array.cpp new file mode 100644 index 000000000..6d60296f2 --- /dev/null +++ b/libs/numeric/odeint/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; +} diff --git a/libs/numeric/odeint/performance/odeint_rk4_lorenz_array.cpp b/libs/numeric/odeint/performance/odeint_rk4_lorenz_array.cpp deleted file mode 100644 index 5d73a4d63..000000000 --- a/libs/numeric/odeint/performance/odeint_rk4_lorenz_array.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/odeint_rk4_lorenz_range.cpp b/libs/numeric/odeint/performance/odeint_rk4_lorenz_range.cpp deleted file mode 100644 index 0e7172416..000000000 --- a/libs/numeric/odeint/performance/odeint_rk4_lorenz_range.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/odeint_rk4_phase_lattice.cpp b/libs/numeric/odeint/performance/odeint_rk4_phase_lattice.cpp deleted file mode 100644 index 1703b5c8a..000000000 --- a/libs/numeric/odeint/performance/odeint_rk4_phase_lattice.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/odeint_rk4_phase_lattice_mkl.cpp b/libs/numeric/odeint/performance/odeint_rk4_phase_lattice_mkl.cpp deleted file mode 100644 index c900c0e4c..000000000 --- a/libs/numeric/odeint/performance/odeint_rk4_phase_lattice_mkl.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/openmp/Jamfile.v2 b/libs/numeric/odeint/performance/openmp/Jamfile.v2 deleted file mode 100644 index 098cd53f0..000000000 --- a/libs/numeric/odeint/performance/openmp/Jamfile.v2 +++ /dev/null @@ -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 ; diff --git a/libs/numeric/odeint/performance/openmp/osc_chain_1d.cpp b/libs/numeric/odeint/performance/openmp/osc_chain_1d.cpp deleted file mode 100644 index 8a5a2e857..000000000 --- a/libs/numeric/odeint/performance/openmp/osc_chain_1d.cpp +++ /dev/null @@ -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; -} diff --git a/libs/numeric/odeint/performance/openmp/osc_chain_1d_system.hpp b/libs/numeric/odeint/performance/openmp/osc_chain_1d_system.hpp deleted file mode 100644 index d9a6c222d..000000000 --- a/libs/numeric/odeint/performance/openmp/osc_chain_1d_system.hpp +++ /dev/null @@ -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 diff --git a/libs/numeric/odeint/performance/openmp/osc_chain_speedup.gnu b/libs/numeric/odeint/performance/openmp/osc_chain_speedup.gnu deleted file mode 100755 index 4d6bc95e4..000000000 --- a/libs/numeric/odeint/performance/openmp/osc_chain_speedup.gnu +++ /dev/null @@ -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 diff --git a/libs/numeric/odeint/performance/openmp/osc_chain_speedup.sh b/libs/numeric/odeint/performance/openmp/osc_chain_speedup.sh deleted file mode 100755 index 9d4c91096..000000000 --- a/libs/numeric/odeint/performance/openmp/osc_chain_speedup.sh +++ /dev/null @@ -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 diff --git a/libs/numeric/odeint/performance/performance.py b/libs/numeric/odeint/performance/performance.py deleted file mode 100644 index a51ce32b0..000000000 --- a/libs/numeric/odeint/performance/performance.py +++ /dev/null @@ -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() -""" diff --git a/libs/numeric/odeint/performance/phase_lattice.hpp b/libs/numeric/odeint/performance/phase_lattice.hpp deleted file mode 100644 index 6657edeea..000000000 --- a/libs/numeric/odeint/performance/phase_lattice.hpp +++ /dev/null @@ -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] ); - } - -}; diff --git a/libs/numeric/odeint/performance/phase_lattice_mkl.hpp b/libs/numeric/odeint/performance/phase_lattice_mkl.hpp deleted file mode 100644 index 4feea47c1..000000000 --- a/libs/numeric/odeint/performance/phase_lattice_mkl.hpp +++ /dev/null @@ -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_ */ diff --git a/libs/numeric/odeint/performance/plot_result.py b/libs/numeric/odeint/performance/plot_result.py index f0968e51a..f39e49fce 100644 --- a/libs/numeric/odeint/performance/plot_result.py +++ b/libs/numeric/odeint/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 + +plt.rc("font", size=16) + + +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]) + + +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")] + +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")] -from pylab import * +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")] -#toolset = "gcc-4.5" -toolset = "intel-11.1" -#toolset = "msvc" -#toolset = "msvc-10.0" +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")] -#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 = "" +print t_c_intel -res = loadtxt( bin_path + "rk4_lorenz.dat" ) -res = 100*res[0]/res +ind = np.arange(3) # the x locations for the groups +width = 0.15 # the width of the bars -bar_width = 0.6 +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") -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) +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}) -savefig( bin_path + "rk4_lorenz.png" ) +plt.savefig("perf.pdf") +plt.savefig("perf.png", dpi=50) -show() +plt.show() diff --git a/libs/numeric/odeint/performance/rk4_lorenz.f b/libs/numeric/odeint/performance/rk4_lorenz.f deleted file mode 100644 index d47e53298..000000000 --- a/libs/numeric/odeint/performance/rk4_lorenz.f +++ /dev/null @@ -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 - diff --git a/libs/numeric/odeint/performance/rk_performance_test_case.hpp b/libs/numeric/odeint/performance/rk_performance_test_case.hpp deleted file mode 100644 index fa9988b7d..000000000 --- a/libs/numeric/odeint/performance/rk_performance_test_case.hpp +++ /dev/null @@ -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; -} diff --git a/libs/numeric/odeint/performance/rt_algebra.hpp b/libs/numeric/odeint/performance/rt_algebra.hpp deleted file mode 100644 index e4805aac8..000000000 --- a/libs/numeric/odeint/performance/rt_algebra.hpp +++ /dev/null @@ -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]; - } - } -}; diff --git a/libs/numeric/odeint/performance/rt_explicit_rk.hpp b/libs/numeric/odeint/performance/rt_explicit_rk.hpp deleted file mode 100644 index 3ac05f562..000000000 --- a/libs/numeric/odeint/performance/rt_explicit_rk.hpp +++ /dev/null @@ -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_ */ diff --git a/libs/numeric/odeint/performance/rt_generic_rk4_lorenz.cpp b/libs/numeric/odeint/performance/rt_generic_rk4_lorenz.cpp deleted file mode 100644 index ce4e93b55..000000000 --- a/libs/numeric/odeint/performance/rt_generic_rk4_lorenz.cpp +++ /dev/null @@ -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 ); -} diff --git a/libs/numeric/odeint/performance/rt_generic_rk4_phase_lattice.cpp b/libs/numeric/odeint/performance/rt_generic_rk4_phase_lattice.cpp deleted file mode 100644 index 0464f6914..000000000 --- a/libs/numeric/odeint/performance/rt_generic_rk4_phase_lattice.cpp +++ /dev/null @@ -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 ); -} |