Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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"
37#include "Teuchos_CommandLineProcessor.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#include "adolc/taping.h"
65#endif
66
67using std::cout;
68using std::endl;
69
70// A simple performance test that computes a Taylor expansion of a simple
71// expression using many variants of Taylor polynomial classes.
72
73template <typename T>
74inline void
75func(const T& x1, const T& x2, T& y) {
76 y = x1*x2 + sin(x1)/x2;
77}
78
79template <typename TaylorType>
80double
81do_time(int degree, int nloop)
82{
83 TaylorType x1, x2, y;
84 Sacado::Random<double> urand(0.0, 1.0);
85
86 x1 = TaylorType(degree, urand.number());
87 x2 = TaylorType(degree, urand.number());
88 y = 0.0;
89 for (int j=0; j<=degree; j++) {
90 x1.fastAccessCoeff(j) = urand.number();
91 x2.fastAccessCoeff(j) = urand.number();
92 }
93
94 Teuchos::Time timer("mult", false);
95 timer.start(true);
96 for (int j=0; j<nloop; j++) {
97 func(x1, x2, y);
98 }
99 timer.stop();
100
101 return timer.totalElapsedTime() / nloop;
102}
103
104#ifdef HAVE_ADOLC
105double
106do_time_adolc(int degree, int nloop)
107{
108 Sacado::Random<double> urand(0.0, 1.0);
109 double **X, **Y;
110
111 X = new double*[2];
112 X[0] = new double[degree+1];
113 X[1] = new double[degree+1];
114 Y = new double*[1];
115 Y[0] = new double[degree+1];
116 for (int j=0; j<=degree; j++) {
117 X[0][j] = urand.number();
118 X[1][j] = urand.number();
119 Y[0][j] = 0.0;
120 }
121
122 trace_on(0);
123 adouble x1, x2, y;
124 x1 <<= X[0][0];
125 x2 <<= X[1][0];
126 func(x1, x2, y);
127 y >>= Y[0][0];
128 trace_off();
129
130 Teuchos::Time timer("mult", false);
131 timer.start(true);
132 for (int j=0; j<nloop; j++) {
133 forward(0,1,2,degree,0,X,Y);
134 }
135 timer.stop();
136
137 delete [] X[1];
138 delete [] X[0];
139 delete [] X;
140
141 delete [] Y[0];
142 delete [] Y;
143
144 return timer.totalElapsedTime() / nloop;
145}
146
147double
148do_time_adolc_retape(int degree, int nloop)
149{
150 Sacado::Random<double> urand(0.0, 1.0);
151 double **X, **Y;
152
153 X = new double*[2];
154 X[0] = new double[degree+1];
155 X[1] = new double[degree+1];
156 Y = new double*[1];
157 Y[0] = new double[degree+1];
158 for (int j=0; j<=degree; j++) {
159 X[0][j] = urand.number();
160 X[1][j] = urand.number();
161 Y[0][j] = 0.0;
162 }
163 adouble x1, x2, y;
164
165 Teuchos::Time timer("mult", false);
166 timer.start(true);
167 for (int j=0; j<nloop; j++) {
168 trace_on(0);
169 x1 <<= X[0][0];
170 x2 <<= X[1][0];
171 func(x1, x2, y);
172 y >>= Y[0][0];
173 trace_off();
174 forward(0,1,2,degree,0,X,Y);
175 }
176 timer.stop();
177
178 delete [] X[1];
179 delete [] X[0];
180 delete [] X;
181
182 delete [] Y[0];
183 delete [] Y;
184
185 return timer.totalElapsedTime() / nloop;
186}
187#endif
188
189int main(int argc, char* argv[]) {
190 int ierr = 0;
191
192 try {
193 double t;
194 int p = 2;
195 int w = p+7;
196
197 // Set up command line options
198 Teuchos::CommandLineProcessor clp;
199 clp.setDocString("This program tests the speed of various forward mode AD implementations for a single multiplication operation");
200 int degree = 10;
201 clp.setOption("degree", &degree, "Polynomial degree");
202 int nloop = 1000000;
203 clp.setOption("nloop", &nloop, "Number of loops");
204
205 // Parse options
206 Teuchos::CommandLineProcessor::EParseCommandLineReturn
207 parseReturn= clp.parse(argc, argv);
208 if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
209 return 1;
210
211 std::cout.setf(std::ios::scientific);
212 std::cout.precision(p);
213 std::cout << "Times (sec) for degree = " << degree
214 << " nloop = " << nloop << ": " << std::endl;
215
216 t = do_time< Sacado::Tay::Taylor<double> >(degree, nloop);
217 std::cout << "Taylor: " << std::setw(w) << t << std::endl;
218
219 t = do_time< Sacado::Tay::CacheTaylor<double> >(degree, nloop);
220 std::cout << "CacheTaylor: " << std::setw(w) << t << std::endl;
221
222#ifdef HAVE_ADOLC
223 t = do_time_adolc(degree, nloop);
224 std::cout << "ADOL-C: " << std::setw(w) << t << std::endl;
225 t = do_time_adolc_retape(degree, nloop);
226 std::cout << "ADOL-C (retape): " << std::setw(w) << t << std::endl;
227#endif
228
229 }
230 catch (std::exception& e) {
231 cout << e.what() << endl;
232 ierr = 1;
233 }
234 catch (const char *s) {
235 cout << s << endl;
236 ierr = 1;
237 }
238 catch (...) {
239 cout << "Caught unknown exception!" << endl;
240 ierr = 1;
241 }
242
243 return ierr;
244}
Sacado::Tay::Taylor< double > TaylorType
sin(expr.val())
int main()
A random number generator that generates random numbers uniformly distributed in the interval (a,...
ScalarT number()
Get random number.
T & fastAccessCoeff(int i)
Returns degree i term without bounds checking.
const double y
const char * p
void func(const T &x1, const T &x2, T &y)
double do_time(int degree, int nloop)