1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-18 Bradley M. Bell
3
4 CppAD is distributed under the terms of the
5 Eclipse Public License Version 2.0.
6
7 This Source Code may also be made available under the following
8 Secondary License when the conditions for such availability set forth
9 in the Eclipse Public License, Version 2.0 are satisfied:
10 GNU General Public License, Version 2.0 or later.
11 ---------------------------------------------------------------------------- */
12 # include <cppad/cppad.hpp>
13 # include <omp.h>
14
15 namespace { // BEGIN_EMPTY_NAMESPACE
16
17 using CppAD::vector;
18
19 // ----------------------------------------------------------------------------
20 // prefer reverse mode during computation of Jacobians
21
22 // example_tmb_atomic
23 class example_tmb_atomic : public CppAD::atomic_base<double> {
24 public:
25 // constructor
example_tmb_atomic(const std::string & name)26 example_tmb_atomic(const std::string& name)
27 : CppAD::atomic_base<double>(name)
28 { }
29 // forward (only implement zero order)
forward(size_t p,size_t q,const vector<bool> & vx,vector<bool> & vy,const vector<double> & tx,vector<double> & ty)30 virtual bool forward(
31 size_t p ,
32 size_t q ,
33 const vector<bool>& vx ,
34 vector<bool>& vy ,
35 const vector<double>& tx ,
36 vector<double>& ty )
37 {
38 // check for errors in usage
39 bool ok = p == 0 && q == 0;
40 ok &= tx.size() == 1;
41 ok &= ty.size() == 1;
42 ok &= vx.size() <= 1;
43 if( ! ok )
44 return false;
45
46 // variable information
47 if( vx.size() > 0 )
48 vy[0] = vx[0];
49
50 // y = 1 / x
51 ty[0] = 1.0 / tx[0];
52
53 return ok;
54 }
55 // reverse (implement first order)
reverse(size_t q,const vector<double> & tx,const vector<double> & ty,vector<double> & px,const vector<double> & py)56 virtual bool reverse(
57 size_t q ,
58 const vector<double>& tx ,
59 const vector<double>& ty ,
60 vector<double>& px ,
61 const vector<double>& py )
62 {
63 // check for errors in usage
64 bool ok = q == 0;
65 ok &= tx.size() == 1;
66 ok &= ty.size() == 1;
67 ok &= px.size() == 1;
68 ok &= py.size() == 1;
69 if( ! ok )
70 return false;
71
72 // y = 1 / x
73 // dy/dx = - 1 / (x * x)
74 double dy_dx = -1.0 / ( tx[0] * tx[0] );
75 px[0] = py[0] * dy_dx;
76
77 return ok;
78 }
79 };
80
81 } // END_EMPTY_NAMESPACE
82
prefer_reverse(void)83 bool prefer_reverse(void)
84 { bool ok = true;
85 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
86
87 // Create atomic functions
88 example_tmb_atomic afun("reciprocal");
89
90 // Declare independent variables
91 size_t n = 1;
92 CPPAD_TESTVECTOR( CppAD::AD<double> ) ax(n);
93 ax[0] = 5.0;
94 CppAD::Independent(ax);
95
96 // Compute dependent variables
97 size_t m = 1;
98 CPPAD_TESTVECTOR( CppAD::AD<double> ) ay(m);
99 afun(ax, ay);
100
101 // Create f(x) = 1 / x
102 CppAD::ADFun<double> f(ax, ay);
103
104 // Use Jacobian to compute f'(x) = - 1 / (x * x).
105 // This would fail with the normal CppAD distribution because it would use
106 // first order forward mode for the calculation.
107 CPPAD_TESTVECTOR(double) x(n), dy_dx(m);
108 x[0] = 2.0;
109 dy_dx = f.Jacobian(x);
110
111 // check the result
112 double check = -1.0 / (x[0] * x[0]);
113 ok &= CppAD::NearEqual(dy_dx[0], check, eps99, eps99);
114
115 return ok;
116 }
117