1 /*
2 * Copyright 2012 Karsten Ahnert
3 * Copyright 2012 Mario Mulansky
4 *
5 * Distributed under the Boost Software License, Version 1.0.
6 * (See accompanying file LICENSE_1_0.txt or
7 * copy at http://www.boost.org/LICENSE_1_0.txt)
8 */
9
10
11 #include <iostream>
12 #include <fstream>
13 #include <utility>
14 #include "time.h"
15
16 #include <boost/numeric/odeint.hpp>
17 #include <boost/phoenix/phoenix.hpp>
18 #include <boost/numeric/mtl/mtl.hpp>
19
20 #include <boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp>
21
22 using namespace std;
23 using namespace boost::numeric::odeint;
24
25 namespace phoenix = boost::phoenix;
26
27
28
29 typedef mtl::dense_vector< double > vec_mtl4;
30 typedef mtl::compressed2D< double > mat_mtl4;
31
32 typedef boost::numeric::ublas::vector< double > vec_ublas;
33 typedef boost::numeric::ublas::matrix< double > mat_ublas;
34
35
36 // two systems defined 1 & 2 both are mostly sparse with the number of element variable
37 struct system1_mtl4
38 {
39
operator ()system1_mtl440 void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
41 {
42 int size = mtl::size(x);
43
44 dxdt[ 0 ] = -0.06*x[0];
45
46 for (int i =1; i< size ; ++i){
47
48 dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
49 }
50
51 }
52 };
53
54 struct jacobi1_mtl4
55 {
operator ()jacobi1_mtl456 void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
57 {
58 int size = mtl::size(x);
59 mtl::matrix::inserter<mat_mtl4> ins(J);
60
61 ins[0][0]=-0.06;
62
63 for (int i =1; i< size ; ++i)
64 {
65 ins[i][i-1] = + 4.2;
66 ins[i][i] = -4.2*x[i];
67 }
68 }
69 };
70
71
72
73 struct system1_ublas
74 {
75
operator ()system1_ublas76 void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
77 {
78 int size = x.size();
79
80 dxdt[ 0 ] = -0.06*x[0];
81
82 for (int i =1; i< size ; ++i){
83
84 dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
85 }
86
87 }
88 };
89
90 struct jacobi1_ublas
91 {
operator ()jacobi1_ublas92 void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
93 {
94 int size = x.size();
95 // mtl::matrix::inserter<mat_mtl4> ins(J);
96
97 J(0,0)=-0.06;
98
99 for (int i =1; i< size ; ++i){
100 //ins[i][0]=120.0*x[i];
101 J(i,i-1) = + 4.2;
102 J(i,i) = -4.2*x[i];
103
104 }
105 }
106 };
107
108 struct system2_mtl4
109 {
110
operator ()system2_mtl4111 void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
112 {
113 int size = mtl::size(x);
114
115
116 for (int i =0; i< size/5 ; i+=5){
117
118 dxdt[ i ] = -0.5*x[i];
119 dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
120 dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
121 dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
122 dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
123
124 }
125
126 }
127 };
128
129 struct jacobi2_mtl4
130 {
operator ()jacobi2_mtl4131 void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
132 {
133 int size = mtl::size(x);
134 mtl::matrix::inserter<mat_mtl4> ins(J);
135
136 for (int i =0; i< size/5 ; i+=5){
137
138 ins[ i ][i] = -0.5;
139 ins[i+1][i+1]=25*x[i+2];
140 ins[i+1][i+2] = 25*x[i+1];
141 ins[i+1][i+3] = -740*2*x[i+3];
142 ins[i+1][i] =+4.2e-2;
143
144 ins[i+2][i]= 50*x[i];
145 ins[i+2][i+3]= -740*2*x[i+3];
146 ins[i+3][i+1] = -25*x[i+2];
147 ins[i+3][i+2] = -25*x[i+1];
148 ins[i+3][i+3] = +740*2*x[i+3];
149 ins[i+4][i] = 0.25*x[i+1];
150 ins[i+4][i+1] =0.25*x[i];
151 ins[i+4][i+3]=-44.5;
152
153
154
155 }
156 }
157 };
158
159
160
161 struct system2_ublas
162 {
163
operator ()system2_ublas164 void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
165 {
166 int size = x.size();
167 for (int i =0; i< size/5 ; i+=5){
168
169 dxdt[ i ] = -4.2e-2*x[i];
170 dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
171 dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
172 dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
173 dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
174
175 }
176
177 }
178 };
179
180 struct jacobi2_ublas
181 {
operator ()jacobi2_ublas182 void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
183 {
184 int size = x.size();
185
186 for (int i =0; i< size/5 ; i+=5){
187
188 J(i ,i) = -4.2e-2;
189 J(i+1,i+1)=25*x[i+2];
190 J(i+1,i+2) = 25*x[i+1];
191 J(i+1,i+3) = -740*2*x[i+3];
192 J(i+1,i) =+4.2e-2;
193
194 J(i+2,i)= 50*x[i];
195 J(i+2,i+3)= -740*2*x[i+3];
196 J(i+3,i+1) = -25*x[i+2];
197 J(i+3,i+2) = -25*x[i+1];
198 J(i+3,i+3) = +740*2*x[i+3];
199 J(i+4,i) = 0.25*x[i+1];
200 J(i+4,i+1) =0.25*x[i];
201 J(i+4,i+3)=-44.5;
202
203
204
205 }
206
207
208 }
209 };
210
211
212
213
testRidiculouslyMassiveArray(int size)214 void testRidiculouslyMassiveArray( int size )
215 {
216 typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
217 typedef boost::numeric::odeint::implicit_euler< double > booststepper;
218
219 vec_mtl4 x(size , 0.0);
220 x[0]=1;
221
222
223 double dt = 0.02;
224 double endtime = 10.0;
225
226 clock_t tstart_mtl4 = clock();
227 size_t num_of_steps_mtl4 = integrate_const(
228 mtl4stepper() ,
229 make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
230 x , 0.0 , endtime , dt );
231 clock_t tend_mtl4 = clock() ;
232
233 clog << x[0] << endl;
234 clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
235
236 vec_ublas x_ublas(size , 0.0);
237 x_ublas[0]=1;
238
239 clock_t tstart_boost = clock();
240 size_t num_of_steps_ublas = integrate_const(
241 booststepper() ,
242 make_pair( system1_ublas() , jacobi1_ublas() ) ,
243 x_ublas , 0.0 , endtime , dt );
244 clock_t tend_boost = clock() ;
245
246 clog << x_ublas[0] << endl;
247 clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
248
249 clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
250 return ;
251 }
252
253
254
testRidiculouslyMassiveArray2(int size)255 void testRidiculouslyMassiveArray2( int size )
256 {
257 typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
258 typedef boost::numeric::odeint::implicit_euler< double > booststepper;
259
260
261 vec_mtl4 x(size , 0.0);
262 x[0]=100;
263
264
265 double dt = 0.01;
266 double endtime = 10.0;
267
268 clock_t tstart_mtl4 = clock();
269 size_t num_of_steps_mtl4 = integrate_const(
270 mtl4stepper() ,
271 make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
272 x , 0.0 , endtime , dt );
273
274
275 clock_t tend_mtl4 = clock() ;
276
277 clog << x[0] << endl;
278 clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
279
280 vec_ublas x_ublas(size , 0.0);
281 x_ublas[0]=100;
282
283 clock_t tstart_boost = clock();
284 size_t num_of_steps_ublas = integrate_const(
285 booststepper() ,
286 make_pair( system1_ublas() , jacobi1_ublas() ) ,
287 x_ublas , 0.0 , endtime , dt );
288
289
290 clock_t tend_boost = clock() ;
291
292 clog << x_ublas[0] << endl;
293 clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
294
295 clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
296 return ;
297 }
298
299
300
301
main(int argc,char ** argv)302 int main( int argc , char **argv )
303 {
304 std::vector< size_t > length;
305 length.push_back( 8 );
306 length.push_back( 16 );
307 length.push_back( 32 );
308 length.push_back( 64 );
309 length.push_back( 128 );
310 length.push_back( 256 );
311
312 for( size_t i=0 ; i<length.size() ; ++i )
313 {
314 clog << "Testing with size " << length[i] << endl;
315 testRidiculouslyMassiveArray( length[i] );
316 }
317 clog << endl << endl;
318
319 for( size_t i=0 ; i<length.size() ; ++i )
320 {
321 clog << "Testing with size " << length[i] << endl;
322 testRidiculouslyMassiveArray2( length[i] );
323 }
324 }
325