Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
ad_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// Brief demo of Fad and Rad for computing first derivatives
31
32#include <Sacado_No_Kokkos.hpp> // for FAD and RAD
33#include <cstdio> // nicer than streams in some respects
34using std::printf;
35
36using namespace std;
37
38// Typedefs reduce gobbledygook below
39
40typedef Sacado::Fad::DFad<double> F; // FAD with # of ind. vars given later
41typedef Sacado::Fad::SFad<double,2> F2; // FAD with # of ind. vars fixed at 2
42typedef Sacado::Rad::ADvar<double> R; // for RAD
43
44template <typename T>
45const T func2(T &a, T &b) // sample function of 2 variables
46{ return sqrt(a*a + 2*b*b); }
47
48template <typename T>
49const T func(int n, T *x) // sample function of n variables
50 // == func2 when n == 2
51{
52 int i;
53 T t = 0;
54 for(i = 1; i < n; i++)
55 t += i*x[i]*x[i];
56 return sqrt(t);
57 }
58
59// demo of FAD (forward-mode AD), with general number of ind. vars
60
61 void
63{
64 F a, b, x[5], y;
65 int i, n;
66
67 printf("Fad_demo...\n\n");
68
69 // first try n == 2
70 a = 1.;
71 b = 2.;
72 // indicate the independent variables, and initialize their partials to 1:
73 a.diff(0,2); // 0 ==> this is the first independent var., of 2
74 b.diff(1,2); // 1 ==> this is the second ind. var.
75
76 y = func2(a,b);
77
78 printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
79
80 printf("partials of func2 = %g, %g\n", y.dx(0), y.dx(1));
81
82 // When we know the result was not constant (i.e., did involve ind. vars)
83 // or when hasFastAccess() is true, we access partials more quickly
84 // by using member function fastAccessDx rather than dx
85
86 if (y.hasFastAccess())
87 printf("Repeat with fastAccess: partials of func2 = %g, %g\n",
88 y.fastAccessDx(0), y.fastAccessDx(1));
89
90 // Similar exercise with general n, in this case n == 5
91 n = 5;
92 for(i = 0; i < n; i++) {
93 x[i] = i;
94 x[i].diff(i, n);
95 }
96 y = func(n, x);
97 printf("\nfunc(5,x) for x = (0,1,2,3,4) = %g\n", y.val());
98 for(i = 0; i < n; i++)
99 printf("d func / d x[%d] = %g == %g\n", i, y.dx(i), y.fastAccessDx(i));
100 }
101
102// Fad_demo2 == repeat first part of Fad_Demo with type F2 instead of F
103// i.e., with fixed-size allocations
104
105 void
107{
108 F2 a, b, y;
109
110 printf("\n\nFad2_demo...\n\n");
111
112 a = 1.;
113 b = 2.;
114 // indicate the independent variables, and initialize their partials to 1:
115 a.diff(0,2); // 0 ==> this is the first independent var., of 2
116 b.diff(1,2); // 1 ==> this is the second ind. var.
117
118 y = func2(a,b);
119
120 printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
121
122 printf("partials of func2 = %g, %g\n", y.dx(0), y.dx(1));
123
124 if (y.hasFastAccess())
125 printf("Repeat with fastAccess: partials of func2 = %g, %g\n",
126 y.fastAccessDx(0), y.fastAccessDx(1));
127 }
128
129// Fad_demo3 == repeat of Fad_Demo2 with a different constructor, one that
130// indicates the independent variables and their initial values
131// and removes the need to invoke .diff()
132
133 void
135{
136 F2 a(2,0,1.), b(2,1,2.), y;
137
138 printf("\n\nFad3_demo...\n\n");
139
140 y = func2(a,b);
141
142 printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
143
144 printf("partials of func2 = %g, %g\n", y.dx(0), y.dx(1));
145
146 if (y.hasFastAccess())
147 printf("Repeat with fastAccess: partials of func2 = %g, %g\n",
148 y.fastAccessDx(0), y.fastAccessDx(1));
149 }
150
151// Rad_demo == repeat of Fad_Demo with type R instead of F,
152// i.e., with reverse-mode rather than forward-mode AD
153
154 void
156{
157 R a, b, x[5], y;
158 int i, n;
159
160 printf("\n\nRad_demo...\n\n");
161
162 // first try n == 2
163 a = 1.;
164 b = 2.;
165
166 y = func2(a,b);
167
168 R::Gradcomp(); // do the reverse sweep
169
170 printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
171
172 printf("partials of func2 = %g, %g\n", a.adj(), b.adj());
173
174 // Similar exercise with general n, in this case n == 5
175 n = 5;
176 for(i = 0; i < n; i++)
177 x[i] = i;
178 y = func(n, x);
179 printf("\nfunc(5,x) for x = (0,1,2,3,4) = %g\n", y.val());
180
181 // the .val() values are always available; we must call Gradcomp
182 // before accessing the adjoints, i.e., the .adj() values...
183
184 R::Gradcomp();
185
186 for(i = 0; i < n; i++)
187 printf("d func / d x[%d] = %g\n", i, x[i].adj());
188 }
189
190 int
192{
193 Fad_demo();
194 Fad2_demo();
195 Fad3_demo();
196 Rad_demo();
197 return 0;
198 }
sqrt(expr.val())
#define T
Definition: Sacado_rad.hpp:573
Sacado::Fad::DFad< double > F
Definition: ad_example.cpp:40
void Fad2_demo()
Definition: ad_example.cpp:106
Sacado::Rad::ADvar< double > R
Definition: ad_example.cpp:42
void Fad_demo()
Definition: ad_example.cpp:62
const T func2(T &a, T &b)
Definition: ad_example.cpp:45
int main()
Definition: ad_example.cpp:191
void Rad_demo()
Definition: ad_example.cpp:155
const T func(int n, T *x)
Definition: ad_example.cpp:49
void Fad3_demo()
Definition: ad_example.cpp:134
Sacado::Fad::SFad< double, 2 > F2
Definition: ad_example.cpp:41
const double y