← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/sacado/test/performance/fad_fe_jac_fill.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/sacado/test/performance/fad_fe_jac_fill.cpp
 
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 
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.hpp"
 
33
#include "Sacado_Fad_SimpleFad.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
// ADOL-C includes
 
43
#ifdef HAVE_ADOLC
 
44
#ifdef PACKAGE
 
45
#undef PACKAGE
 
46
#endif
 
47
#ifdef PACKAGE_NAME
 
48
#undef PACKAGE_NAME
 
49
#endif
 
50
#ifdef PACKAGE_BUGREPORT
 
51
#undef PACKAGE_BUGREPORT
 
52
#endif
 
53
#ifdef PACKAGE_STRING
 
54
#undef PACKAGE_STRING
 
55
#endif
 
56
#ifdef PACKAGE_TARNAME
 
57
#undef PACKAGE_TARNAME
 
58
#endif
 
59
#ifdef PACKAGE_VERSION
 
60
#undef PACKAGE_VERSION
 
61
#endif
 
62
#ifdef VERSION
 
63
#undef VERSION
 
64
#endif
 
65
#define ADOLC_TAPELESS
 
66
#define NUMBER_DIRECTIONS 100
 
67
#include "adolc/adouble.h"
 
68
#include "adolc/drivers/drivers.h"
 
69
#include "adolc/interfaces.h"
 
70
#endif
 
71
 
 
72
// A performance test that computes a finite-element-like Jacobian using
 
73
// several Fad variants
 
74
 
 
75
template <>
 
76
Sacado::Fad::MemPool* Sacado::Fad::MemPoolStorage<double>::defaultPool_ = NULL;
 
77
 
 
78
void FAD::error(const char *msg) {
 
79
  std::cout << msg << std::endl;
 
80
}
 
81
 
 
82
struct ElemData {
 
83
  static const unsigned int nqp = 2;
 
84
  static const unsigned int nnode = 2;
 
85
  double w[nqp], jac[nqp], phi[nqp][nnode], dphi[nqp][nnode];
 
86
  unsigned int gid[nnode];
 
87
 
 
88
  ElemData(double mesh_spacing) {
 
89
    // Quadrature points
 
90
    double xi[nqp];
 
91
    xi[0] = -1.0 / std::sqrt(3.0);
 
92
    xi[1] =  1.0 / std::sqrt(3.0);
 
93
    
 
94
    for (unsigned int i=0; i<nqp; i++) {
 
95
      // Weights
 
96
      w[i] = 1.0;
 
97
      
 
98
      // Element transformation Jacobian
 
99
      jac[i] = 0.5*mesh_spacing;
 
100
      
 
101
      // Shape functions & derivatives
 
102
      phi[i][0] = 0.5*(1.0 - xi[i]);
 
103
      phi[i][1] = 0.5*(1.0 + xi[i]);
 
104
      dphi[i][0] = -0.5;
 
105
      dphi[i][1] = 0.5;
 
106
    }
 
107
  }
 
108
};
 
109
 
 
110
template <class FadType>
 
111
void fad_init_fill(const ElemData& e,
 
112
                   unsigned int neqn,
 
113
                   const std::vector<double>& x, 
 
114
                   std::vector<FadType>& x_fad) {
 
115
  for (unsigned int node=0; node<e.nnode; node++)
 
116
    for (unsigned int eqn=0; eqn<neqn; eqn++)
 
117
      x_fad[node*neqn+eqn] = FadType(e.nnode*neqn, node*neqn+eqn, 
 
118
                                     x[e.gid[node]*neqn+eqn]);
 
119
}
 
120
 
 
121
#ifdef HAVE_ADOLC
 
122
#ifndef ADOLC_TAPELESS
 
