# Presentation: 137389

## Automatic Differentiation of C++ Codes With Sacado

- Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company,for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

- Eric Phipps and David Gay
- Sandia National Laboratories
- Trilinos User’s Group Meeting
- November 7, 2006

**Automatic Differentiation of C++ Codes With Sacado**

## Outline

- Introduction to automatic differentiation
- Forward mode
- Reverse mode
- Taylor polynomial mode
- Software implementations
- Source Transformation
- Operator Overloading
- Sacado
- Differentiating large-scale element-based codes

## What Is Automatic Differentation (AD) ?

- Analytic derivatives without hand-coding
- All differentiable computations are composition of simple operations
- sin(), log(), +, *, /, etc...
- We know the derivatives of these simple operations
- We have the chain rule from calculus
- Systematic application of the chain rule through your computation differentiating each statement line-by-line.

## A Simple Example

**A Simple Example**

|  |  |
| --- | --- |
| 2.000 | 1.000 |
| 7.389 | 7.389 |
| 0.693 | 0.500 |
| 1.386 | 1.693 |
| 8.775 | 9.082 |
| 0.605 | -7.233 |

|  |
| --- |
| 2.000 |
| 7.389 |
| 0.693 |
| 1.386 |
| 8.775 |
| 0.605 |

- Analytic derivative evaluated to machine precision

## Related Methods

- Automatic Differentiation

- Symbolic Differentiation

- Finite Differencing

## Why is this important?

- We need analytic first & higher derivatives for predictive simulations
- Computational design, optimization and parameter estimation
- Stability analysis
- Uncertainty quantification
- Verification and validation
- Analytic derivatives improve robustness and efficiency
- Infeasible to expect application developers to code analytic derivatives
- Time consuming, error prone, and difficult to verify
- Thousands of possible parameters in a large code
- Developers must understand what derivatives are needed
- Automatic differentiation solves these problems

**Notes:**
So why does Sandia need analytic derivatives to achieve its missions?
Major emphasis on developing predictive simulation capabilities
Which means we focus more and more on optimization, ..., and derivatives are a critical piece of the “glue” that ties these methods to simulation codes.
Typically, Sandia has relied on FD for much of this work
But this is inefficient, especially when derivatives w.r.t. large number of quantities are needed, and not robust, particularly in the presence of ill-conditioning or when higher derivatives are needed.
To make progress towards predictive simulations, we need analytic derivatives.
Not practical to expect code developers to code analytic derivatives
Takes too long and it is too easy to make mistakes
Especially considering that there can be thousands of parameters in a code you might want sensitivities of for UQ or validation, for example
But most importantly because it forces app developers to understand what derivatives are needed for a diverse set of algorithms, both now and in the future
Automatic differentiation is a way to get these analytic derivatives without these problems

Transition:  How is this actually implemented

## Tangent Propagation

**Tangent Propagation**

- Tangents
- For each intermediate operation
- Tangents map forward through evaluation

| Operation | Tangent Rule |
| --- | --- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

**Notes:**
Y_dot depends only only on x_0, x_dot.

## A Simple Tangent Example

**Notes:**
Tangents can be evaluated along with original function evaluation, and each tangent computation only involves values, tangents for its arguments

## Forward Mode AD via Tangent Propagation

- Choice of space curve        is arbitrary
- Tangent    depends only on     ,
- Given     and   :
- Propagate     vectors                  simultaneously
- Forward mode AD:
-     is called the seed matrix.  Setting equal to identity matrix yields full Jacobian
- Computational cost
- Jacobian-vector products, directional derivatives, Jacobians for

- Jacobian vector product

- Jacobian matrix product

**Notes:**
x_i is the ith indep. var, tangent components are initialized as ith row of V

## Gradient Propagation

- Gradients
- For each intermediate operation
- Gradients map backward through evaluation

| Operation | Gradient Rule |
| --- | --- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

## A Simple Gradient Example

**Notes:**
Gradient evaluation is fundamentally different than tangent propagation
Gradient for s_1 depends on gradient for s_3, but that can’t be evaluated until the next operation
Requires buffering results of all intermediate operations (all of the s_i’s) to sum buffere (called the tape) 

