Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
high_order_example.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Sacado Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25 // (etphipp@sandia.gov).
26 //
27 // ***********************************************************************
28 // @HEADER
29 
30 // high_order_example
31 //
32 // usage:
33 // high_order_example
34 //
35 // output:
36 // computes a high-order derivative of a simple function along a small set
37 // of directions/modes using DFad, and then computes the same thing through
38 // the multivariate chain rule.
39 
40 #include <iostream>
41 
42 #include "Sacado.hpp"
43 
45 
46 // Function to compute derivatives of
47 template <typename Scalar>
48 Scalar func(const int m, const Scalar x[]) {
49  using std::exp;
50  Scalar z = 1.0;
51  for (int i=0; i<m; ++i)
52  z *= x[i];
53  return exp(z);
54 }
55 
56 // Build nested Fad type for computing derivatives up to order N
57 template <int N>
58 struct MakeFad {
59  // Replace double in FadType with MakeFad<N-1>::type
60  typedef typename MakeFad<N-1>::type nested_type;
62 
63  // Initialize Fad for computation of full derivative
64  static type apply(const int n, const int i, const double x) {
65  return type(n,i,MakeFad<N-1>::apply(n,i,x));
66  }
67 
68  // Initialize Fad object for derivative-matrix product
69  static type apply(const int n, const double x, const double v[]) {
70  type x_fad(n,MakeFad<N-1>::apply(n,x,v));
71  for (int i=0; i<n; ++i)
72  x_fad.fastAccessDx(i) = v[i];
73  return x_fad;
74  }
75 };
76 template <>
77 struct MakeFad<1> {
78  typedef FadType type;
79 
80  // Initialize Fad for computation of full derivative
81  static type apply(const int n, const int i, const double x) {
82  return type(n,i,x);
83  }
84 
85  // Initialize Fad object for derivative-matrix product
86  static type apply(const int n, const double x, const double v[]) {
87  type x_fad(n,x);
88  for (int i=0; i<n; ++i)
89  x_fad.fastAccessDx(i) = v[i];
90  return x_fad;
91  }
92 };
93 
94 // Extract d^r/(dx_k_1 ... dx_k_r) where the sequence k_1 ... k_r is indicated
95 // by the passed iterator (e.g., iterator into a std::vector)
96 template <typename TermIterator>
97 double
98 extract_derivative(const double& x, TermIterator term, TermIterator term_end)
99 {
100  return x;
101 }
102 template <typename T, typename TermIterator>
103 double
104 extract_derivative(const T& x, TermIterator term, TermIterator term_end)
105 {
106  // If there are no terms, return value
107  if (term == term_end)
109 
110  // Get the first term
111  auto k = *term;
112 
113  // Extract remaining terms
114  return extract_derivative(x.fastAccessDx(k), ++term, term_end);
115 }
116 
117 int main(int argc, char **argv)
118 {
119  const int deg = 6; // Order of derivative to compute
120  const int m = 3; // Number of coordinates
121  const int n = 2; // Number of modes
122  const double x0[m] = { 1.0, 1.0, 1.0 }; // Expansion point
123  const double V[m][n] = { {0.1, 0.2},
124  {0.3, 0.4},
125  {0.5, 0.6} }; // Mode coordinates
126 
127  // Derivative term we wish to extract. This corresponds to the derivative
128  // d^6/(ds_0 ds_0 ds_0 ds_1 ds_1 ds_1) = d^6/(ds_0^3 ds_1^3)
129  std::vector<int> term = { 0, 0, 0, 1, 1, 1 };
130 
131  // For y = f(x(s)), x(s) = x0 + V*s, compute
132  // d^r y / (ds_i_1 ... ds_i_r) where (i_1,...,i_r) is indicated by term
133  typedef typename MakeFad<deg>::type NestedFadType;
134  NestedFadType x_fad[m];
135  for (int i=0; i<m; ++i)
136  x_fad[i] = MakeFad<deg>::apply(n,x0[i],V[i]);
137  const NestedFadType f_fad = func(m,x_fad);
138  const double z1 = extract_derivative(f_fad, term.begin(), term.end());
139 
140  // Now compute the same thing by first computing
141  // d^k y / (dx_j_1 .. dx_j_k) for 0 <= k <= deg and j_1,...,j_k = 0,...,m-1
142  NestedFadType x_fad2[m];
143  for (int i=0; i<m; ++i)
144  x_fad2[i] = MakeFad<deg>::apply(m,i,x0[i]);
145  const NestedFadType f_fad2 = func(m,x_fad2);
146 
147  // Now compute multivariate chain-rule:
148  // d^r y / (ds_i_1 ... ds_i_r) =
149  // \sum_{j_1,...,j_r=0}^{m-1} d^r y/(dx_j_1 .. dx_j_r) *
150  // V[j_1][i_1]...V[j_r][i_r].
151  // This requires iterating (j_1,...,j_r) over the m^r-dimensional tensor
152  // product space [0,m-1] x ... x [0,m-1]
153  double z2 = 0.0;
154  const int r = term.size();
155  std::vector<int> t(r, 0);
156  bool finished = false;
157  while (!finished) {
158  const double z = extract_derivative(f_fad2, t.begin(), t.end());
159  double c = 1.0;
160  for (int i=0; i<r; ++i) {
161  c *= V[t[i]][term[i]];
162  }
163  z2 += z*c;
164 
165  // Increment t to next term in tensor product space
166  ++t[0];
167  int j=0;
168  while (j<r-1 && t[j] >= m) {
169  t[j] = 0;
170  ++j;
171  ++t[j];
172  }
173  if (t[r-1] == m)
174  finished = true;
175  }
176 
177  const double error = std::abs(z1-z2)/std::abs(z2);
178  std::cout << "z (reduced) = " << z1 << " z (full) = " << z2
179  << " error = " << error << std::endl;
180 
181  const double tol = 1.0e-14;
182  if (error < tol) {
183  std::cout << "\nExample passed!" << std::endl;
184  return 0;
185  }
186 
187  std::cout <<"\nSomething is wrong, example failed!" << std::endl;
188  return 1;
189 }
static type apply(const int n, const double x, const double v[])
abs(expr.val())
static type apply(const int n, const int i, const double x)
MakeFad< N-1 >::type nested_type
Sacado::Fad::DFad< double > FadType
F::template apply< A1, A2, A3, A4, A5 >::type type
static type apply(const int n, const double x, const double v[])
#define T
Definition: Sacado_rad.hpp:573
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
int main()
Definition: ad_example.cpp:191
std::string error
Definition: fad_expr.cpp:48
double extract_derivative(const double &x, TermIterator term, TermIterator term_end)
const double tol
const T func(int n, T *x)
Definition: ad_example.cpp:49
Sacado::mpl::apply< FadType, nested_type >::type type
static type apply(const int n, const int i, const double x)
exp(expr.val())
int n
static KOKKOS_INLINE_FUNCTION const T & eval(const T &x)