123
void adolc_init_fill(bool retape,
 
124
                     int tag,
 
125
                     const ElemData& e,
 
126
                     unsigned int neqn,
 
127
                     const std::vector<double>& x, 
 
128
                     std::vector<double>& x_local,
 
129
                     std::vector<adouble>& x_ad) {
 
130
  if (retape)
 
131
    trace_on(tag);
 
132
  for (unsigned int node=0; node<e.nnode; node++)
 
133
    for (unsigned int eqn=0; eqn<neqn; eqn++) {
 
134
      x_local[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
 
135
      if (retape)
 
136
        x_ad[node*neqn+eqn] <<= x_local[node*neqn+eqn];
 
137
    }
 
138
}
 
139
 
 
140
#else
 
141
 
 
142
void adolc_tapeless_init_fill(const ElemData& e,
 
143
                              unsigned int neqn,
 
144
                              const std::vector<double>& x, 
 
145
                              std::vector<adtl::adouble>& x_ad) {
 
146
  for (unsigned int node=0; node<e.nnode; node++)
 
147
    for (unsigned int eqn=0; eqn<neqn; eqn++) {
 
148
      x_ad[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
 
149
      for (unsigned int i=0; i<neqn*e.nnode; i++)
 
150
        x_ad[node*neqn+eqn].setADValue(i, 0.0);
 
151
      x_ad[node*neqn+eqn].setADValue(node*neqn+eqn, 1.0);
 
152
    }
 
153
  
 
154
}
 
155
#endif
 
156
#endif
 
157
 
 
158
void analytic_init_fill(const ElemData& e,
 
159
                        unsigned int neqn,
 
160
                        const std::vector<double>& x, 
 
161
                        std::vector<double>& x_local) {
 
162
  for (unsigned int node=0; node<e.nnode; node++) 
 
163
    for (unsigned int eqn=0; eqn<neqn; eqn++)
 
164
      x_local[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
 
165
}
 
166
 
 
167
template <class T>
 
168
void template_element_fill(const ElemData& e, 
 
169
                           unsigned int neqn,
 
170
                           const std::vector<T>& x, 
 
171
                           std::vector<T>& u, 
 
172
                           std::vector<T>& du, 
 
173
                           std::vector<T>& f) {
 
174
  // Construct element solution, derivative
 
175
  for (unsigned int qp=0; qp<e.nqp; qp++) {
 
176
    for (unsigned int eqn=0; eqn<neqn; eqn++) {
 
177
      u[qp*neqn+eqn] = 0.0;
 
178
      du[qp*neqn+eqn] = 0.0;
 
179
      for (unsigned int node=0; node<e.nnode; node++) {
 
180
        u[qp*neqn+eqn] += x[node*neqn+eqn] * e.phi[qp][node];
 
181
        du[qp*neqn+eqn] += x[node*neqn+eqn] * e.dphi[qp][node];
 
182
      }
 
183
    }
 
184
  }
 
185
 
 
186
  // Compute sum of equations for coupling
 
187
  std::vector<T> s(e.nqp*neqn);
 
188
  for (unsigned int qp=0; qp<e.nqp; qp++) {
 
189
    for (unsigned int eqn=0; eqn<neqn; eqn++) {
 
190
      s[qp*neqn+eqn] = 0.0;
 
191
      for (unsigned int j=0; j<neqn; j++) {
 
192
        if (j != eqn)
 
193
          s[qp*neqn+eqn] += u[qp*neqn+j]; 
 
194
      }
 
195
    }
 
196
  }
 
197
 
 
198
  // Evaluate element residual
 
199
  for (unsigned int node=0; node<e.nnode; node++) {
 
200
    for (unsigned int eqn=0; eqn<neqn; eqn++) {
 
201
      unsigned int row = node*neqn+eqn;
 
202
      f[row] = 0.0;
 
203
      for (unsigned int qp=0; qp<e.nqp; qp++) {
 
204
        f[row] += 
 
205
          e.w[qp]*e.jac[qp]*(-(1.0/(e.jac[qp]*e.jac[qp]))*
 
206
                             du[qp*neqn+eqn]*e.dphi[qp][node] + 
 
207
                             e.phi[qp][node]*(exp(u[qp*neqn+eqn]) + 
 
208
                                              u[qp*neqn+eqn]*s[qp*neqn+eqn]));
 
209
      }
 
210
    }
 
211
  }
 
212
}
 
213
 
 
214
void analytic_element_fill(const ElemData& e, 
 
215
                           unsigned int neqn,
 
216
                           const std::vector<double>& x, 
 
217
                           std::vector<double>& u, 
 
218
                           std::vector<double>& du, 
 
219
                           std::vector<double>& f,
 
220
                           std::vector< std::vector<double> >& jac) {
 
221
  // Construct element solution, derivative
 
222
  for (unsigned int qp=0; qp<e.nqp; qp++) {
 
223
    for (unsigned int eqn=0; eqn<neqn; eqn++) {
 
224
      u[qp*neqn+eqn] = 0.0;
 
225
      du[qp*neqn+eqn] = 0.0;
 
226
      for (unsigned int node=0; node<e.nnode; node++) {
 
227
        u[qp*neqn+eqn] += x[node*neqn+eqn] * e.phi[qp][node];
 
228
        du[qp*neqn+eqn] += x[node*neqn+eqn] * e.dphi[qp][node];
 
229
      }
 
230
    }
 
231
  }
 
232
 
 
233
  // Compute sum of equations for coupling
 
234
  std::vector<double> s(e.nqp*neqn);
 
235
  for (unsigned int qp=0; qp<e.nqp; qp++) {
 
236
    for (unsigned int eqn=0; eqn<neqn; eqn++) {
 
237
      s[qp*neqn+eqn] = 0.0;
 
238
      for (unsigned int j=0; j<neqn; j++) {
 
239
        if (j != eqn)
 
240
          s[qp*neqn+eqn] += u[qp*neqn+j]; 
 
241
      }
 
242
    }
 
243
  }
 
244
 
 
245
  // Evaluate element residual
 
246
  for (unsigned int node_row=0; node_row<e.nnode; node_row++) {
 
247
    for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
 
248
      unsigned int row = node_row*neqn+eqn_row;
 
249
      f[row] = 0.0;
 
250
      for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
 
251
        for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
 
252
          unsigned int col = node_col*neqn+eqn_col;
 
253
          jac[row][col] = 0.0;
 
254
        }
 
255
      }
 
256
      for (unsigned int qp=0; qp<e.nqp; qp++) {
 
257
        double c1 = e.w[qp]*e.jac[qp];
 
258
        double c2 = -(1.0/(e.jac[qp]*e.jac[qp]))*e.dphi[qp][node_row];
 
259
        double c3 = e.phi[qp][node_row]*(exp(u[qp*neqn+eqn_row]));
 
260
        double c4 = e.phi[qp][node_row]*u[qp*neqn+eqn_row];
 
261
        double c5 = e.phi[qp][node_row]*s[qp*neqn+eqn_row];
 
262
        double c35 = c3+c5;
 
263
        double c14 = c1*c4;
 
264
        f[row] += c1*(c2*du[qp*neqn+eqn_row] + c3 + c4*s[qp*neqn+eqn_row]);
 
265
        for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
 
266
          jac[row][node_col*neqn+eqn_row] += 
 
267
            c1*(c2*e.dphi[qp][node_col] + c35*e.phi[qp][node_col]);
 
268
          for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
 
269
            if (eqn_col != eqn_row)
 
270
              jac[row][node_col*neqn+eqn_col] += c14*e.phi[qp][node_col];
 
271
          }
 
272
        }
 
273
      }
 
274
    }
 
