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 // Sacado Package
4 //
5 // Copyright 2006 NTESS and the Sacado contributors.
6 // SPDX-License-Identifier: LGPL-2.1-or-later
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Sacado_Random.hpp"
11 
12 #include "Teuchos_Time.hpp"
13 #include "Teuchos_Array.hpp"
14 #include <fstream>
15 
16 #include "fad_expr_funcs.hpp"
17 
18 // A simple performance test that computes the derivative of expressions of
19 // various depths.
20 
21 template <typename T, typename F>
22 double do_time(const T x[], int nloop, const F& f)
23 {
24  T y = 0.0;
25  Teuchos::Time timer("F", false);
26 
27  timer.start(true);
28  for (int j=0; j<nloop; j++)
29  f(x, y);
30  timer.stop();
31  return timer.totalElapsedTime() / nloop;
32 }
33 
34 template <typename T, template <typename, int> class F>
35 void do_times(const T x[], int nloop, Teuchos::Array<double>& times)
36 {
37  int i = 0;
38  times[i++] = do_time(x, nloop, F<T,1>());
39  times[i++] = do_time(x, nloop, F<T,2>());
40  times[i++] = do_time(x, nloop, F<T,3>());
41  times[i++] = do_time(x, nloop, F<T,4>());
42  times[i++] = do_time(x, nloop, F<T,5>());
43  times[i++] = do_time(x, nloop, F<T,10>());
44  times[i++] = do_time(x, nloop, F<T,15>());
45  times[i++] = do_time(x, nloop, F<T,20>());
46 }
47 
48 #ifdef HAVE_ADOLC
49 template <typename F>
50 double do_time_adolc(double *x, double **seed, int d, int nloop,
51  int tag, const F& f)
52 {
53  Teuchos::Time timer("F", false);
54  int n = F::n;
55  trace_on(tag);
56  adouble *x_ad = new adouble[n];
57  for (int i=0; i<n; i++)
58  x_ad[i] <<= x[i];
59  adouble y_ad = 0.0;
60  f(x_ad, y_ad);
61  double y;
62  y_ad >>= y;
63  delete [] x_ad;
64  trace_off();
65 
66  double **jac = new double*[1];
67  jac[0] = new double[d];
68  timer.start(true);
69  for (int j=0; j<nloop; j++)
70  forward(tag, 1, n, d, x, seed, &y, jac);
71  timer.stop();
72 
73  delete [] jac[0];
74  delete [] jac;
75 
76  return timer.totalElapsedTime() / nloop;
77 }
78 
79 template <template <typename,int> class F>
80 void do_times_adolc(double *x, double **seed, int d, int nloop,
81  int& tag, Teuchos::Array<double>& times)
82 {
83  int i = 0;
84  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,1>());
85  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,2>());
86  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,3>());
87  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,4>());
88  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,5>());
89  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,10>());
90  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,15>());
91  times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,20>());
92 }
93 #endif
94 
95 template <typename FadType>
96 void print_times(const std::string& screen_name, const std::string& file_name)
97 {
98  const int nderiv = 10;
99  int deriv_dim[nderiv] = { 0, 1, 3, 5, 10, 15, 20, 30, 40, 50 };
100  Sacado::Random<double> urand(0.0, 1.0);
101  int p = 1;
102  int w = p+7;
103  std::ofstream file(file_name.c_str(), std::ios::out);
104 
105  {
106  std::cout.setf(std::ios::scientific);
107  std::cout.precision(p);
108  std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
109  << std::endl;
110  std::cout << std::setw(5) << "deriv" << " ";
111  for (int i=0; i<ExprFuncs::nfunc; i++)
112  std::cout << std::setw(w) << ExprFuncs::mult_names[i] << " ";
113  std::cout << std::endl;
114  std::cout << "===== ";
115  for (int i=0; i<ExprFuncs::nfunc; i++) {
116  for (int j=0; j<w; j++)
117  std::cout << '=';
118  std::cout << " ";
119  }
120  std::cout << std::endl;
121 
122  // Get function evaluation times
123  double x[ExprFuncs::nx_max];
124  for (int i=0; i<ExprFuncs::nx_max; i++)
125  x[i] = urand.number();
126  int nloop_func = 10000000;
127  Teuchos::Array<double> times_func(ExprFuncs::nfunc);
128  do_times<double,ExprFuncs::mult>(x, nloop_func, times_func);
129 
130  // Get times for each derivative dimension
131  for (int i=0; i<nderiv; i++) {
133  for (int k=0; k<ExprFuncs::nx_max; k++) {
134  fx[k] = FadType(deriv_dim[i], urand.number());
135  for (int j=0; j<deriv_dim[i]; j++) {
136  fx[k].fastAccessDx(j) = urand.number();
137  }
138  }
139 
140  //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
141  int nloop = static_cast<int>(100000.0);
142  Teuchos::Array<double> times(ExprFuncs::nfunc);
143  do_times<FadType,ExprFuncs::mult>(fx, nloop, times);
144 
145  // Print times
146  int d = deriv_dim[i];
147  if (d == 0)
148  d = 1;
149  std::cout << std::setw(5) << deriv_dim[i] << " ";
150  file << deriv_dim[i] << " ";
151  for (int j=0; j<times.size(); j++) {
152  double rel_time = times[j]/(times_func[j]*d);
153  std::cout << std::setw(w) << rel_time << " ";
154  file << rel_time << " ";
155  }
156  std::cout << std::endl;
157  file << std::endl;
158  }
159  }
160 
161  {
162  std::cout.setf(std::ios::scientific);
163  std::cout.precision(p);
164  std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
165  << std::endl;
166  std::cout << std::setw(5) << "deriv" << " ";
167  for (int i=0; i<ExprFuncs::nfunc; i++)
168  std::cout << std::setw(w) << ExprFuncs::add_names[i] << " ";
169  std::cout << std::endl;
170  std::cout << "===== ";
171  for (int i=0; i<ExprFuncs::nfunc; i++) {
172  for (int j=0; j<w; j++)
173  std::cout << '=';
174  std::cout << " ";
175  }
176  std::cout << std::endl;
177 
178  // Get function evaluation times
179  double x[ExprFuncs::nx_max];
180  for (int i=0; i<ExprFuncs::nx_max; i++)
181  x[i] = urand.number();
182  int nloop_func = 10000000;
183  Teuchos::Array<double> times_func(ExprFuncs::nfunc);
184  do_times<double,ExprFuncs::add>(x, nloop_func, times_func);
185 
186  // Get times for each derivative dimension
187  for (int i=0; i<nderiv; i++) {
189  for (int k=0; k<ExprFuncs::nx_max; k++) {
190  fx[k] = FadType(deriv_dim[i], urand.number());
191  for (int j=0; j<deriv_dim[i]; j++) {
192  fx[k].fastAccessDx(j) = urand.number();
193  }
194  }
195 
196  //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
197  int nloop = static_cast<int>(100000.0);
198  Teuchos::Array<double> times(ExprFuncs::nfunc);
199  do_times<FadType,ExprFuncs::add>(fx, nloop, times);
200 
201  // Print times
202  int d = deriv_dim[i];
203  if (d == 0)
204  d = 1;
205  std::cout << std::setw(5) << deriv_dim[i] << " ";
206  file << deriv_dim[i] << " ";
207  for (int j=0; j<times.size(); j++) {
208  double rel_time = times[j]/(times_func[j]*d);
209  std::cout << std::setw(w) << rel_time << " ";
210  file << rel_time << " ";
211  }
212  std::cout << std::endl;
213  file << std::endl;
214  }
215  }
216 
217  {
218  std::cout.setf(std::ios::scientific);
219  std::cout.precision(p);
220  std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
221  << std::endl;
222  std::cout << std::setw(5) << "deriv" << " ";
223  for (int i=0; i<ExprFuncs::nfunc; i++)
224  std::cout << std::setw(w) << ExprFuncs::nest_names[i] << " ";
225  std::cout << std::endl;
226  std::cout << "===== ";
227  for (int i=0; i<ExprFuncs::nfunc; i++) {
228  for (int j=0; j<w; j++)
229  std::cout << '=';
230  std::cout << " ";
231  }
232  std::cout << std::endl;
233 
234  // Get function evaluation times
235  double x[ExprFuncs::nx_max];
236  for (int i=0; i<ExprFuncs::nx_max; i++)
237  x[i] = urand.number();
238  int nloop_func = 10000000;
239  Teuchos::Array<double> times_func(ExprFuncs::nfunc);
240  do_times<double,ExprFuncs::nest>(x, nloop_func, times_func);
241 
242  // Get times for each derivative dimension
243  for (int i=0; i<nderiv; i++) {
245  for (int k=0; k<ExprFuncs::nx_max; k++) {
246  fx[k] = FadType(deriv_dim[i], urand.number());
247  for (int j=0; j<deriv_dim[i]; j++) {
248  fx[k].fastAccessDx(j) = urand.number();
249  }
250  }
251 
252  //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
253  int nloop = static_cast<int>(100000.0);
254  Teuchos::Array<double> times(ExprFuncs::nfunc);
255  do_times<FadType,ExprFuncs::nest>(fx, nloop, times);
256 
257  // Print times
258  int d = deriv_dim[i];
259  if (d == 0)
260  d = 1;
261  std::cout << std::setw(5) << deriv_dim[i] << " ";
262  file << deriv_dim[i] << " ";
263  for (int j=0; j<times.size(); j++) {
264  double rel_time = times[j]/(times_func[j]*d);
265  std::cout << std::setw(w) << rel_time << " ";
266  file << rel_time << " ";
267  }
268  std::cout << std::endl;
269  file << std::endl;
270  }
271  }
272 }
273 
274 #ifdef HAVE_ADOLC
275 void print_times_adolc(const std::string& screen_name,
276  const std::string& file_name)
277 {
278  const int nderiv = 10;
279  int deriv_dim[nderiv] = { 0, 1, 3, 5, 10, 15, 20, 30, 40, 50 };
280  const int deriv_max = 50;
281  Sacado::Random<double> urand(0.0, 1.0);
282  int p = 1;
283  int w = p+7;
284  std::ofstream file(file_name.c_str(), std::ios::out);
285  int tag = 0;
286 
287  double **seed = new double*[ExprFuncs::nx_max];
288  for (int i=0; i<ExprFuncs::nx_max; i++) {
289  seed[i] = new double[deriv_max];
290  for (int j=0; j<deriv_max; j++)
291  seed[i][j] = urand.number();
292  }
293 
294  {
295  std::cout.setf(std::ios::scientific);
296  std::cout.precision(p);
297  std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
298  << std::endl;
299  std::cout << std::setw(5) << "deriv" << " ";
300  for (int i=0; i<ExprFuncs::nfunc; i++)
301  std::cout << std::setw(w) << ExprFuncs::mult_names[i] << " ";
302  std::cout << std::endl;
303  std::cout << "===== ";
304  for (int i=0; i<ExprFuncs::nfunc; i++) {
305  for (int j=0; j<w; j++)
306  std::cout << '=';
307  std::cout << " ";
308  }
309  std::cout << std::endl;
310 
311  // Get function evaluation times
312  double x[ExprFuncs::nx_max];
313  for (int i=0; i<ExprFuncs::nx_max; i++)
314  x[i] = urand.number();
315  int nloop_func = 10000000;
316  Teuchos::Array<double> times_func(ExprFuncs::nfunc);
317  do_times<double,ExprFuncs::mult>(x, nloop_func, times_func);
318 
319  // Get times for each derivative dimension
320  for (int i=0; i<nderiv; i++) {
321  //int nloop = static_cast<int>(100000.0/(deriv_dim[i]+1));
322  int nloop = static_cast<int>(10000.0);
323  Teuchos::Array<double> times(ExprFuncs::nfunc);
324  do_times_adolc<ExprFuncs::mult>(x, seed, deriv_dim[i], nloop, tag, times);
325 
326  // Print times
327  int d = deriv_dim[i];
328  if (d == 0)
329  d = 1;
330  std::cout << std::setw(5) << deriv_dim[i] << " ";
331  file << deriv_dim[i] << " ";
332  for (int j=0; j<times.size(); j++) {
333  double rel_time = times[j]/(times_func[j]*d);
334  std::cout << std::setw(w) << rel_time << " ";
335  file << rel_time << " ";
336  }
337  std::cout << std::endl;
338  file << std::endl;
339  }
340  }
341 
342  {
343  std::cout.setf(std::ios::scientific);
344  std::cout.precision(p);
345  std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
346  << std::endl;
347  std::cout << std::setw(5) << "deriv" << " ";
348  for (int i=0; i<ExprFuncs::nfunc; i++)
349  std::cout << std::setw(w) << ExprFuncs::add_names[i] << " ";
350  std::cout << std::endl;
351  std::cout << "===== ";
352  for (int i=0; i<ExprFuncs::nfunc; i++) {
353  for (int j=0; j<w; j++)
354  std::cout << '=';
355  std::cout << " ";
356  }
357  std::cout << std::endl;
358 
359  // Get function evaluation times
360  double x[ExprFuncs::nx_max];
361  for (int i=0; i<ExprFuncs::nx_max; i++)
362  x[i] = urand.number();
363  int nloop_func = 10000000;
364  Teuchos::Array<double> times_func(ExprFuncs::nfunc);
365  do_times<double,ExprFuncs::add>(x, nloop_func, times_func);
366 
367  // Get times for each derivative dimension
368  for (int i=0; i<nderiv; i++) {
369  //int nloop = static_cast<int>(100000.0/(deriv_dim[i]+1));
370  int nloop = static_cast<int>(10000.0);
371  Teuchos::Array<double> times(ExprFuncs::nfunc);
372  do_times_adolc<ExprFuncs::add>(x, seed, deriv_dim[i], nloop, tag, times);
373 
374  // Print times
375  int d = deriv_dim[i];
376  if (d == 0)
377  d = 1;
378  std::cout << std::setw(5) << deriv_dim[i] << " ";
379  file << deriv_dim[i] << " ";
380  for (int j=0; j<times.size(); j++) {
381  double rel_time = times[j]/(times_func[j]*d);
382  std::cout << std::setw(w) << rel_time << " ";
383  file << rel_time << " ";
384  }
385  std::cout << std::endl;
386  file << std::endl;
387  }
388  }
389 
390  {
391  std::cout.setf(std::ios::scientific);
392  std::cout.precision(p);
393  std::cout << screen_name << " Relative times (time/(func_time*nderiv)): "
394  << std::endl;
395  std::cout << std::setw(5) << "deriv" << " ";
396  for (int i=0; i<ExprFuncs::nfunc; i++)
397  std::cout << std::setw(w) << ExprFuncs::nest_names[i] << " ";
398  std::cout << std::endl;
399  std::cout << "===== ";
400  for (int i=0; i<ExprFuncs::nfunc; i++) {
401  for (int j=0; j<w; j++)
402  std::cout << '=';
403  std::cout << " ";
404  }
405  std::cout << std::endl;
406 
407  // Get function evaluation times
408  double x[ExprFuncs::nx_max];
409  for (int i=0; i<ExprFuncs::nx_max; i++)
410  x[i] = urand.number();
411  int nloop_func = 10000000;
412  Teuchos::Array<double> times_func(ExprFuncs::nfunc);
413  do_times<double,ExprFuncs::nest>(x, nloop_func, times_func);
414 
415  // Get times for each derivative dimension
416  for (int i=0; i<nderiv; i++) {
417  //int nloop = static_cast<int>(10000.0/(deriv_dim[i]+1));
418  int nloop = static_cast<int>(1000.0);
419  Teuchos::Array<double> times(ExprFuncs::nfunc);
420  do_times_adolc<ExprFuncs::nest>(x, seed, deriv_dim[i], nloop, tag, times);
421 
422  // Print times
423  int d = deriv_dim[i];
424  if (d == 0)
425  d = 1;
426  std::cout << std::setw(5) << deriv_dim[i] << " ";
427  file << deriv_dim[i] << " ";
428  for (int j=0; j<times.size(); j++) {
429  double rel_time = times[j]/(times_func[j]*d);
430  std::cout << std::setw(w) << rel_time << " ";
431  file << rel_time << " ";
432  }
433  std::cout << std::endl;
434  file << std::endl;
435  }
436  }
437 
438  delete [] seed;
439 }
440 #endif
441 
442 int main() {
443  print_times< Sacado::Fad::DFad<double> >(
444  "Sacado::Fad::DFad", "fad_expr_depth_dfad.txt");
445  print_times< Sacado::ELRFad::DFad<double> >(
446  "Sacado::ELRFad::DFad", "fad_expr_depth_elr_dfad.txt");
447  print_times< Sacado::CacheFad::DFad<double> >(
448  "Sacado::CacheFad::DFad", "fad_expr_depth_cache_dfad.txt");
449  print_times< Sacado::ELRCacheFad::DFad<double> >(
450  "Sacado::ELRCacheFad::DFad", "fad_expr_depth_elr_cache_dfad.txt");
451  print_times< Sacado::Fad::SimpleFad<double> >(
452  "Sacado::Fad::SimpleFad", "fad_expr_depth_simple_fad.txt");
453 #ifdef HAVE_ADOLC
454  print_times_adolc("ADOL-C", "fad_expr_depth_adolc.txt");
455 #endif
456 
457  return 0;
458 }
const char * p
static const int nx_max
void do_times(const T x[], int nloop, Teuchos::Array< double > &times)
Sacado::Fad::DFad< double > FadType
ScalarT number()
Get random number.
#define T
Definition: Sacado_rad.hpp:553
void start(bool reset=false)
static const char * mult_names[nfunc]
static const char * add_names[nfunc]
double stop()
void print_times(const std::string &screen_name, const std::string &file_name)
int main()
Definition: ad_example.cpp:171
double do_time(int nderiv, int nloop)
Definition: fad_expr.cpp:48
size_type size() const
double totalElapsedTime(bool readCurrentTime=false) const
static const int nfunc
static const char * nest_names[nfunc]
int n
const double y