Sacado Package Browser (Single Doxygen Collection)  Version of the Day
fad_lj_grad.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"
34 #include "Sacado_CacheFad_DFad.hpp"
35 
36 #include "Fad/fad.h"
37 #include "TinyFadET/tfad.h"
38 
39 #include "Teuchos_Time.hpp"
40 #include "Teuchos_CommandLineProcessor.hpp"
41 
42 // A simple performance test that computes the derivative of a simple
43 // expression using many variants of Fad.
44 
45 template <>
47 
48 void FAD::error(const char *msg) {
49  std::cout << msg << std::endl;
50 }
51 
52 namespace {
53  double xi[3], xj[3], pa[4], f[3], delr[3];
54 }
55 
56 template <typename T>
57 inline T
58 vec3_distsq(const T xi[], const double xj[]) {
59  T delr0 = xi[0]-xj[0];
60  T delr1 = xi[1]-xj[1];
61  T delr2 = xi[2]-xj[2];
62  return delr0*delr0 + delr1*delr1 + delr2*delr2;
63 }
64 
65 template <typename T>
66 inline T
67 vec3_distsq(const T xi[], const double xj[], T delr[]) {
68  delr[0] = xi[0]-xj[0];
69  delr[1] = xi[1]-xj[1];
70  delr[2] = xi[2]-xj[2];
71  return delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
72 }
73 
74 template <typename T>
75 inline void
76 lj(const T xi[], const double xj[], T& energy) {
77  T delr2 = vec3_distsq(xi,xj);
78  T delr_2 = 1.0/delr2;
79  T delr_6 = delr_2*delr_2*delr_2;
80  energy = (pa[1]*delr_6 - pa[2])*delr_6 - pa[3];
81 }
82 
83 inline void
84 lj_and_grad(const double xi[], const double xj[], double& energy,
85  double f[]) {
86  double delr2 = vec3_distsq(xi,xj,delr);
87  double delr_2 = 1.0/delr2;
88  double delr_6 = delr_2*delr_2*delr_2;
89  energy = (pa[1]*delr_6 - pa[2])*delr_6 - pa[3];
90  double tmp = (-12.0*pa[1]*delr_6 - 6.0*pa[2])*delr_6*delr_2;
91  f[0] = delr[0]*tmp;
92  f[1] = delr[1]*tmp;
93  f[2] = delr[2]*tmp;
94 }
95 
96 template <typename FadType>
97 double
98 do_time(int nloop)
99 {
100  Teuchos::Time timer("lj", false);
101  FadType xi_fad[3], energy;
102 
103  for (int i=0; i<3; i++) {
104  xi_fad[i] = FadType(3, i, xi[i]);
105  }
106 
107  timer.start(true);
108  for (int j=0; j<nloop; j++) {
109 
110  lj(xi_fad, xj, energy);
111 
112  for (int i=0; i<3; i++)
113  f[i] += -energy.fastAccessDx(i);
114  }
115  timer.stop();
116 
117  return timer.totalElapsedTime() / nloop;
118 }
119 
120 double
122 {
123  Teuchos::Time timer("lj", false);
124  double energy, ff[3];
125 
126  timer.start(true);
127  for (int j=0; j<nloop; j++) {
128 
129  lj_and_grad(xi, xj, energy, ff);
130 
131  for (int i=0; i<3; i++)
132  f[i] += -ff[i];
133 
134  }
135  timer.stop();
136 
137  return timer.totalElapsedTime() / nloop;
138 }
139 
140 int main(int argc, char* argv[]) {
141  int ierr = 0;
142 
143  try {
144  double t, ta;
145  int p = 2;
146  int w = p+7;
147 
148  // Set up command line options
149  Teuchos::CommandLineProcessor clp;
150  clp.setDocString("This program tests the speed of various forward mode AD implementations for a single multiplication operation");
151  int nloop = 1000000;
152  clp.setOption("nloop", &nloop, "Number of loops");
153 
154  // Parse options
155  Teuchos::CommandLineProcessor::EParseCommandLineReturn
156  parseReturn= clp.parse(argc, argv);
157  if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
158  return 1;
159 
160  // Memory pool & manager
162  Sacado::Fad::MemPool* pool = poolManager.getMemoryPool(3);
164 
165  std::cout.setf(std::ios::scientific);
166  std::cout.precision(p);
167  std::cout << "Times (sec) nloop = " << nloop << ": " << std::endl;
168 
169  Sacado::Random<double> urand(0.0, 1.0);
170  for (int i=0; i<3; i++) {
171  xi[i] = urand.number();
172  xj[i] = urand.number();
173  pa[i] = urand.number();
174  }
175  pa[3] = urand.number();
176 
177  ta = do_time_analytic(nloop);
178  std::cout << "Analytic: " << std::setw(w) << ta << std::endl;
179 
180  t = do_time< FAD::TFad<3,double> >(nloop);
181  std::cout << "TFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
182 
183  t = do_time< FAD::Fad<double> >(nloop);
184  std::cout << "Fad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
185 
186  t = do_time< Sacado::Fad::SFad<double,3> >(nloop);
187  std::cout << "SFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
188 
189  t = do_time< Sacado::Fad::SLFad<double,3> >(nloop);
190  std::cout << "SLFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
191 
192  t = do_time< Sacado::Fad::DFad<double> >(nloop);
193  std::cout << "DFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
194 
195  t = do_time< Sacado::Fad::DMFad<double> >(nloop);
196  std::cout << "DMFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
197 
198  t = do_time< Sacado::ELRFad::SFad<double,3> >(nloop);
199  std::cout << "ELRSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
200 
201  t = do_time< Sacado::ELRFad::SLFad<double,3> >(nloop);
202  std::cout << "ELRSLFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
203 
204  t = do_time< Sacado::ELRFad::DFad<double> >(nloop);
205  std::cout << "ELRDFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
206 
207  t = do_time< Sacado::CacheFad::DFad<double> >(nloop);
208  std::cout << "CacheFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
209 
210  t = do_time< Sacado::Fad::DVFad<double> >(nloop);
211  std::cout << "DVFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << std::endl;
212 
213  }
214  catch (std::exception& e) {
215  std::cout << e.what() << std::endl;
216  ierr = 1;
217  }
218  catch (const char *s) {
219  std::cout << s << std::endl;
220  ierr = 1;
221  }
222  catch (...) {
223  std::cout << "Caught unknown exception!" << std::endl;
224  ierr = 1;
225  }
226 
227  return ierr;
228 }
MemPool * getMemoryPool(unsigned int dim)
Get memory pool for supplied dimension dim.
Sacado::Fad::DFad< double > FadType
double do_time_analytic(int nloop)
ScalarT number()
Get random number.
double do_time(int nloop)
Definition: fad_lj_grad.cpp:98
#define T
Definition: Sacado_rad.hpp:573
int main(int argc, char *argv[])
Derivative array storage class using dynamic memory allocation.
void lj_and_grad(const double xi[], const double xj[], double &energy, double f[])
Definition: fad_lj_grad.cpp:84
Forward-mode AD class using dynamic memory allocation and expression templates.
void lj(const T xi[], const double xj[], T &energy)
Definition: fad_lj_grad.cpp:76
T vec3_distsq(const T xi[], const double xj[])
Definition: fad_lj_grad.cpp:58