275
  }
 
276
}
 
277
 
 
278
template <class FadType>
 
279
void fad_process_fill(const ElemData& e,
 
280
                      unsigned int neqn,
 
281
                      const std::vector<FadType>& f_fad, 
 
282
                      std::vector<double>& f,
 
283
                      std::vector< std::vector<double> >& jac) {
 
284
  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
 
285
    f[e.gid[0]*neqn+eqn_row] += f_fad[0*neqn+eqn_row].val();
 
286
    f[e.gid[1]*neqn+eqn_row] += f_fad[1*neqn+eqn_row].val();
 
287
    for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
 
288
      for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
 
289
        unsigned int col = node_col*neqn+eqn_col;
 
290
        unsigned int next_col = (node_col+1)*neqn+eqn_col;
 
291
        jac[e.gid[0]*neqn+eqn_row][next_col] += 
 
292
          f_fad[0*neqn+eqn_row].fastAccessDx(col);
 
293
        jac[e.gid[1]*neqn+eqn_row][col] += 
 
294
          f_fad[1*neqn+eqn_row].fastAccessDx(col);
 
295
      }
 
296
    }
 
297
  }
 
298
}
 
299
 
 
300
#ifdef HAVE_ADOLC
 
301
#ifndef ADOLC_TAPELESS
 
302
void adolc_process_fill(bool retape,
 
303
                        int tag, 
 
304
                        const ElemData& e,
 
305
                        unsigned int neqn,
 
306
                        std::vector<double>& x_local,
 
307
                        std::vector<adouble>& f_ad, 
 
308
                        std::vector<double>& f,
 
309
                        std::vector<double>& f_local,
 
310
                        std::vector< std::vector<double> >& jac,
 
311
                        double **seed,
 
312
                        double **jac_local) {
 
313
  if (retape) {
 
314
    for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++)
 
315
      f_ad[0*neqn+eqn_row] >>= f_local[0*neqn+eqn_row];
 
316
    for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++)
 
317
      f_ad[1*neqn+eqn_row] >>= f_local[1*neqn+eqn_row];
 
318
    trace_off();
 
319
  }
 
320
 
 
321
  //jacobian(tag, neqn*e.nnode, neqn*e.nnode, &x_local[0], jac_local);
 
322
  forward(tag, neqn*e.nnode, neqn*e.nnode, neqn*e.nnode, &x_local[0],
 
323
          seed, &f_local[0], jac_local);
 
324
  
 
325
  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
 
326
    f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
 
327
    f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
 
328
    for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
 
329
      for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
 
330
        unsigned int col = node_col*neqn+eqn_col;
 
331
        unsigned int next_col = (node_col+1)*neqn+eqn_col;
 
