Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
47template <typename Scalar>
48Scalar 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
57template <int N>
58struct 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};
76template <>
77struct 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)
96template <typename TermIterator>
97double
98extract_derivative(const double& x, TermIterator term, TermIterator term_end)
99{
100 return x;
101}
102template <typename T, typename TermIterator>
103double
104extract_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
117int 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}
exp(expr.val())
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
const int N
int main()
Fad specializations for Teuchos::BLAS wrappers.
Uncopyable z
double extract_derivative(const double &x, TermIterator term, TermIterator term_end)
Sacado::Fad::DFad< double > FadType
Scalar func(const int m, const Scalar x[])
static type apply(const int n, const double x, const double v[])
static type apply(const int n, const int i, const double x)
static type apply(const int n, const double x, const double v[])
MakeFad< N-1 >::type nested_type
static type apply(const int n, const int i, const double x)
Sacado::mpl::apply< FadType, nested_type >::type type
static SACADO_INLINE_FUNCTION const T & eval(const T &x)
F::template apply< A1, A2, A3, A4, A5 >::type type
const double tol