summaryrefslogtreecommitdiff
path: root/libs/numeric/odeint/performance
diff options
context:
space:
mode:
Diffstat (limited to 'libs/numeric/odeint/performance')
-rw-r--r--libs/numeric/odeint/performance/Jamfile.v250
-rw-r--r--libs/numeric/odeint/performance/Makefile43
-rw-r--r--libs/numeric/odeint/performance/SIMD/Makefile33
-rwxr-xr-xlibs/numeric/odeint/performance/SIMD/perf_roessler.sh22
-rw-r--r--libs/numeric/odeint/performance/SIMD/roessler.cpp125
-rw-r--r--libs/numeric/odeint/performance/SIMD/roessler_simd.cpp149
-rw-r--r--libs/numeric/odeint/performance/c_lorenz.c57
-rw-r--r--libs/numeric/odeint/performance/fortran_lorenz.f9060
-rw-r--r--libs/numeric/odeint/performance/fusion_algebra.hpp202
-rw-r--r--libs/numeric/odeint/performance/fusion_explicit_error_rk.hpp63
-rw-r--r--libs/numeric/odeint/performance/fusion_explicit_rk_new.hpp217
-rw-r--r--libs/numeric/odeint/performance/generic_odeint_rk4_lorenz.cpp84
-rw-r--r--libs/numeric/odeint/performance/gsl_rk4_lorenz.cpp71
-rw-r--r--libs/numeric/odeint/performance/lorenz.hpp13
-rw-r--r--libs/numeric/odeint/performance/lorenz_gsl.hpp30
-rw-r--r--libs/numeric/odeint/performance/mpi/Jamfile.v220
-rw-r--r--libs/numeric/odeint/performance/mpi/osc_chain_1d.cpp112
-rw-r--r--libs/numeric/odeint/performance/mpi/osc_chain_1d_system.hpp88
-rwxr-xr-xlibs/numeric/odeint/performance/mpi/osc_chain_speedup.gnu38
-rwxr-xr-xlibs/numeric/odeint/performance/mpi/osc_chain_speedup.sh24
-rw-r--r--libs/numeric/odeint/performance/nr_rk4_lorenz.cpp84
-rw-r--r--libs/numeric/odeint/performance/nr_rk4_phase_lattice.cpp85
-rw-r--r--libs/numeric/odeint/performance/odeint_rk4_array.cpp63
-rw-r--r--libs/numeric/odeint/performance/odeint_rk4_lorenz_array.cpp71
-rw-r--r--libs/numeric/odeint/performance/odeint_rk4_lorenz_range.cpp59
-rw-r--r--libs/numeric/odeint/performance/odeint_rk4_phase_lattice.cpp64
-rw-r--r--libs/numeric/odeint/performance/odeint_rk4_phase_lattice_mkl.cpp69
-rw-r--r--libs/numeric/odeint/performance/openmp/Jamfile.v221
-rw-r--r--libs/numeric/odeint/performance/openmp/osc_chain_1d.cpp127
-rw-r--r--libs/numeric/odeint/performance/openmp/osc_chain_1d_system.hpp98
-rwxr-xr-xlibs/numeric/odeint/performance/openmp/osc_chain_speedup.gnu50
-rwxr-xr-xlibs/numeric/odeint/performance/openmp/osc_chain_speedup.sh38
-rw-r--r--libs/numeric/odeint/performance/performance.py70
-rw-r--r--libs/numeric/odeint/performance/phase_lattice.hpp43
-rw-r--r--libs/numeric/odeint/performance/phase_lattice_mkl.hpp57
-rw-r--r--libs/numeric/odeint/performance/plot_result.py74
-rw-r--r--libs/numeric/odeint/performance/rk4_lorenz.f53
-rw-r--r--libs/numeric/odeint/performance/rk_performance_test_case.hpp68
-rw-r--r--libs/numeric/odeint/performance/rt_algebra.hpp32
-rw-r--r--libs/numeric/odeint/performance/rt_explicit_rk.hpp87
-rw-r--r--libs/numeric/odeint/performance/rt_generic_rk4_lorenz.cpp81
-rw-r--r--libs/numeric/odeint/performance/rt_generic_rk4_phase_lattice.cpp83
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 );
-}