332
        jac[e.gid[0]*neqn+eqn_row][next_col] += jac_local[0*neqn+eqn_row][col];
 
333
        jac[e.gid[1]*neqn+eqn_row][col] += jac_local[1*neqn+eqn_row][col];
 
334
      }
 
335
    }
 
336
  }
 
337
}
 
338
 
 
339
#else
 
340
 
 
341
void adolc_tapeless_process_fill(const ElemData& e,
 
342
                                 unsigned int neqn,
 
343
                                 const std::vector<adtl::adouble>& f_ad, 
 
344
                                 std::vector<double>& f,
 
345
                                 std::vector< std::vector<double> >& jac) {
 
346
  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
 
347
    f[e.gid[0]*neqn+eqn_row] += f_ad[0*neqn+eqn_row].getValue();
 
348
    f[e.gid[1]*neqn+eqn_row] += f_ad[1*neqn+eqn_row].getValue();
 
349
    for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
 
350
      for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
 
351
        unsigned int col = node_col*neqn+eqn_col;
 
352
        unsigned int next_col = (node_col+1)*neqn+eqn_col;
 
353
        jac[e.gid[0]*neqn+eqn_row][next_col] += 
 
354
          f_ad[0*neqn+eqn_row].getADValue(col);
 
355
        jac[e.gid[1]*neqn+eqn_row][col] += 
 
356
          f_ad[1*neqn+eqn_row].getADValue(col);
 
357
      }
 
358
    }
 
359
  }
 
360
}
 
361
#endif
 
362
#endif
 
363
 
 
364
void analytic_process_fill(const ElemData& e,
 
365
                           unsigned int neqn,
 
366
                           const std::vector<double>& f_local, 
 
367
                           const std::vector< std::vector<double> >& jac_local, 
 
368
                           std::vector<double>& f,
 
369
                           std::vector< std::vector<double> >& jac) {
 
370
  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
 
371
    f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
 
372
    f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
 
373
    for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
 
374
      for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
 
375
        unsigned int col = node_col*neqn+eqn_col;
 
376
        unsigned int next_col = (node_col+1)*neqn+eqn_col;
 
377
        jac[e.gid[0]*neqn+eqn_row][next_col] += jac_local[0*neqn+eqn_row][col];
 
378
        jac[e.gid[1]*neqn+eqn_row][col] += jac_local[1*neqn+eqn_row][col];
 
379
      }
 
380
    }
 
381
  }
 
382
}
 
383
 
 
384
void residual_process_fill(const ElemData& e,
 
385
                           unsigned int neqn,
 
386
                           const std::vector<double>& f_local, 
 
387
                           std::vector<double>& f) {
 
388
  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
 
389
    f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
 
390
    f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
 
391
  }
 
392
}
 
393
 
 
394
template <class FadType>
 
395
double fad_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
 
396
                    double mesh_spacing) {
 
397
  ElemData e(mesh_spacing);
 
398
 
 
399
  // Solution vector, residual, jacobian
 
400
  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
 
401
  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
 
402
  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
 
403
    for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
 
404
      unsigned int row = node_row*num_eqns + eqn_row;
 
405
      x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
 
406
      f[row] = 0.0;
 
407
      jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
 
408
      for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
 
409
        for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
 
410
          unsigned int col = node_col*num_eqns + eqn_col;
 
411
          jac[row][col] = 0.0;
 
412
        }
 
413
      }
 
414
    }
 
415
  }
 
416
 
 
417
  Teuchos::Time timer("FE Fad Jacobian Fill", false);
 
418
  timer.start(true);
 
419
  std::vector<FadType> x_fad(e.nnode*num_eqns), f_fad(e.nnode*num_eqns);
 
420
  std::vector<FadType> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
 
421
  for (unsigned int i=0; i<num_nodes-1; i++) {
 
422
    e.gid[0] = i;
 
423
    e.gid[1] = i+1;
 
424
    
 
425
    fad_init_fill(e, num_eqns, x, x_fad);
 
426
    template_element_fill(e, num_eqns, x_fad, u, du, f_fad);
 
427
    fad_process_fill(e, num_eqns, f_fad, f, jac);
 
428
  }
 
429
  timer.stop();
 
430
 
 
431
  // std::cout << "Fad Residual = " << std::endl;
 
432
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
 
433
  //   std::cout << "\t" << f[i] << std::endl;
 
434
 
 
435
  // std::cout.precision(8);
 
436
  // std::cout << "Fad Jacobian = " << std::endl;
 
437
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
 
438
  //   std::cout << "\t";
 
439
  //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
 
440
  //     std::cout << jac[i][j] << "\t";
 
441
  //   std::cout << std::endl;
 
442
  // }
 
