1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
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;
}
|
