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 void FAD::error(const char *msg) {
41  std::cout << msg << std::endl;
42 }
43 
44 int main(int argc, char* argv[]) {
45  int ierr = 0;
46 
47  try {
48  double t, ta, tr;
49  int p = 2;
50  int w = p+7;
51 
52  // Maximum number of derivative components for SLFad
53  const int slfad_max = 130;
54 
55  // Set up command line options
56  Teuchos::CommandLineProcessor clp;
57  clp.setDocString("This program tests the speed of various forward mode AD implementations for a finite-element-like Jacobian fill");
58  int num_nodes = 100000;
59  int num_eqns = 2;
60  int rt = 0;
61  clp.setOption("n", &num_nodes, "Number of nodes");
62  clp.setOption("p", &num_eqns, "Number of equations");
63  clp.setOption("rt", &rt, "Include ADOL-C retaping test");
64 
65  // Parse options
66  Teuchos::CommandLineProcessor::EParseCommandLineReturn
67  parseReturn= clp.parse(argc, argv);
68  if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
69  return 1;
70 
71  double mesh_spacing = 1.0 / static_cast<double>(num_nodes - 1);
72 
73  std::cout.setf(std::ios::scientific);
74  std::cout.precision(p);
75  std::cout << "num_nodes = " << num_nodes
76  << ", num_eqns = " << num_eqns << ": " << std::endl
77  << " " << " Time " << "\t"<< "Time/Analytic" << "\t"
78  << "Time/(2*p*Residual)" << std::endl;
79 
80  ta = 1.0;
81  tr = 1.0;
82 
83  tr = residual_fill(num_nodes, num_eqns, mesh_spacing);
84 
85  ta = analytic_jac_fill(num_nodes, num_eqns, mesh_spacing);
86  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;
87 
88 #ifdef HAVE_ADOLC
89 #ifndef ADOLC_TAPELESS
90  t = adolc_jac_fill(num_nodes, num_eqns, mesh_spacing);
91  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;
92 
93  if (rt != 0) {
94  t = adolc_retape_jac_fill(num_nodes, num_eqns, mesh_spacing);
95  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;
96  }
97 
98 #else
99  t = adolc_tapeless_jac_fill(num_nodes, num_eqns, mesh_spacing);
100  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;
101 #endif
102 #endif
103 
104 #ifdef HAVE_ADIC
105  t = adic_jac_fill(num_nodes, num_eqns, mesh_spacing);
106  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;
107 #endif
108 
109  if (num_eqns*2 == 4) {
110  t = fad_jac_fill< FAD::TFad<16,double> >(num_nodes, num_eqns, mesh_spacing);
111  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;
112  }
113  else if (num_eqns*2 == 16) {
114  t = fad_jac_fill< FAD::TFad<16,double> >(num_nodes, num_eqns, mesh_spacing);
115  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;
116  }
117  else if (num_eqns*2 == 32) {
118  t = fad_jac_fill< FAD::TFad<32,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 == 64) {
122  t = fad_jac_fill< FAD::TFad<64,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 
126  t = fad_jac_fill< FAD::Fad<double> >(num_nodes, num_eqns, mesh_spacing);
127  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;
128 
129  if (num_eqns*2 == 4) {
130  t = fad_jac_fill< Sacado::Fad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
131  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;
132  }
133  else if (num_eqns*2 == 16) {
134  t = fad_jac_fill< Sacado::Fad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
135  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;
136  }
137  else if (num_eqns*2 == 32) {
138  t = fad_jac_fill< Sacado::Fad::SFad<double,32> >(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 == 64) {
142  t = fad_jac_fill< Sacado::Fad::SFad<double,64> >(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 
146  if (num_eqns*2 < slfad_max) {
147  t = fad_jac_fill< Sacado::Fad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
148  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;
149  }
150 
151  t = fad_jac_fill< Sacado::Fad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
152  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;
153 
154  t = fad_jac_fill< Sacado::Fad::SimpleFad<double> >(num_nodes, num_eqns, mesh_spacing);
155  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;
156 
157  if (num_eqns*2 == 4) {
158  t = fad_jac_fill< Sacado::ELRFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
159  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;
160  }
161  else if (num_eqns*2 == 16) {
162  t = fad_jac_fill< Sacado::ELRFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
163  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;
164  }
165  else if (num_eqns*2 == 32) {
166  t = fad_jac_fill< Sacado::ELRFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
167  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;
168  }
169  else if (num_eqns*2 == 64) {
170  t = fad_jac_fill< Sacado::ELRFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
171  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;
172  }
173 
174  if (num_eqns*2 < slfad_max) {
175  t = fad_jac_fill< Sacado::ELRFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
176  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;
177  }
178 
179  t = fad_jac_fill< Sacado::ELRFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
180  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;
181 
182  if (num_eqns*2 == 4) {
183  t = fad_jac_fill< Sacado::CacheFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
184  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;
185  }
186  else if (num_eqns*2 == 16) {
187  t = fad_jac_fill< Sacado::CacheFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
188  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;
189  }
190  else if (num_eqns*2 == 32) {
191  t = fad_jac_fill< Sacado::CacheFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
192  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;
193  }
194  else if (num_eqns*2 == 64) {
195  t = fad_jac_fill< Sacado::CacheFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
196  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;
197  }
198 
199  if (num_eqns*2 < slfad_max) {
200  t = fad_jac_fill< Sacado::CacheFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
201  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;
202  }
203 
204  t = fad_jac_fill< Sacado::CacheFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
205  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;
206 
207  if (num_eqns*2 == 4) {
208  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
209  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;
210  }
211  else if (num_eqns*2 == 16) {
212  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
213  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;
214  }
215  else if (num_eqns*2 == 32) {
216  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
217  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;
218  }
219  else if (num_eqns*2 == 64) {
220  t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
221  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;
222  }
223 
224  if (num_eqns*2 < slfad_max) {
225  t = fad_jac_fill< Sacado::ELRCacheFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
226  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;
227  }
228 
229  t = fad_jac_fill< Sacado::ELRCacheFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
230  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;
231 
232  t = fad_jac_fill< Sacado::Fad::DVFad<double> >(num_nodes, num_eqns, mesh_spacing);
233  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;
234 
235  }
236  catch (std::exception& e) {
237  std::cout << e.what() << std::endl;
238  ierr = 1;
239  }
240  catch (const char *s) {
241  std::cout << s << std::endl;
242  ierr = 1;
243  }
244  catch (...) {
245  std::cout << "Caught unknown exception!" << std::endl;
246  ierr = 1;
247  }
248 
249  return ierr;
250 }
residual_fill
double residual_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
Definition: fe_jac_fill_funcs.cpp:183
analytic_jac_fill
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
Definition: fe_jac_fill_funcs.cpp:116
Teuchos::CommandLineProcessor::setOption
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
fe_jac_fill_funcs.hpp
main
int main()
Definition: ad_example.cpp:191
Teuchos::CommandLineProcessor::parse
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
Teuchos::CommandLineProcessor::setDocString
void setDocString(const char doc_string[])
Generated by   doxygen 1.8.5