443
 
 
444
  return timer.totalElapsedTime();
 
445
}
 
446
 
 
447
#ifdef HAVE_ADOLC
 
448
#ifndef ADOLC_TAPELESS
 
449
double adolc_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
 
450
                      double mesh_spacing) {
 
451
  ElemData e(mesh_spacing);
 
452
 
 
453
  // Solution vector, residual, jacobian
 
454
  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
 
455
  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
 
456
  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
 
457
    for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
 
458
      unsigned int row = node_row*num_eqns + eqn_row;
 
459
      x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
 
460
      f[row] = 0.0;
 
461
      jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
 
462
      for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
 
463
        for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
 
464
          unsigned int col = node_col*num_eqns + eqn_col;
 
465
          jac[row][col] = 0.0;
 
466
        }
 
467
      }
 
468
    }
 
469
  }
 
470
 
 
471
  Teuchos::Time timer("FE ADOL-C Jacobian Fill", false);
 
472
  timer.start(true);
 
473
  std::vector<adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
 
474
  std::vector<adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
 
475
  std::vector<double> x_local(e.nnode*num_eqns);
 
476
  std::vector<double> f_local(e.nnode*num_eqns);
 
477
  double **jac_local = new double*[e.nnode*num_eqns];
 
478
  double **seed = new double*[e.nnode*num_eqns];
 
479
  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
 
480
    jac_local[i] = new double[e.nnode*num_eqns];
 
481
    seed[i] = new double[e.nnode*num_eqns];
 
482
    for (unsigned int j=0; j<e.nnode*num_eqns; j++)
 
483
      seed[i][j] = 0.0;
 
484
    seed[i][i] = 1.0;
 
485
  }
 
486
  
 
487
  // Tape first element
 
488
  e.gid[0] = 0;
 
489
  e.gid[1] = 1;
 
490
  adolc_init_fill(true, 0, e, num_eqns, x, x_local, x_ad);
 
491
  template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
 
492
  adolc_process_fill(true, 0, e, num_eqns, x_local, f_ad, f, f_local,
 
493
                     jac, seed, jac_local);
 
494
 
 
495
  // Now do remaining fills reusing tape
 
496
  for (unsigned int i=1; i<num_nodes-1; i++) {
 
497
    e.gid[0] = i;
 
498
    e.gid[1] = i+1;
 
499
    
 
500
    adolc_init_fill(false, 0, e, num_eqns, x, x_local, x_ad);
 
501
    adolc_process_fill(false, 0, e, num_eqns, x_local, f_ad, f, f_local,
 
502
                       jac, seed, jac_local);
 
503
  }
 
504
  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
 
505
    delete [] jac_local[i];
 
506
    delete [] seed[i];
 
507
  }
 
508
  delete [] jac_local;
 
509
  delete [] seed;
 
510
  timer.stop();
 
511
 
 
512
  // std::cout << "ADOL-C Residual = " << std::endl;
 
513
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
 
514
  //   std::cout << "\t" << f[i] << std::endl;
 
515
 
 
516
  // std::cout.precision(8);
 
517
  // std::cout << "ADOL-C Jacobian = " << std::endl;
 
518
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
 
519
  //   std::cout << "\t";
 
520
  //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
 
521
  //     std::cout << jac[i][j] << "\t";
 
522
  //   std::cout << std::endl;
 
523
  // }
 
524
 
 
525
  return timer.totalElapsedTime();
 
526
}
 
527
 
 
528
double adolc_retape_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
 
529
                             double mesh_spacing) {
 
530
  ElemData e(mesh_spacing);
 
531
 
 
532
  // Solution vector, residual, jacobian
 
533
  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
 
534
  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
 
535
  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
 
536
    for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
 
537
      unsigned int row = node_row*num_eqns + eqn_row;
 
538
      x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
 
539
      f[row] = 0.0;
 
540
      jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
 
541
      for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
 
542
        for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
 
543
          unsigned int col = node_col*num_eqns + eqn_col;
 
544
          jac[row][col] = 0.0;
 
545
        }
 
546
      }
 
547
    }
 
548
  }
 
549
 
 
550
  Teuchos::Time timer("FE ADOL-C Retape Jacobian Fill", false);
 
551
  timer.start(true);
 
552
  std::vector<adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
 
553
  std::vector<adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
 
554
  std::vector<double> x_local(e.nnode*num_eqns);
 
555
  std::vector<double> f_local(e.nnode*num_eqns);
 
556
  double **jac_local = new double*[e.nnode*num_eqns];
 
557
  double **seed = new double*[e.nnode*num_eqns];
 
558
  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
 
559
    jac_local[i] = new double[e.nnode*num_eqns];
 
560
    seed[i] = new double[e.nnode*num_eqns];
 