## Reverse Mode AD via Gradient Propagation

- Choice of normal     is arbitrary
- Gradient    depends only on     ,
- Given     and   :
- Propagate     vectors                   simultaneously
- Reverse mode AD:
-     is called the seed matrix.  Setting equal to identity matrix yields full Jacobian
- Computational cost
- Jacobian-transpose products, gradients, Jacobians for

- Jacobian-transpose vector product

- Jacobian-transpose matrix product

## Taylor Polynomial Propagation

- Extension of tangent propagation to higher degree
- Given            coefficients
- Computational cost

| Operation | Taylor Rule |
| --- | --- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

## Software Implementations

- Source transformation
- Preprocessor reads code to be differentiated, uses AD to generate derivative code, writes-out differentiated code in original source language which is then compiled using a standard compiler
- Resulting derivative computation is usually very efficient
- Works well for simple languages (FORTRAN, some C)
- ADIFOR, ADIC out of Argonne
- Extremely difficult for C++ (no existing tool)
- Operator overloading
- New data types are created for forward, reverse, Taylor modes, and intrinsic operations/elementary operations are overloaded to compute derivatives as a side-effect
- Generally easy to incorporate into C++ codes
- Generally slower than source transformation due to function call overhead
- Requires changing data types from floats/doubles to AD types
- C++ templates greatly help
- ADOL-C (slow), FAD/TFAD (fast)

## ADIFOR* Example

**      ****subroutine g_func(g_p_, x, g_x, ldg_x, y, g_y, ldg_y)**
**C       Initializations removed for clarity...**
**        **
**        ****d2_v = exp(x(1))**
**        ****d1_p =  d2_v**
**        ****do g_i_ = 1, g_p_**
**          ****g_u(g_i_) = d1_p * g_x(g_i_, 1)**
**        ****enddo**
**        ****u = d2_v**
**C--------**
**        ****do g_i_ = 1, g_p_**
**          ****g_v(g_i_) = x(1) * g_x(g_i_, 2) + x(2) * g_x(g_i_, 1)**
**        ****enddo**
**        ****v = x(1) * x(2)**
**C--------**
**        ****do g_i_ = 1, g_p_**
**          ****g_w(g_i_) = g_v(g_i_) + g_u(g_i_)**
**        ****enddo**
**        ****w = u + v**
**C--------**
**        ****d2_v = sin(w)**
**        ****d1_p = cos(w)**
**        ****do g_i_ = 1, g_p_**
**          ****g_y(g_i_, 1) = d1_p * g_w(g_i_)**
**        ****enddo**
**        ****y(1) = d2_v**
**C      continues...**

**ADIFOR********* Example**

**       ****subroutine func(x, y)**
**C      **
**       ****double precision x(2), y(2)**
**       ****double precision u, v, w**
**C      **
**       ****u = exp(x(1))**
**       ****v = x(1)*x(2)**
**       ****w = u+v**
**       ****y(1) = sin(w)**
**C      **
**       ****u = x(1)**2**
**       ****v = y(1) + u**
**       ****y(2) = y(1)/v**
**C      **
**       ****return**
**       ****end**

- *ADIFOR 2.0D
- www-unix.mcs.anl.gov/autodiff/ADIFOR/

## (Naive) Operator Overloading Example

**class Tangent {**
**public:**
**  ****static const int N = 2;**
**  ****double val;**
**  ****double dot[N];**
**};**
**Tangent operator+(const Tangent& a, const Tangent& b) {**
**  ****Tangent c;**
**  ****c.val = a.val + b.val;**
**  ****for (int i=0; i<Tangent::N; i++) **
**    ****c.dot[i] = a.dot[i] + b.dot[i];**
**	****return c;**
**}**
**Tangent operator*(const Tangent& a, const Tangent& b) {**
**  ****Tangent c;**
**  ****c.val = a.val * b.val;**
**  ****for (int i=0; i<Tangent::N; i++) **
**    ****c.dot[i] = a.val * b.dot[i] + a.dot[i]*b.val;**
**	****return c;**
**}**
**Tangent exp(const Tangent& a) {**
**  ****Tangent c;**
**  ****c.val = exp(a.val);**
**  ****for (int i=0; i<Tangent::N; i++) **
**    ****c.dot[i] = c.val * a.dot[i];**
**	****return c;**
**}**

