ceph/src/boost/libs/numeric/odeint/performance/SIMD/roessler.cpp

   1 /*
   2  * Simulation of an ensemble of Roessler attractors
   3  *
   4  * Copyright 2014 Mario Mulansky
   5  *
   6  * Distributed under the Boost Software License, Version 1.0.
   7  * (See accompanying file LICENSE_1_0.txt or
   8  * copy at http://www.boost.org/LICENSE_1_0.txt)
   9  *
  10  */
  11
  12
  13 #include <iostream>
  14 #include <vector>
  15 #include <random>
  16
  17 #include <boost/timer.hpp>
  18 #include <boost/array.hpp>
  19
  20 #include <boost/numeric/odeint.hpp>
  21
  22 namespace odeint = boost::numeric::odeint;
  23
  24 typedef boost::timer timer_type;
  25
  26 typedef double fp_type;
  27 //typedef float fp_type;
  28
  29 typedef boost::array<fp_type, 3> state_type;
  30 typedef std::vector<state_type> state_vec;
  31
  32 //---------------------------------------------------------------------------
  33 struct roessler_system {
  34     const fp_type m_a, m_b, m_c;
  35
  36     roessler_system(const fp_type a, const fp_type b, const fp_type c)
  37         : m_a(a), m_b(b), m_c(c)
  38     {}
  39
  40     void operator()(const state_type &x, state_type &dxdt, const fp_type t) const
  41     {
  42         dxdt[0] = -x[1] - x[2];
  43         dxdt[1] = x[0] + m_a * x[1];
  44         dxdt[2] = m_b + x[2] * (x[0] - m_c);
  45     }
  46 };
  47
  48 //---------------------------------------------------------------------------
  49 int main(int argc, char *argv[]) {
  50 if(argc<3)
  51 {
  52     std::cerr << "Expected size and steps as parameter" << std::endl;
  53     exit(1);
  54 }
  55 const size_t n = atoi(argv[1]);
  56 const size_t steps = atoi(argv[2]);
  57 //const size_t steps = 50;
  58
  59 const fp_type dt = 0.01;
  60
  61 const fp_type a = 0.2;
  62 const fp_type b = 1.0;
  63 const fp_type c = 9.0;
  64
  65 // random initial conditions on the device
  66 std::vector<fp_type> x(n), y(n), z(n);
  67 std::default_random_engine generator;
  68 std::uniform_real_distribution<fp_type> distribution_xy(-8.0, 8.0);
  69 std::uniform_real_distribution<fp_type> distribution_z(0.0, 20.0);
  70 auto rand_xy = std::bind(distribution_xy, std::ref(generator));
  71 auto rand_z = std::bind(distribution_z, std::ref(generator));
  72 std::generate(x.begin(), x.end(), rand_xy);
  73 std::generate(y.begin(), y.end(), rand_xy);
  74 std::generate(z.begin(), z.end(), rand_z);
  75
  76 state_vec state(n);
  77 for(size_t i=0; i<n; ++i)
  78 {
  79     state[i][0] = x[i];
  80     state[i][1] = y[i];
  81     state[i][2] = z[i];
  82 }
  83
  84 std::cout.precision(16);
  85
  86 std::cout << "# n: " << n << std::endl;
  87
  88 std::cout << x[0] << std::endl;
  89
  90
  91 // Stepper type - use never_resizer for slight performance improvement
  92 odeint::runge_kutta4_classic<state_type, fp_type, state_type, fp_type,
  93                              odeint::array_algebra,
  94                              odeint::default_operations,
  95                              odeint::never_resizer> stepper;
  96
  97 roessler_system sys(a, b, c);
  98
  99 timer_type timer;
 100
 101 fp_type t = 0.0;
 102
 103 for (int step = 0; step < steps; step++)
 104 {
 105     for(size_t i=0; i<n; ++i)
 106     {
 107         stepper.do_step(sys, state[i], t, dt);
 108     }
 109     t += dt;
 110 }
 111
 112 std::cout << "Integration finished, runtime for " << steps << " steps: ";
 113 std::cout << timer.elapsed() << " s" << std::endl;
 114
 115 // compute some accumulation to make sure all results have been computed
 116 fp_type s = 0.0;
 117 for(size_t i = 0; i < n; ++i)
 118 {
 119     s += state[i][0];
 120 }
 121
 122 std::cout << state[0][0] << std::endl;
 123 std::cout << s/n << std::endl;
 124
 125 }