561
    for (unsigned int j=0; j<e.nnode*num_eqns; j++)
 
562
      seed[i][j] = 0.0;
 
563
    seed[i][i] = 1.0;
 
564
  }
 
565
  for (unsigned int i=0; i<num_nodes-1; i++) {
 
566
    e.gid[0] = i;
 
567
    e.gid[1] = i+1;
 
568
    
 
569
    adolc_init_fill(true, 1, e, num_eqns, x, x_local, x_ad);
 
570
    template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
 
571
    adolc_process_fill(true, 1, e, num_eqns, x_local, f_ad, f, f_local,
 
572
                       jac, seed, jac_local);
 
573
  }
 
574
  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
 
575
    delete [] jac_local[i];
 
576
    delete [] seed[i];
 
577
  }
 
578
  delete [] jac_local;
 
579
  delete [] seed;
 
580
  timer.stop();
 
581
 
 
582
  // std::cout << "ADOL-C Residual (retaped) = " << std::endl;
 
583
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
 
584
  //   std::cout << "\t" << f[i] << std::endl;
 
585
 
 
586
  // std::cout.precision(8);
 
587
  // std::cout << "ADOL-C Jacobian (retaped) = " << std::endl;
 
588
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
 
589
  //   std::cout << "\t";
 
590
  //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
 
591
  //     std::cout << jac[i][j] << "\t";
 
592
  //   std::cout << std::endl;
 
593
  // }
 
594
 
 
595
  return timer.totalElapsedTime();
 
596
}
 
597
 
 
598
#else
 
599
 
 
600
double adolc_tapeless_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
 
601
                               double mesh_spacing) {
 
602
  ElemData e(mesh_spacing);
 
603
 
 
604
  // Solution vector, residual, jacobian
 
605
  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
 
606
  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
 
607
  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
 
608
    for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
 
609
      unsigned int row = node_row*num_eqns + eqn_row;
 
610
      x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
 
611
      f[row] = 0.0;
 
612
      jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
 
613
      for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
 
614
        for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
 
615
          unsigned int col = node_col*num_eqns + eqn_col;
 
616
          jac[row][col] = 0.0;
 
617
        }
 
618
      }
 
619
    }
 
620
  }
 
621
 
 
622
  Teuchos::Time timer("FE Tapeless ADOL-C Jacobian Fill", false);
 
623
  timer.start(true);
 
624
  std::vector<adtl::adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
 
625
  std::vector<adtl::adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
 
626
  for (unsigned int i=0; i<num_nodes-1; i++) {
 
627
    e.gid[0] = i;
 
628
    e.gid[1] = i+1;
 
629
    
 
630
    adolc_tapeless_init_fill(e, num_eqns, x, x_ad);
 
631
    template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
 
632
    adolc_tapeless_process_fill(e, num_eqns, f_ad, f, jac);
 
633
  }
 
634
  timer.stop();
 
635
 
 
636
  // std::cout << "Tapeless ADOL-C Residual = " << std::endl;
 
637
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
 
638
  //   std::cout << "\t" << f[i] << std::endl;
 
639
 
 
640
  // std::cout.precision(8);
 
641
  // std::cout << "Tapeless ADOL-C Jacobian = " << std::endl;
 
642
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
 
643
  //   std::cout << "\t";
 
644
  //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
 
645
  //     std::cout << jac[i][j] << "\t";
 
646
  //   std::cout << std::endl;
 
647
  // }
 
648
 
 
649
  return timer.totalElapsedTime();
 
650
}
 
651
#endif
 
652
#endif
 
653
 
 
654
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
 
655
                         double mesh_spacing) {
 
656
  ElemData e(mesh_spacing);
 
657
 
 
658
  // Solution vector, residual, jacobian
 
659
  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
 
660
  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
 
661
  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
 
662
    for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
 
663
      unsigned int row = node_row*num_eqns + eqn_row;
 
664
      x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
 
665
      f[row] = 0.0;
 
666
      jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
 
667
      for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
 
668
        for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
 
669
          unsigned int col = node_col*num_eqns + eqn_col;
 
670
          jac[row][col] = 0.0;
 
671
        }
 
672
      }
 
673
    }
 
674
  }
 
675
 
 
676
  Teuchos::Time timer("FE Analytic Jacobian Fill", false);
 
677
  timer.start(true);
 
678
  std::vector<double> x_local(e.nnode*num_eqns), f_local(e.nnode*num_eqns);
 
679
  std::vector< std::vector<double> > jac_local(e.nnode*num_eqns);
 
680
  std::vector<double> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
 
681
  for (unsigned int i=0; i<e.nnode*num_eqns; i++)
 
682
    jac_local[i] = std::vector<double>(e.nnode*num_eqns);
 
