Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
40void FAD::error(const char *msg) {
41 std::cout << msg << std::endl;
42}
43
44int 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}
int main()
Definition: ad_example.cpp:191
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)
const char * p