Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fad_fe_jac_fill.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 "fe_jac_fill_funcs.hpp"
33 
34 #include "Fad/fad.h"
35 #include "TinyFadET/tfad.h"
36 
37 // A performance test that computes a finite-element-like Jacobian using
38 // several Fad variants
39 
40 template <>
42 
43 void FAD::error(const char *msg) {
44  std::cout << msg << std::endl;
45 }
46 
47 int main(int argc, char* argv[]) {
48  int ierr = 0;
49 
50  try {
51  double t, ta, tr;
52  int p = 2;
53  int w = p+7;
54 
55  // Maximum number of derivative components for SLFad
56  const int slfad_max = 130;
57 
58  // Set up command line options
60  clp.setDocString("This program tests the speed of various forward mode AD implementations for a finite-element-like Jacobian fill");
61  int num_nodes = 100000;
62  int num_eqns = 2;
63  int rt = 0;
64  clp.setOption("n", &num_nodes, "Number of nodes");
65  clp.setOption("p", &num_eqns, "Number of equations");
66  clp.setOption("rt", &rt, "Include ADOL-C retaping test");
67 
68  // Parse options
70  parseReturn= clp.parse(argc, argv);
72  return 1;
73 
74  double mesh_spacing = 1.0 / static_cast<double>(num_nodes - 1);
75 
76  // Memory pool & manager
77  Sacado::Fad::MemPoolManager<double> poolManager(num_nodes*num_eqns);
78  Sacado::Fad::MemPool* pool = poolManager.getMemoryPool(num_nodes*num_eqns);
80 
81  std::cout.setf(std::ios::scientific);
82  std::cout.precision(p);
83  std::cout << "num_nodes = " << num_nodes
84  << ", num_eqns = " << num_eqns << ": " << std::endl
85  << " " << " Time " << "\t"<< "Time/Analytic" << "\t"
86  << "Time/(2*p*Residual)" << std::endl;
87 
88  ta = 1.0;
89  tr = 1.0;
90 
91  tr = residual_fill(num_nodes, num_eqns, mesh_spacing);
92 
93  ta = analytic_jac_fill(num_nodes, num_eqns, mesh_spacing);
94  std::cout << "Analytic: " << std::setw(w) << ta << "\t" << std::setw(w) << ta/ta << "\t" << std::setw(w) << ta/(2.0*num_eqns*tr) << std::endl;
95 
96 #ifdef HAVE_ADOLC
97 #ifndef ADOLC_TAPELESS
98  t = adolc_jac_fill(num_nodes, num_eqns, mesh_spacing);
99  std::cout << "ADOL-C: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
100 
101  if (rt != 0) {
102  t = adolc_retape_jac_fill(num_nodes, num_eqns, mesh_spacing);
103  std::cout << "ADOL-C(rt): " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
104  }
105 
106 #else
107  t = adolc_tapeless_jac_fill(num_nodes, num_eqns, mesh_spacing);
108  std::cout << "ADOL-C(tl): " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
109 #endif
110 #endif
111 
112 #ifdef HAVE_ADIC
113  t = adic_jac_fill(num_nodes, num_eqns, mesh_spacing);
114  std::cout << "ADIC: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
115 #endif
116 
117  if (num_eqns*2 == 4) {
118  t = fad_jac_fill< FAD::TFad<16,double> >(num_nodes, num_eqns, mesh_spacing);
119  std::cout << "TFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
120  }
121  else if (num_eqns*2 == 16) {
122  t = fad_jac_fill< FAD::TFad<16,double> >(num_nodes, num_eqns, mesh_spacing);
123  std::cout << "TFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
124  }
125  else if (num_eqns*2 == 32) {
126  t = fad_jac_fill< FAD::TFad<32,double> >(num_nodes, num_eqns, mesh_spacing);
127  std::cout << "TFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
128  }
129  else if (num_eqns*2 == 64) {
130  t = fad_jac_fill< FAD::TFad<64,double> >(num_nodes, num_eqns, mesh_spacing);
131  std::cout << "TFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
132  }
133 
134  t = fad_jac_fill< FAD::Fad<double> >(num_nodes, num_eqns, mesh_spacing);
135  std::cout << "Fad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
136 
137  if (num_eqns*2 == 4) {
138  t = fad_jac_fill< Sacado::Fad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
139  std::cout << "SFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
140  }
141  else if (num_eqns*2 == 16) {
142  t = fad_jac_fill< Sacado::Fad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
143  std::cout << "SFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
144  }
145  else if (num_eqns*2 == 32) {
146  t = fad_jac_fill< Sacado::Fad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
147  std::cout << "SFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
148  }
149  else if (num_eqns*2 == 64) {
150  t = fad_jac_fill< Sacado::Fad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
151  std::cout << "SFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
152  }
153 
154  if (num_eqns*2 < slfad_max) {
155  t = fad_jac_fill< Sacado::Fad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
156  std::cout << "SLFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
157  }
158 
159  t = fad_jac_fill< Sacado::Fad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
160  std::cout << "DFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
161 
162  t = fad_jac_fill< Sacado::Fad::SimpleFad<double> >(num_nodes, num_eqns, mesh_spacing);
163  std::cout << "SimpleFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
164 
165  t = fad_jac_fill< Sacado::Fad::DMFad<double> >(num_nodes, num_eqns, mesh_spacing);
166  std::cout << "DMFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
167 
168  if (num_eqns*2 == 4) {
169  t = fad_jac_fill< Sacado::ELRFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
170  std::cout << "ELRSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
171  }
172  else if (num_eqns*2 == 16) {
173  t = fad_jac_fill< Sacado::ELRFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
174  std::cout << "ELRSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
175  }
176  else if (num_eqns*2 == 32) {
177  t = fad_jac_fill< Sacado::ELRFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
178  std::cout << "ELRSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
179  }
180  else if (num_eqns*2 == 64) {
181  t = fad_jac_fill< Sacado::ELRFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
182  std::cout << "ELRSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
183  }
184 
185  if (num_eqns*2 < slfad_max) {
186  t = fad_jac_fill< Sacado::ELRFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
187  std::cout << "ELRSLFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
188  }
189 
190  t = fad_jac_fill< Sacado::ELRFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
191  std::cout << "ELRDFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
192 
193  if (num_eqns*2 == 4) {
194  t = fad_jac_fill< Sacado::CacheFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
195  std::cout << "CacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
196  }
197  else if (num_eqns*2 == 16) {
198  t = fad_jac_fill< Sacado::CacheFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
199  std::cout << "CacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
200  }
201  else if (num_eqns*2 == 32) {
202  t = fad_jac_fill< Sacado::CacheFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
203  std::cout << "CacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
204  }
205  else if (num_eqns*2 == 64) {
206  t = fad_jac_fill< Sacado::CacheFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
207  std::cout << "CacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
208  }
209 
210  if (num_eqns*2 < slfad_max) {
211  t = fad_jac_fill< Sacado::CacheFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
212  std::cout << "CacheSLFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
213  }
214 
215  t = fad_jac_fill< Sacado::CacheFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
216  std::cout << "CacheFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
217 
218  if (num_eqns*2 == 4) {
219  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
220  std::cout << "ELRCacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
221  }
222  else if (num_eqns*2 == 16) {
223  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
224  std::cout << "ELRCacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
225  }
226  else if (num_eqns*2 == 32) {
227  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
228  std::cout << "ELRCacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
229  }
230  else if (num_eqns*2 == 64) {
231  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
232  std::cout << "ELRCacheSFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
233  }
234 
235  if (num_eqns*2 < slfad_max) {
236  t = fad_jac_fill< Sacado::ELRCacheFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
237  std::cout << "ELRCacheSLFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
238  }
239 
240  t = fad_jac_fill< Sacado::ELRCacheFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
241  std::cout << "ELRCacheFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
242 
243  t = fad_jac_fill< Sacado::Fad::DVFad<double> >(num_nodes, num_eqns, mesh_spacing);
244  std::cout << "DVFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
245 
246  }
247  catch (std::exception& e) {
248  std::cout << e.what() << std::endl;
249  ierr = 1;
250  }
251  catch (const char *s) {
252  std::cout << s << std::endl;
253  ierr = 1;
254  }
255  catch (...) {
256  std::cout << "Caught unknown exception!" << std::endl;
257  ierr = 1;
258  }
259 
260  return ierr;
261 }
MemPool * getMemoryPool(unsigned int dim)
Get memory pool for supplied dimension dim.
double residual_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
Derivative array storage class using dynamic memory allocation.
int main()
Definition: ad_example.cpp:191
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
void setDocString(const char doc_string[])
Forward-mode AD class using dynamic memory allocation and expression templates.