Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
taylor_expr.cpp
Go to the documentation of this file.
1 // $Id$
2 // $Source$
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Sacado Package
7 // Copyright (2006) Sandia Corporation
8 //
9 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // This library is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Lesser General Public License as
14 // published by the Free Software Foundation; either version 2.1 of the
15 // License, or (at your option) any later version.
16 //
17 // This library is distributed in the hope that it will be useful, but
18 // WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 // Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License along with this library; if not, write to the Free Software
24 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
25 // USA
26 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
27 // (etphipp@sandia.gov).
28 //
29 // ***********************************************************************
30 // @HEADER
31 
32 #include "Sacado_Random.hpp"
33 #include "Sacado_No_Kokkos.hpp"
35 
36 #include "Teuchos_Time.hpp"
38 
39 // ADOL-C includes
40 #ifdef HAVE_ADOLC
41 #ifdef PACKAGE
42 #undef PACKAGE
43 #endif
44 #ifdef PACKAGE_NAME
45 #undef PACKAGE_NAME
46 #endif
47 #ifdef PACKAGE_BUGREPORT
48 #undef PACKAGE_BUGREPORT
49 #endif
50 #ifdef PACKAGE_STRING
51 #undef PACKAGE_STRING
52 #endif
53 #ifdef PACKAGE_TARNAME
54 #undef PACKAGE_TARNAME
55 #endif
56 #ifdef PACKAGE_VERSION
57 #undef PACKAGE_VERSION
58 #endif
59 #ifdef VERSION
60 #undef VERSION
61 #endif
62 #include "adolc/adouble.h"
63 #include "adolc/interfaces.h"
64 #endif
65 
66 using std::cout;
67 using std::endl;
68 
69 // A simple performance test that computes a Taylor expansion of a simple
70 // expression using many variants of Taylor polynomial classes.
71 
72 template <typename T>
73 inline void
74 func(const T& x1, const T& x2, T& y) {
75  y = x1*x2 + sin(x1)/x2;
76 }
77 
78 template <typename TaylorType>
79 double
80 do_time(int degree, int nloop)
81 {
82  TaylorType x1, x2, y;
83  Sacado::Random<double> urand(0.0, 1.0);
84 
85  x1 = TaylorType(degree, urand.number());
86  x2 = TaylorType(degree, urand.number());
87  y = 0.0;
88  for (int j=0; j<=degree; j++) {
89  x1.fastAccessCoeff(j) = urand.number();
90  x2.fastAccessCoeff(j) = urand.number();
91  }
92 
93  Teuchos::Time timer("mult", false);
94  timer.start(true);
95  for (int j=0; j<nloop; j++) {
96  func(x1, x2, y);
97  }
98  timer.stop();
99 
100  return timer.totalElapsedTime() / nloop;
101 }
102 
103 #ifdef HAVE_ADOLC
104 double
105 do_time_adolc(int degree, int nloop)
106 {
107  Sacado::Random<double> urand(0.0, 1.0);
108  double **X, **Y;
109 
110  X = new double*[2];
111  X[0] = new double[degree+1];
112  X[1] = new double[degree+1];
113  Y = new double*[1];
114  Y[0] = new double[degree+1];
115  for (int j=0; j<=degree; j++) {
116  X[0][j] = urand.number();
117  X[1][j] = urand.number();
118  Y[0][j] = 0.0;
119  }
120 
121  trace_on(0);
122  adouble x1, x2, y;
123  x1 <<= X[0][0];
124  x2 <<= X[1][0];
125  func(x1, x2, y);
126  y >>= Y[0][0];
127  trace_off();
128 
129  Teuchos::Time timer("mult", false);
130  timer.start(true);
131  for (int j=0; j<nloop; j++) {
132  forward(0,1,2,degree,0,X,Y);
133  }
134  timer.stop();
135 
136  delete [] X[1];
137  delete [] X[0];
138  delete [] X;
139 
140  delete [] Y[0];
141  delete [] Y;
142 
143  return timer.totalElapsedTime() / nloop;
144 }
145 
146 double
147 do_time_adolc_retape(int degree, int nloop)
148 {
149  Sacado::Random<double> urand(0.0, 1.0);
150  double **X, **Y;
151 
152  X = new double*[2];
153  X[0] = new double[degree+1];
154  X[1] = new double[degree+1];
155  Y = new double*[1];
156  Y[0] = new double[degree+1];
157  for (int j=0; j<=degree; j++) {
158  X[0][j] = urand.number();
159  X[1][j] = urand.number();
160  Y[0][j] = 0.0;
161  }
162  adouble x1, x2, y;
163 
164  Teuchos::Time timer("mult", false);
165  timer.start(true);
166  for (int j=0; j<nloop; j++) {
167  trace_on(0);
168  x1 <<= X[0][0];
169  x2 <<= X[1][0];
170  func(x1, x2, y);
171  y >>= Y[0][0];
172  trace_off();
173  forward(0,1,2,degree,0,X,Y);
174  }
175  timer.stop();
176 
177  delete [] X[1];
178  delete [] X[0];
179  delete [] X;
180 
181  delete [] Y[0];
182  delete [] Y;
183 
184  return timer.totalElapsedTime() / nloop;
185 }
186 #endif
187 
188 int main(int argc, char* argv[]) {
189  int ierr = 0;
190 
191  try {
192  double t;
193  int p = 2;
194  int w = p+7;
195 
196  // Set up command line options
198  clp.setDocString("This program tests the speed of various forward mode AD implementations for a single multiplication operation");
199  int degree = 10;
200  clp.setOption("degree", &degree, "Polynomial degree");
201  int nloop = 1000000;
202  clp.setOption("nloop", &nloop, "Number of loops");
203 
204  // Parse options
206  parseReturn= clp.parse(argc, argv);
208  return 1;
209 
210  std::cout.setf(std::ios::scientific);
211  std::cout.precision(p);
212  std::cout << "Times (sec) for degree = " << degree
213  << " nloop = " << nloop << ": " << std::endl;
214 
215  t = do_time< Sacado::Tay::Taylor<double> >(degree, nloop);
216  std::cout << "Taylor: " << std::setw(w) << t << std::endl;
217 
218  t = do_time< Sacado::Tay::CacheTaylor<double> >(degree, nloop);
219  std::cout << "CacheTaylor: " << std::setw(w) << t << std::endl;
220 
221 #ifdef HAVE_ADOLC
222  t = do_time_adolc(degree, nloop);
223  std::cout << "ADOL-C: " << std::setw(w) << t << std::endl;
224  t = do_time_adolc_retape(degree, nloop);
225  std::cout << "ADOL-C (retape): " << std::setw(w) << t << std::endl;
226 #endif
227 
228  }
229  catch (std::exception& e) {
230  cout << e.what() << endl;
231  ierr = 1;
232  }
233  catch (const char *s) {
234  cout << s << endl;
235  ierr = 1;
236  }
237  catch (...) {
238  cout << "Caught unknown exception!" << endl;
239  ierr = 1;
240  }
241 
242  return ierr;
243 }
ScalarT number()
Get random number.
Sacado::Tay::Taylor< double > TaylorType
#define T
Definition: Sacado_rad.hpp:573
void start(bool reset=false)
T & fastAccessCoeff(int i)
Returns degree i term without bounds checking.
double stop()
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
int main()
Definition: ad_example.cpp:191
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
sin(expr.val())
double do_time(int nderiv, int nloop)
Definition: fad_expr.cpp:73
void setDocString(const char doc_string[])
double totalElapsedTime(bool readCurrentTime=false) const
const T func(int n, T *x)
Definition: ad_example.cpp:49