**(Naive) Operator Overloading Example**

**void func(const double x[], double y[]) {**
**  ****double u, v, w;**
**  ****u = exp(x[0]);**
**  ****v = x[0]*x[1];**
**  ****w = u+v;**
**  ****y[0] = sin(w);**
**     **
**  ****u = x[0]*x[0];**
**  ****v = y[0] + u;**
**  ****y[1] = y[0]/v;**
**}**
**void func(const Tangent x[], Tangent y[]) {**
**  ****Tangent u, v, w;**
**  ****u = exp(x[0]);**
**  ****v = x[0]*x[1];**
**  ****w = u+v;**
**  ****y[0] = sin(w);**
**     **
**  ****u = x[0]*x[0];**
**  ****v = y[0] + u;**
**  ****y[1] = y[0]/v;**
**}**

**Notes:**
Chose static arrays for simplicity
Difference between versions of func is just types
Not very efficient due to copies

## Expression Template Operator Overloading

**template <class E1, E2> class PlusExpr {**
**  ****double val() const { return e1.val() + e2.val(); }**
**  ****double dx(int i) const { return e1.dx(i) + e2.dx(i); }**
**  ****const E1& e1;**
**  ****const E2& e2;**
**};**
**template<class E1, class E2> PlusExpr<E1,E2> **
**operator+(const E1& a, const E2& b) {**
**  ****return PlusExpr<E1,E2>(a,b); **
**}**
**template <class E1> class SinExpr {**
**  ****double val() const { return sin(e1.val())] }**
**  ****double dx(int i) const { return cos(e1.val())*e1.dx(i); }**
**  ****const E1& e1;**
**};**
**template<class E1> SinExpr<E1> sin(const E1& a) {**
**  ****return SinExpr<E1>(a); **
**}**
**class Tangent {**
**public:**
**  ****double val() const { return val; }**
**  ****double dx(int i) const { return dot[i]; }**
**  ****template <class E> Tangent& operator=(const E& e) {**
**    ****val = e.val();**
**    ****for (int i=0; i<N; i++) **
**      ****dot[i] = e.dx(i); **
**   ****}**
**};**

**Expression Template Operator Overloading**