683
  for (unsigned int i=0; i<num_nodes-1; i++) {
 
684
    e.gid[0] = i;
 
685
    e.gid[1] = i+1;
 
686
    
 
687
    analytic_init_fill(e, num_eqns, x, x_local);
 
688
    analytic_element_fill(e, num_eqns, x_local, u, du, f_local, jac_local);
 
689
    analytic_process_fill(e, num_eqns, f_local, jac_local, f, jac);
 
690
  }
 
691
  timer.stop();
 
692
 
 
693
  // std::cout.precision(8);
 
694
  // std::cout << "Analytic Residual = " << std::endl;
 
695
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
 
696
  //   std::cout << "\t" << f[i] << std::endl;
 
697
 
 
698
  // std::cout.precision(8);
 
699
  // std::cout << "Analytic Jacobian = " << std::endl;
 
700
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
 
701
  //   std::cout << "\t";
 
702
  //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
 
703
  //     std::cout << jac[i][j] << "\t";
 
704
  //   std::cout << std::endl;
 
705
  // }
 
706
 
 
707
  return timer.totalElapsedTime();
 
708
}
 
709
 
 
710
double residual_fill(unsigned int num_nodes, unsigned int num_eqns,
 
711
                         double mesh_spacing) {
 
712
  ElemData e(mesh_spacing);
 
713
 
 
714
  // Solution vector, residual, jacobian
 
715
  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
 
716
  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
 
717
    for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
 
718
      unsigned int row = node_row*num_eqns + eqn_row;
 
719
      x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
 
720
      f[row] = 0.0;
 
721
    }
 
722
  }
 
723
 
 
724
  Teuchos::Time timer("FE Residual Fill", false);
 
725
  timer.start(true);
 
726
  std::vector<double> x_local(e.nnode*num_eqns), f_local(e.nnode*num_eqns);
 
727
  std::vector<double> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
 
728
  for (unsigned int i=0; i<num_nodes-1; i++) {
 
729
    e.gid[0] = i;
 
730
    e.gid[1] = i+1;
 
731
    
 
732
    analytic_init_fill(e, num_eqns, x, x_local);
 
733
    template_element_fill(e, num_eqns, x_local, u, du, f_local);
 
734
    residual_process_fill(e, num_eqns, f_local, f);
 
735
  }
 
736
  timer.stop();
 
737
 
 
738
  // std::cout.precision(8);
 
739
  // std::cout << "Analytic Residual = " << std::endl;
 
740
  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
 
741
  //   std::cout << "\t" << f[i] << std::endl;
 
742
 
 
743
  return timer.totalElapsedTime();
 
744
}
 
