Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
fad_expr_depth.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Sacado Package
5// Copyright (2006) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25// (etphipp@sandia.gov).
26//
27// ***********************************************************************
28// @HEADER
29
30#include "Sacado_Random.hpp"
31
32#include "Teuchos_Time.hpp"
33#include "Teuchos_Array.hpp"
34#include <fstream>
35
36#include "fad_expr_funcs.hpp"
37
38// A simple performance test that computes the derivative of expressions of
39// various depths.
40
41template <typename T, typename F>
42double do_time(const T x[], int nloop, const F& f)
43{
44 T y = 0.0;
45 Teuchos::Time timer("F", false);
46
47 timer.start(true);
48 for (int j=0; j<nloop; j++)
49 f(x, y);
50 timer.stop();
51 return timer.totalElapsedTime() / nloop;
52}
53
54template <typename T, template <typename, int> class F>
55void do_times(const T x[], int nloop, Teuchos::Array<double>& times)
56{
57 int i = 0;
58 times[i++] = do_time(x, nloop, F<T,1>());
59 times[i++] = do_time(x, nloop, F<T,2>());
60 times[i++] = do_time(x, nloop, F<T,3>());
61 times[i++] = do_time(x, nloop, F<T,4>());
62 times[i++] = do_time(x, nloop, F<T,5>());
63 times[i++] = do_time(x, nloop, F<T,10>());
64 times[i++] = do_time(x, nloop, F<T,15>());
65 times[i++] = do_time(x, nloop, F<T,20>());
66}
67
68#ifdef HAVE_ADOLC
69template <typename F>
70double do_time_adolc(double *x, double **seed, int d, int nloop,
71 int tag, const F& f)
72{
73 Teuchos::Time timer("F", false);
74 int n = F::n;
75 trace_on(tag);
76 adouble *x_ad = new adouble[n];
77 for (int i=0; i<n; i++)
78 x_ad[i] <<= x[i];
79 adouble y_ad = 0.0;
80 f(x_ad, y_ad);
81 double y;
82 y_ad >>= y;
83 delete [] x_ad;
84 trace_off();
85
86 double **jac = new double*[1];
87 jac[0] = new double[d];
88 timer.start(true);
89 for (int j=0; j<nloop; j++)
90 forward(tag, 1, n, d, x, seed, &y, jac);
91 timer.stop();
92
93 delete [] jac[0];
94 delete [] jac;
95
96 return timer.totalElapsedTime() / nloop;
97}
98
99template <template <typename,int> class F>
100void do_times_adolc(double *x, double **seed, int d, int nloop,
101 int& tag, Teuchos::Array<double>& times)
102{
103 int i = 0;
104 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,1>());
105 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,2>());
106 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,3>());
107 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,4>());
108 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,5>());
109 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,10>());
110 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,15>());
111 times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,20>());
112}
113#endif
114
115template <typename FadType>
116void print_times(const std::string& screen_name, const std::string& file_name)
117{
118 const int nderiv = 10;
119 int deriv_dim[nderiv] = { 0, 1, 3, 5, 10, 15, 20, 30, 40, 50 };
120 Sacado::Random<double> urand(0.0, 1.0);
121 int p = 1;
122 int w = p+7;
123 std::ofstream file(file_name.c_str(), std::ios::out);
124
125 {
126 std::cout.setf(std::ios::scientific);
127 std::cout.precision(p);
128 std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
129 << std::endl;
130 std::cout << std::setw(5) << "deriv" << " ";
131 for (int i=0; i<ExprFuncs::nfunc; i++)
132 std::cout << std::setw(w) << ExprFuncs::mult_names[i] << " ";
133 std::cout << std::endl;
134 std::cout << "===== ";
135 for (int i=0; i<ExprFuncs::nfunc; i++) {
136 for (int j=0; j<w; j++)
137 std::cout << '=';
138 std::cout << " ";
139 }
140 std::cout << std::endl;
141
142 // Get function evaluation times
143 double x[ExprFuncs::nx_max];
144 for (int i=0; i<ExprFuncs::nx_max; i++)
145 x[i] = urand.number();
146 int nloop_func = 10000000;
147 Teuchos::Array<double> times_func(ExprFuncs::nfunc);
148 do_times<double,ExprFuncs::mult>(x, nloop_func, times_func);
149
150 // Get times for each derivative dimension
151 for (int i=0; i<nderiv; i++) {
153 for (int k=0; k<ExprFuncs::nx_max; k++) {
154 fx[k] = FadType(deriv_dim[i], urand.number());
155 for (int j=0; j<deriv_dim[i]; j++) {
156 fx[k].fastAccessDx(j) = urand.number();
157 }
158 }
159
160 //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
161 int nloop = static_cast<int>(100000.0);
162 Teuchos::Array<double> times(ExprFuncs::nfunc);
163 do_times<FadType,ExprFuncs::mult>(fx, nloop, times);
164
165 // Print times
166 int d = deriv_dim[i];
167 if (d == 0)
168 d = 1;
169 std::cout << std::setw(5) << deriv_dim[i] << " ";
170 file << deriv_dim[i] << " ";
171 for (int j=0; j<times.size(); j++) {
172 double rel_time = times[j]/(times_func[j]*d);
173 std::cout << std::setw(w) << rel_time << " ";
174 file << rel_time << " ";
175 }
176 std::cout << std::endl;
177 file << std::endl;
178 }
179 }
180
181 {
182 std::cout.setf(std::ios::scientific);
183 std::cout.precision(p);
184 std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
185 << std::endl;
186 std::cout << std::setw(5) << "deriv" << " ";
187 for (int i=0; i<ExprFuncs::nfunc; i++)
188 std::cout << std::setw(w) << ExprFuncs::add_names[i] << " ";
189 std::cout << std::endl;
190 std::cout << "===== ";
191 for (int i=0; i<ExprFuncs::nfunc; i++) {
192 for (int j=0; j<w; j++)
193 std::cout << '=';
194 std::cout << " ";
195 }
196 std::cout << std::endl;
197
198 // Get function evaluation times
199 double x[ExprFuncs::nx_max];
200 for (int i=0; i<ExprFuncs::nx_max; i++)
201 x[i] = urand.number();
202 int nloop_func = 10000000;
203 Teuchos::Array<double> times_func(ExprFuncs::nfunc);
204 do_times<double,ExprFuncs::add>(x, nloop_func, times_func);
205
206 // Get times for each derivative dimension
207 for (int i=0; i<nderiv; i++) {
209 for (int k=0; k<ExprFuncs::nx_max; k++) {
210 fx[k] = FadType(deriv_dim[i], urand.number());
211 for (int j=0; j<deriv_dim[i]; j++) {
212 fx[k].fastAccessDx(j) = urand.number();
213 }
214 }
215
216 //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
217 int nloop = static_cast<int>(100000.0);
218 Teuchos::Array<double> times(ExprFuncs::nfunc);
219 do_times<FadType,ExprFuncs::add>(fx, nloop, times);
220
221 // Print times
222 int d = deriv_dim[i];
223 if (d == 0)
224 d = 1;
225 std::cout << std::setw(5) << deriv_dim[i] << " ";
226 file << deriv_dim[i] << " ";
227 for (int j=0; j<times.size(); j++) {
228 double rel_time = times[j]/(times_func[j]*d);
229 std::cout << std::setw(w) << rel_time << " ";
230 file << rel_time << " ";
231 }
232 std::cout << std::endl;
233 file << std::endl;
234 }
235 }
236
237 {
238 std::cout.setf(std::ios::scientific);
239 std::cout.precision(p);
240 std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
241 << std::endl;
242 std::cout << std::setw(5) << "deriv" << " ";
243 for (int i=0; i<ExprFuncs::nfunc; i++)
244 std::cout << std::setw(w) << ExprFuncs::nest_names[i] << " ";
245 std::cout << std::endl;
246 std::cout << "===== ";
247 for (int i=0; i<ExprFuncs::nfunc; i++) {
248 for (int j=0; j<w; j++)
249 std::cout << '=';
250 std::cout << " ";
251 }
252 std::cout << std::endl;
253
254 // Get function evaluation times
255 double x[ExprFuncs::nx_max];
256 for (int i=0; i<ExprFuncs::nx_max; i++)
257 x[i] = urand.number();
258 int nloop_func = 10000000;
259 Teuchos::Array<double> times_func(ExprFuncs::nfunc);
260 do_times<double,ExprFuncs::nest>(x, nloop_func, times_func);
261
262 // Get times for each derivative dimension
263 for (int i=0; i<nderiv; i++) {
265 for (int k=0; k<ExprFuncs::nx_max; k++) {
266 fx[k] = FadType(deriv_dim[i], urand.number());
267 for (int j=0; j<deriv_dim[i]; j++) {
268 fx[k].fastAccessDx(j) = urand.number();
269 }
270 }
271
272 //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
273 int nloop = static_cast<int>(100000.0);
274 Teuchos::Array<double> times(ExprFuncs::nfunc);
275 do_times<FadType,ExprFuncs::nest>(fx, nloop, times);
276
277 // Print times
278 int d = deriv_dim[i];
279 if (d == 0)
280 d = 1;
281 std::cout << std::setw(5) << deriv_dim[i] << " ";
282 file << deriv_dim[i] << " ";
283 for (int j=0; j<times.size(); j++) {
284 double rel_time = times[j]/(times_func[j]*d);
285 std::cout << std::setw(w) << rel_time << " ";
286 file << rel_time << " ";
287 }
288 std::cout << std::endl;
289 file << std::endl;
290 }
291 }
292}
293
294#ifdef HAVE_ADOLC
295void print_times_adolc(const std::string& screen_name,
296 const std::string& file_name)
297{
298 const int nderiv = 10;
299 int deriv_dim[nderiv] = { 0, 1, 3, 5, 10, 15, 20, 30, 40, 50 };
300 const int deriv_max = 50;
301 Sacado::Random<double> urand(0.0, 1.0);
302 int p = 1;
303 int w = p+7;
304 std::ofstream file(file_name.c_str(), std::ios::out);
305 int tag = 0;
306
307 double **seed = new double*[ExprFuncs::nx_max];
308 for (int i=0; i<ExprFuncs::nx_max; i++) {
309 seed[i] = new double[deriv_max];
310 for (int j=0; j<deriv_max; j++)
311 seed[i][j] = urand.number();
312 }
313
314 {
315 std::cout.setf(std::ios::scientific);
316 std::cout.precision(p);
317 std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
318 << std::endl;
319 std::cout << std::setw(5) << "deriv" << " ";
320 for (int i=0; i<ExprFuncs::nfunc; i++)
321 std::cout << std::setw(w) << ExprFuncs::mult_names[i] << " ";
322 std::cout << std::endl;
323 std::cout << "===== ";
324 for (int i=0; i<ExprFuncs::nfunc; i++) {
325 for (int j=0; j<w; j++)
326 std::cout << '=';
327 std::cout << " ";
328 }
329 std::cout << std::endl;
330
331 // Get function evaluation times
332 double x[ExprFuncs::nx_max];
333 for (int i=0; i<ExprFuncs::nx_max; i++)
334 x[i] = urand.number();
335 int nloop_func = 10000000;
336 Teuchos::Array<double> times_func(ExprFuncs::nfunc);
337 do_times<double,ExprFuncs::mult>(x, nloop_func, times_func);
338
339 // Get times for each derivative dimension
340 for (int i=0; i<nderiv; i++) {
341 //int nloop = static_cast<int>(100000.0/(deriv_dim[i]+1));
342 int nloop = static_cast<int>(10000.0);
343 Teuchos::Array<double> times(ExprFuncs::nfunc);
344 do_times_adolc<ExprFuncs::mult>(x, seed, deriv_dim[i], nloop, tag, times);
345
346 // Print times
347 int d = deriv_dim[i];
348 if (d == 0)
349 d = 1;
350 std::cout << std::setw(5) << deriv_dim[i] << " ";
351 file << deriv_dim[i] << " ";
352 for (int j=0; j<times.size(); j++) {
353 double rel_time = times[j]/(times_func[j]*d);
354 std::cout << std::setw(w) << rel_time << " ";
355 file << rel_time << " ";
356 }
357 std::cout << std::endl;
358 file << std::endl;
359 }
360 }
361
362 {
363 std::cout.setf(std::ios::scientific);
364 std::cout.precision(p);
365 std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
366 << std::endl;
367 std::cout << std::setw(5) << "deriv" << " ";
368 for (int i=0; i<ExprFuncs::nfunc; i++)
369 std::cout << std::setw(w) << ExprFuncs::add_names[i] << " ";
370 std::cout << std::endl;
371 std::cout << "===== ";
372 for (int i=0; i<ExprFuncs::nfunc; i++) {
373 for (int j=0; j<w; j++)
374 std::cout << '=';
375 std::cout << " ";
376 }
377 std::cout << std::endl;
378
379 // Get function evaluation times
380 double x[ExprFuncs::nx_max];
381 for (int i=0; i<ExprFuncs::nx_max; i++)
382 x[i] = urand.number();
383 int nloop_func = 10000000;
384 Teuchos::Array<double> times_func(ExprFuncs::nfunc);
385 do_times<double,ExprFuncs::add>(x, nloop_func, times_func);
386
387 // Get times for each derivative dimension
388 for (int i=0; i<nderiv; i++) {
389 //int nloop = static_cast<int>(100000.0/(deriv_dim[i]+1));
390 int nloop = static_cast<int>(10000.0);
391 Teuchos::Array<double> times(ExprFuncs::nfunc);
392 do_times_adolc<ExprFuncs::add>(x, seed, deriv_dim[i], nloop, tag, times);
393
394 // Print times
395 int d = deriv_dim[i];
396 if (d == 0)
397 d = 1;
398 std::cout << std::setw(5) << deriv_dim[i] << " ";
399 file << deriv_dim[i] << " ";
400 for (int j=0; j<times.size(); j++) {
401 double rel_time = times[j]/(times_func[j]*d);
402 std::cout << std::setw(w) << rel_time << " ";
403 file << rel_time << " ";
404 }
405 std::cout << std::endl;
406 file << std::endl;
407 }
408 }
409
410 {
411 std::cout.setf(std::ios::scientific);
412 std::cout.precision(p);
413 std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
414 << std::endl;
415 std::cout << std::setw(5) << "deriv" << " ";
416 for (int i=0; i<ExprFuncs::nfunc; i++)
417 std::cout << std::setw(w) << ExprFuncs::nest_names[i] << " ";
418 std::cout << std::endl;
419 std::cout << "===== ";
420 for (int i=0; i<ExprFuncs::nfunc; i++) {
421 for (int j=0; j<w; j++)
422 std::cout << '=';
423 std::cout << " ";
424 }
425 std::cout << std::endl;
426
427 // Get function evaluation times
428 double x[ExprFuncs::nx_max];
429 for (int i=0; i<ExprFuncs::nx_max; i++)
430 x[i] = urand.number();
431 int nloop_func = 10000000;
432 Teuchos::Array<double> times_func(ExprFuncs::nfunc);
433 do_times<double,ExprFuncs::nest>(x, nloop_func, times_func);
434
435 // Get times for each derivative dimension
436 for (int i=0; i<nderiv; i++) {
437 //int nloop = static_cast<int>(10000.0/(deriv_dim[i]+1));
438 int nloop = static_cast<int>(1000.0);
439 Teuchos::Array<double> times(ExprFuncs::nfunc);
440 do_times_adolc<ExprFuncs::nest>(x, seed, deriv_dim[i], nloop, tag, times);
441
442 // Print times
443 int d = deriv_dim[i];
444 if (d == 0)
445 d = 1;
446 std::cout << std::setw(5) << deriv_dim[i] << " ";
447 file << deriv_dim[i] << " ";
448 for (int j=0; j<times.size(); j++) {
449 double rel_time = times[j]/(times_func[j]*d);
450 std::cout << std::setw(w) << rel_time << " ";
451 file << rel_time << " ";
452 }
453 std::cout << std::endl;
454 file << std::endl;
455 }
456 }
457
458 delete [] seed;
459}
460#endif
461
462int main() {
463 print_times< Sacado::Fad::DFad<double> >(
464 "Sacado::Fad::DFad", "fad_expr_depth_dfad.txt");
465 print_times< Sacado::ELRFad::DFad<double> >(
466 "Sacado::ELRFad::DFad", "fad_expr_depth_elr_dfad.txt");
467 print_times< Sacado::CacheFad::DFad<double> >(
468 "Sacado::CacheFad::DFad", "fad_expr_depth_cache_dfad.txt");
469 print_times< Sacado::ELRCacheFad::DFad<double> >(
470 "Sacado::ELRCacheFad::DFad", "fad_expr_depth_elr_cache_dfad.txt");
471 print_times< Sacado::Fad::SimpleFad<double> >(
472 "Sacado::Fad::SimpleFad", "fad_expr_depth_simple_fad.txt");
473#ifdef HAVE_ADOLC
474 print_times_adolc("ADOL-C", "fad_expr_depth_adolc.txt");
475#endif
476
477 return 0;
478}
Sacado::Fad::DFad< double > FadType
A random number generator that generates random numbers uniformly distributed in the interval (a,...
ScalarT number()
Get random number.
void do_times(const T x[], int nloop, Teuchos::Array< double > &times)
double do_time(const T x[], int nloop, const F &f)
void print_times(const std::string &screen_name, const std::string &file_name)
int main()
const double y
const char * p
static const int nfunc
static const char * nest_names[nfunc]
static const char * add_names[nfunc]
static const int nx_max
static const char * mult_names[nfunc]