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