**void func(const Tangent x[], Tangent y[]) {**
**  ****y[0] = sin(exp(x[0]) + x[0]*x[1]);**
**  ****//...**
**}**
**SinExpr< PlusExp< ExpExpr<Tangent>,**
**                  ****MultExpr<Tangent,Tangent> **
**                 ****> **
**        ****>**
**y[0].val = sin(exp(x[0]) + x[0]*x[1]);**
**for (int i=0; i<N; i++) {**
**  ****y[0].dot[i] = cos(exp(x[0]) + x[0]*x[1])***
**     ****(exp(x[0])*x[0].dot[i] +    **
**      ****x[0]*x[1].dot[i] + x[1]*x[0].dot[i]);**

- Public domain Fad/TFad package

## Introducing Sacado (The Package Formerly Known as ADTools)

- New Trilinos package for automatic differentiation of C++ codes
- Loosely translated as “I have derived” in Spanish
- Developers:  Dave Gay, Eric Phipps (with contributions from Ross Bartlett)
- Not in Trilinos 7, will be released with Trilinos 8 (Spring ‘07)
- Designed for fast differentiation of *small* computations using C++ templates
- Forward AD based on expression-templated public domain Fad package
- Dave Gay’s Reverse AD package Rad
- Univariate Taylor polynomials
- Application support utilities
- Coming soon
- Ross’ ScalarFlopCounter
- Multi-variate Taylor polynomials
- Stochastic polynomial chaos expansions
- Teuchos BLAS/LAPACK wrapper specializations
- Depends on Teuchos only (and currently only through tests)

## The Usual Suspects

- Configure options
- --enable-sacado — Enables Sacado at Trilinos top-level
- --enable-sacado-tests, --enable-tests — Enables unit, regression, and performance tests
- 		--with-cppunit-prefix=[path] — Path to CppUnit for unit tests
- 		--with-adolc=[path] — Enables Taylor polynomial unit tests with ADOL-C
- 		--enable-sacado-alltests — Enables additional tests that take a VERY LONG 			time to compile
- --enable-sacado-examples, --enable-examples — Enables examples
- 		--enable-sacado-fem-example — Enables a 1D FEM example application 			(additional dependencies on Epetra, NOX, LOCA, MOOCHO, Rythmos, ...)
- Mailing lists
- _[Sacado-announce@software.sandia.gov](mailto:Sacado-announce@software.sandia.gov)_
- _[Sacado-checkins@software.sandia.gov](mailto:Sacado-checkins@software.sandia.gov)_
- _[Sacado-developers@software.sandia.gov](mailto:Sacado-developers@software.sandia.gov)_
- _[Sacado-regression@software.sandia.gov](mailto:Sacado-regression@software.sandia.gov)_
- _[Sacado-users@software.sandia.gov](mailto:Sacado-users@software.sandia.gov)_
- Bugzilla:  _[http://software.sandia.gov/bugzilla](http://software.sandia.gov/bugzilla)_
- Bonsai:  _[http://software.sandia.gov/bonsai/cvsqueryform.cgi](http://software.sandia.gov/bonsai/cvsqueryform.cgi)_
- Web:  _[http://software.sandai.gov/Trilinos/packages/sacado](http://software.sandai.gov/Trilinos/packages/sacado)_
- Doxygen documentation

## Using Sacado

- As always:  #include “Sacado.hpp”
- All classes are templated on the Scalar type
- Forward AD classes:
- Sacado::Fad::DFad<ScalarT>:  Derivative array is allocated dynamically
- Sacado::Fad::SFad<ScalarT>:  Derivative array is allocated statically and dimension must be known at compile time
- Sacado::Fad::SLFad<ScalarT>:  Like SFad except allocated length may be greater than “used” length
- Sacado::Fad::DMFad<ScalarT>:  Dynamically allocated array using a memory pool (Doesn’t exist yet)
- Reverse mode AD classes:
- ADvar<ScalarT> (Sacado_trad.h)
- Taylor polynomial classes:
- Sacado::Taylor::DTaylor<ScalarT>

## sacado/example/dfad_example.cpp

- #include "Sacado.hpp"
- // The function to differentiate
- template <typename ScalarT>
- ScalarT func(const ScalarT& a, const ScalarT& b, const ScalarT& c) {
-   ScalarT r = c*std::log(b+1.)/std::sin(a);
-   return r;
- }
- int main(int argc, char **argv) {
-   double a = std::atan(1.0);                        // pi/4
-   double b = 2.0;
-   double c = 3.0;
-   int num_deriv = 2;                                // Number of independent variables
-   // Fad objects
-   Sacado::Fad::DFad<double> afad(num_deriv, 0, a); // First (0) indep. var.  Derivative array is [1.0 0.0]
-   Sacado::Fad::DFad<double> bfad(num_deriv, 1, b); // Second (1) indep. var.  Derivative array is [0.0 1.0]
-   Sacado::Fad::DFad<double> cfad(c);               // Passive variable
-   Sacado::Fad::DFad<double> rfad;                  // Result
-   double r = func(a, b, c);                        // Compute function
-   rfad = func(afad, bfad, cfad);                   // Compute function and derivative with AD
-   // Extract value and derivatives
-   double r_ad = rfad.val();     // r
-   double drda_ad = rfad.dx(0);  // dr/da
-   double drdb_ad = rfad.dx(1);  // dr/db

**sacado/example/dfad_example.cpp**

## Seed matrix example

- #include "Sacado.hpp"
- // The function to differentiate
- template <typename ScalarT>
- ScalarT func(const ScalarT& a, const ScalarT& b, const ScalarT& c) {
-   ScalarT r = c*std::log(b+1.)/std::sin(a);
-   return r;
- }
- int main(int argc, char **argv) {
-   double a = std::atan(1.0);                        // pi/4
-   double b = 2.0;
-   double c = 3.0;
-   int num_deriv = 1;                                // Number of derivative components
-   double v1 = 0.5; double v2 = 0.3;                 // Seed matrix
-   // Fad objects
-   Sacado::Fad::DFad<double> afad(num_deriv, a);  afad.fastAccessDx(0) = v1;
-   Sacado::Fad::DFad<double> bfad(num_deriv, b);  bfad.fastAccessDx(0) = v2;
-   Sacado::Fad::DFad<double> cfad(c);               // Passive variable
-   double r = func(a, b, c);                                  // Compute function
-   Sacado::Fad::DFad<double> rfad = func(afad, bfad, cfad);  // Compute function and derivative with AD
-   // Extract value and derivatives
-   double r_ad = rfad.val();     // r
-   double drda_ad = rfad.dx(0);  // dr/da*v1 + dr/da*v2

**Seed matrix example**

## sacado/example/sfad_example.cpp

- #include "Sacado.hpp"
- // The function to differentiate
- template <typename ScalarT>
- ScalarT func(const ScalarT& a, const ScalarT& b, const ScalarT& c) {
-   ScalarT r = c*std::log(b+1.)/std::sin(a);
-   return r;
- }
- int main(int argc, char **argv) {
-   double a = std::atan(1.0);                          // pi/4
-   double b = 2.0;
-   double c = 3.0;
-   int num_deriv = 2;                                  // Number of independent variables
-   // Fad objects
-   Sacado::Fad::SFad<double,2> afad(num_deriv, 0, a); // First (0) indep. var.  Derivative array is [1.0 0.0]
-   Sacado::Fad::SFad<double,2> bfad(num_deriv, 1, b); // Second (1) indep. var.  Derivative array is [0.0 1.0]
-   Sacado::Fad::SFad<double,2> cfad(c);               // Passive variable
-   Sacado::Fad::SFad<double,2> rfad;                  // Result
-   double r = func(a, b, c);                          // Compute function
-   rfad = func(afad, bfad, cfad);                     // Compute function and derivative with AD
-   // Extract value and derivatives
-   double r_ad = rfad.val();     // r
-   double drda_ad = rfad.dx(0);  // dr/da
-   double drdb_ad = rfad.dx(1);  // dr/db

**sacado/example/sfad_example.cpp**

## sacado/example/trad_example.cpp

- #include "Sacado.hpp"
- // The function to differentiate
- template <typename ScalarT>
- ScalarT func(const ScalarT& a, const ScalarT& b, const ScalarT& c) {
-   ScalarT r = c*std::log(b+1.)/std::sin(a);
-   return r;
- }
- int main(int argc, char **argv) {
-   double a = std::atan(1.0);                          // pi/4
-   double b = 2.0;
-   double c = 3.0;
-   int num_deriv = 2;                                  // Number of independent variables
-   // Rad objects
-   Sacado::Rad::ADvar<double> arad = a;
-   Sacado::Rad::ADvar<double> brad = b;
-   Sacado::Rad::ADvar<double> crad = c;              // Passive variable
-   Sacado::Rad::ADvar<double> rrad;                  // Result
-   double r = func(a, b, c);                         // Compute function
-   rrad = func(arad, brad, crad);                    // Compute function and derivative with AD
-   Sacado::Rad::ADVar<double>::Gradcomp();           // Compute gradients
-   // Extract value and derivatives
-   double r_ad = rrad.val();     // r
-   double drda_ad = arad.adj();  // dr/da
-   double drdb_ad = brad.adj();  // dr/db

**sacado/example/trad_example.cpp**

## Computing Higher Derivatives

- AD classes are templated, so AD classes can be nested to compute higher derivatives
- Forward-forward:
- Reverse-forward:
- Forward-Taylor:
- Reverse-Taylor:
- Etc...

## Forward or Reverse?

- Forward:  Computes derivatives column-wise
- Number of independent variables <= number of dependent variables
- Square Jacobians for Newton’s method
- Sensitivities with small numbers of parameters
- Algorithm naturally calls for Jacobian-vector/matrix products
- (Block) Matrix-free Newton-Krylov
- Reverse:  Computes derivatives row-wise
- Number of independent variables > number of dependent variables
- Gradients of scalar valued functions
- Sensitivities with respect to large numbers of parameters
- Algorithm naturally calls for Jacobian-transpose-vector/matrix products
- (Block) Matrix-free solves of transpose matrix
- Optimization

## Differentiating Element-Based Codes

- Global residual computation (ignoring boundary computations):
- Jacobian computation:
- Jacobian-transpose product computation:
- Hybrid symbolic/AD procedure
- Element-level derivatives computed via AD
- Exactly the same as how you would do this “manually”
- Avoids parallelization issues

**Notes:**
X is, e.g., discretized global PDE solution, f is, e.g., the global finite-element residual

## AD for Element-Level Derivatives

- Element computations are
- Narrow and shallow
- Dense
- Template application’s element code via C++ templates
- App developers code and maintain single templated code base
- Easy to add new AD types
- Not a black-box approach
- Requires very good knowledge of the application code
- Ideas developed within Charon
- Large scale finite element PDE semiconductor device simulation
- 1DFEM Example in Sacado demonstrates how to leverage AD in a finite element application code
- --enable-sacado-fem-example

**Notes:**

Given the simplicity of the element computations, they are well suited to simple applications of AD
Very object-oriented C++ code element code, given lack of source transformation tools, we wanted to investigate whether operator overloading-based would be up to the task

Developed these ideas within the context of Charon, a large-scale semiconductor device simulation code
Templated finite element residual computation in Charon, and the electric current calculation, important later in the talk

## Performance

**Performance**

- Scalability of the element-level derivative computation

- Set of N hypothetical chemical species:

- DOF per element = 4*N

- Forward mode AD
- Faster than FD
- Better scalability in number of PDEs
- Analytic derivative
- Provides Jacobian for all Charon physics
- Reverse mode AD
- Scalable adjoint/gradient

**Notes:**
Recall since derivatives are computed at the element-level, the relevant quantity for the expense of the derivative computation is the number of DOF per element, not the number of elements.

So we looked at the cost of both the Jacobian and Jacobian-tranpose-product fills in terms of both run-time and flops, relative to the cost of the residual fill for a steady-state mass-transfer problem.

Jacobian:  50% more flops for FAD, but more than 3 times faster
Adjoint:  Jacobian-transpose-product in cost of approx. 9 residual fills, independent of DOF

Conclusion:  AD is certainly efficient enough

## How to use Sacado

- Template code to be differentiated:  double -> ScalarT
- Replace independent/dependent variables with AD variables
- Initialize seed matrix
- Forward:  Derivative array of i’th independent variable is i’th row of seed matrix
- Reverse:  Derivative array of i’th dependent variable is i’th row of seed matrix
- Currently Rad only supports one column at a time
- Evaluate function on AD variables
- Instantiates template classes/functions
- Extract derivatives
- Forward:  Derivative components of dependent variables
- Reverse:  Derivative components of independent variables

## How Sacado relates to other packages

- Many Trilinos packages need derivatives
- NOX (nonlinear solves)
- LOCA (stability analysis)
- Rythmos (time integration)
- MOOCHO,  Aristos (optimization)
- Sacado does not provide these derivatives directly
- Sacado is not a black-box AD solution
- Sacado provides low level AD capabilities
- Application codes use Sacado to build derivatives these packages need

## Best Practices

- Don’t differentiate your global function with AD
- Only use AD for the hard, nonlinear parts
- Never differentiate solvers with AD...instead use AD for the derivative of the solution
- Always put Fad inside of Rad, not the other way around

## Where Sacado is at now

- Sacado is ready for use
- --enable-sacado
- Code development:
- Forward AD with memory pool
- Higher derivatives
- Extending reverse mode
- Mostly working on
- Getting ready for release
- Documentation, examples, tutorials, etc...
- Developing 1DFEM example application
- Demonstrates how to leverage Sacado in element-based codes

- Auxiliary Slides

## RAD: Specialized overloading for rf

- Reverse mode AD tool developed by David Gay
- AD data type ADvar
- Forward sweep: evaluate residual fill on Advar type
- Computes values, partials
- Store values, partials, and connectivity in “tape”
- Reverse sweep: ADcontext::Gradcomp();
- Accumulates adjoints
- Reclaims “tape” memory for future gradient computations
- Block memory allocation for efficiency.
- Templated to support higher derivatives (e.g., Hessians)

- void ADcontext::Gradcomp() {
- Derp *d = Derp::LastDerp;
- d->c->aval = 1;
- for(; d; d = d->next)
- 	 d->a->aval += d->da * d->c->aval;
- 	 // ... (arrange to recycle memory)

- Significantly less overhead than other approaches
- Computes and stores partials on forward sweep
- Only memory get/put, +, * on reverse sweep
- ADOL-C
- Stores representation of each operation (name) and computes partial on reverse sweep
- Allows reuse of tape

## RAD Performance

| Relative Times (f + rf) for f = 1 |  |
| --- | --- |
| Handcoded gradient |  1.07 |
| RAD |  9.14 |
| nlc |  1.00 |
| ADOL-C new tape | 55.0 |
| ADOL-C old tape | 15.4 |

## Templating Application Codes for AD

**template <typename ScalarT> **
**void fillResidual(ScalarT *x, ScalarT *f);**
**#include “fillResidualImpl.h”**
**fillResidualImpl.h**
**template <typename ScalarT> **
**void fillResidual(ScalarT *x, ScalarT *f) {**
**	****// Fill elemental residual**
**	****...**
**}**
**main.C**
**#include “fillResidual.h”**
**int main() {**
**	****double *x, *f;**
**  ****Fad<double> *x_fad, *f_fad;**
**	****...**
**	****// Fill residual**
**	****fillResidual<double>(x,f)**
**  ****// Fill Jacobian**
**  ****fillResidual< Fad<double> >(x_fad, f_fad);**
**	****...**
**}**

**Templating Application Codes for AD**

**fillResidual.h**
**void fillResidual(double *x, double *f);**
**fillResidual.C**
**void fillResidual(double *x, double *f) {**
**	****// Fill elemental residual**
**	****...**
**}**
**main.C**
**#include “fillResidual.h”**
**int main() {**
**	****double *x, *f;**
**	****...**
**	****// Fill residual**
**	****fillResidual(x,f)**
**	****...**
**}**

- Our interface to an application code is a templated elemental residual fill
- Templating makes it easy to interchange AD data types
- Developers only need to code residuals

**Notes:**
Our interface to an application code is a residual fill templated on the scalar type (i.e., double)
- Allows easy substitution of doubles with AD types for derivative computations

Templated code is difficult to deal with in a day-to-day setting
- Changes to the implementation of templated functions require recompiling of many source files 
- Very time consuming, especially for optimized builds
- Makes it impractical to template a large application code or framework
- Interfacing templated to non-template code requires separate code to manage each template instantiation
- Makes it difficult to add new AD types to application code build

## Templating Application Codes for AD

**template <typename ScalarT> **
**void fillResidual(ScalarT *x, ScalarT *f);**
**fillResidual.C**
**template <typename ScalarT> **
**void fillResidual(ScalarT *x, ScalarT *f) {**
**	****// Fill elemental residual**
**	****...**
**}**
**template void fillResidual<double>(...);**
**template void fillReisudal< Fad<double> >(...);**
**main.C**
**#include “fillResidual.h”**
**int main() {**
**	****double *x, *f;**
**  ****Fad<double> *x_fad, *f_fad;**
**	****...**
**	****// Fill residual**
**	****fillResidual<double>(x,f)**
**  ****// Fill Jacobian**
**  ****fillResidual< Fad<double> >(x_fad, f_fad);**
**	****...**
**}**

**Templating Application Codes for AD**

- Developed reusable “template infrastructure” components:
- Macros that explicitly instantiate template classes/functions preserving the original layout of the application code/framework into separate translation units
- Sequence containers and iterators storing all instantiations of template classes/functions providing access to instantiations, independent of the number and types of instantiations
- Adapators interfacing templated AD code to non-templated derivative code (e.g., source transformation of FORTRAN)
- Eliminates a significant obstacle to widespread use of AD at Sandia
- Easier to create/manage templated code
- Easy to add AD types for new derivative computations