745
 
 
746
int main(int argc, char* argv[]) {
 
747
  int ierr = 0;
 
748
 
 
749
  try {
 
750
    double t, ta, tr;
 
751
    int p = 2;
 
752
    int w = p+7;
 
753
 
 
754
    // Maximum number of derivative components for SLFad
 
755
    const int slfad_max = 100;
 
756
 
 
757
    // Set up command line options
 
758
    Teuchos::CommandLineProcessor clp;
 
759
    clp.setDocString("This program tests the speed of various forward mode AD implementations for a finite-element-like Jacobian fill");
 
760
    int num_nodes = 100000;
 
761
    int num_eqns = 2;
 
762
    int rt = 0;
 
763
    clp.setOption("n", &num_nodes, "Number of nodes");
 
764
    clp.setOption("p", &num_eqns, "Number of equations");
 
765
    clp.setOption("rt", &rt, "Include ADOL-C retaping test");
 
766
 
 
767
    // Parse options
 
768
    Teuchos::CommandLineProcessor::EParseCommandLineReturn
 
769
      parseReturn= clp.parse(argc, argv);
 
770
    if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
 
771
      return 1;
 
772
 
 
773
    double mesh_spacing = 1.0 / static_cast<double>(num_nodes - 1);
 
774
 
 
775
    // Memory pool & manager
 
776
    Sacado::Fad::MemPoolManager<double> poolManager(num_nodes*num_eqns);
 
777
    Sacado::Fad::MemPool* pool = poolManager.getMemoryPool(num_nodes*num_eqns);
 
778
    Sacado::Fad::DMFad<double>::setDefaultPool(pool);
 
779
 
 
780
    std::cout.setf(std::ios::scientific);
 
781
    std::cout.precision(p);
 
782
    std::cout << "num_nodes =  " << num_nodes 
 
783
              << ", num_eqns = " << num_eqns << ":  " << std::endl
 
784
              << "           " << "   Time" << "\t"<< "Time/Analytic" << "\t"
 
785
              << "Time/Residual" << std::endl;
 
786
 
 
787
    ta = 1.0;
 
788
    tr = 1.0;
 
789
 
 
790
    tr = residual_fill(num_nodes, num_eqns, mesh_spacing);
 
791
 
 
792
    ta = analytic_jac_fill(num_nodes, num_eqns, mesh_spacing);
 
793
    std::cout << "Analytic:  " << std::setw(w) << ta << "\t" << std::setw(w) << ta/ta << "\t" << std::setw(w) << ta/tr << std::endl;
 
794
 
 
795
#ifdef HAVE_ADOLC
 
796
#ifndef ADOLC_TAPELESS
 
797
    t = adolc_jac_fill(num_nodes, num_eqns, mesh_spacing);
 
798
    std::cout << "ADOL-C:    " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
799
 
 
800
    if (rt != 0) {
 
801
      t = adolc_retape_jac_fill(num_nodes, num_eqns, mesh_spacing);
 
802
      std::cout << "ADOL-C(rt):" << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
803
    }
 
804
 
 
805
#else
 
806
    t = adolc_tapeless_jac_fill(num_nodes, num_eqns, mesh_spacing);
 
807
    std::cout << "ADOL-C(tl):" << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
808
#endif
 
809
#endif
 
810
 
 
811
    if (num_eqns*2 == 4) {
 
812
      t = fad_jac_fill< FAD::TFad<4,double> >(num_nodes, num_eqns, mesh_spacing);
 
813
      std::cout << "TFad:      " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
814
    }
 
815
    else if (num_eqns*2 == 40) {
 
816
      t = fad_jac_fill< FAD::TFad<40,double> >(num_nodes, num_eqns, mesh_spacing);
 
817
      std::cout << "TFad:      " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
818
    }
 
819
 
 
820
    t = fad_jac_fill< FAD::Fad<double> >(num_nodes, num_eqns, mesh_spacing);
 
821
    std::cout << "Fad:       " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
822
 
 
823
    if (num_eqns*2 == 4) {
 
824
      t = fad_jac_fill< Sacado::Fad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
 
825
      std::cout << "SFad:      " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
826
    }
 
827
    else if (num_eqns*2 == 40) {
 
828
      t = fad_jac_fill< Sacado::Fad::SFad<double,40> >(num_nodes, num_eqns, mesh_spacing);
 
829
      std::cout << "SFad:      " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
830
    }
 
831
 
 
832
    if (num_eqns*2 < slfad_max) {
 
833
      t = fad_jac_fill< Sacado::Fad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
 
834
      std::cout << "SLFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
835
    }
 
836
    
 
837
    t = fad_jac_fill< Sacado::Fad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
 
838
    std::cout << "DFad:      " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
839
 
 
840
    t = fad_jac_fill< Sacado::Fad::SimpleFad<double> >(num_nodes, num_eqns, mesh_spacing);
 
841
    std::cout << "SimpleFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
842
 
 
843
    t = fad_jac_fill< Sacado::Fad::DMFad<double> >(num_nodes, num_eqns, mesh_spacing);
 
844
    std::cout << "DMFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl; 
 
845
 
 
846
    if (num_eqns*2 == 4) {
 
847
      t = fad_jac_fill< Sacado::ELRFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
 
848
      std::cout << "ELRSFad:   " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
849
    }
 
850
    else if (num_eqns*2 == 40) {
 
851
      t = fad_jac_fill< Sacado::ELRFad::SFad<double,40> >(num_nodes, num_eqns, mesh_spacing);
 
852
      std::cout << "ELRSFad:   " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
853
    }
 
854
 
 
855
    if (num_eqns*2 < slfad_max) {
 
856
      t = fad_jac_fill< Sacado::ELRFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
 
857
      std::cout << "ELRSLFad:  " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
858
    }
 
859
 
 
860
    t = fad_jac_fill< Sacado::ELRFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
 
861
    std::cout << "ELRDFad:   " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
862
    
 
863
    t = fad_jac_fill< Sacado::CacheFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
 
864
    std::cout << "CacheFad:  " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
865
 
 
866
    t = fad_jac_fill< Sacado::Fad::DVFad<double> >(num_nodes, num_eqns, mesh_spacing);
 
867
    std::cout << "DVFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/tr << std::endl;
 
868
    
 
869
  }
 
870
  catch (std::exception& e) {
 
871
    cout << e.what() << endl;
 
872
    ierr = 1;
 
873
  }
 
874
  catch (const char *s) {
 
875
    cout << s << endl;
 
876
    ierr = 1;
 
877
  }
 
878
  catch (...) {
 
879
    cout << "Caught unknown exception!" << endl;
 
880
    ierr = 1;
 
881
  }
 
882
 
 
883
  return ierr;
 
884
}