Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
41 template <typename T, typename F>
42 double 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 
54 template <typename T, template <typename, int> class F>
55 void 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
69 template <typename F>
70 double 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 
99 template <template <typename,int> class F>
100 void 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 
115 template <typename FadType>
116 void 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++) {
152  FadType fx[ExprFuncs::nx_max];
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++) {
208  FadType fx[ExprFuncs::nx_max];
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++) {
264  FadType fx[ExprFuncs::nx_max];
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
295 void 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 
462 int 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 }
ExprFuncs::nx_max
static const int nx_max
Definition: fad_expr_funcs.hpp:71
do_times
void do_times(const T x[], int nloop, Teuchos::Array< double > &times)
Definition: fad_expr_depth.cpp:55
Teuchos_Time.hpp
FadType
Sacado::Fad::DFad< double > FadType
Definition: blas_example.cpp:49
Sacado::Random::number
ScalarT number()
Get random number.
Definition: Sacado_RandomImp.hpp:89
fad_expr_funcs.hpp
T
#define T
Definition: Sacado_rad.hpp:573
Sacado_Random.hpp
Teuchos::Time::start
void start(bool reset=false)
ExprFuncs::mult_names
static const char * mult_names[nfunc]
Definition: fad_expr_funcs.hpp:68
ExprFuncs::add_names
static const char * add_names[nfunc]
Definition: fad_expr_funcs.hpp:70
Teuchos::Time::stop
double stop()
print_times
void print_times(const std::string &screen_name, const std::string &file_name)
Definition: fad_expr_depth.cpp:116
main
int main()
Definition: ad_example.cpp:191
Teuchos_Array.hpp
do_time
double do_time(int nderiv, int nloop)
Definition: fad_expr.cpp:73
Teuchos::Array::size
size_type size() const
Teuchos::Time::totalElapsedTime
double totalElapsedTime(bool readCurrentTime=false) const
ExprFuncs::nfunc
static const int nfunc
Definition: fad_expr_funcs.hpp:67
ExprFuncs::nest_names
static const char * nest_names[nfunc]
Definition: fad_expr_funcs.hpp:69
n
int n
Generated by   doxygen 1.8.5