1 | /* $Id: div_zero_one.cpp 1370 2009-05-31 05:31:50Z bradbell $ */ |
---|
2 | /* -------------------------------------------------------------------------- |
---|
3 | CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell |
---|
4 | |
---|
5 | CppAD is distributed under multiple licenses. This distribution is under |
---|
6 | the terms of the |
---|
7 | Common Public License Version 1.0. |
---|
8 | |
---|
9 | A copy of this license is included in the COPYING file of this distribution. |
---|
10 | Please visit http://www.coin-or.org/CppAD/ for information on other licenses. |
---|
11 | -------------------------------------------------------------------------- */ |
---|
12 | |
---|
13 | /* |
---|
14 | Test the use of the special parameters zero and one with the multiply operator |
---|
15 | */ |
---|
16 | |
---|
17 | # include <cppad/cppad.hpp> |
---|
18 | |
---|
19 | typedef CppAD::AD<double> ADdouble; |
---|
20 | typedef CppAD::AD< ADdouble > ADDdouble; |
---|
21 | |
---|
22 | bool DivZeroOne(void) |
---|
23 | { |
---|
24 | using namespace CppAD; |
---|
25 | |
---|
26 | bool ok = true; |
---|
27 | |
---|
28 | size_t i; |
---|
29 | for(i = 0; i < 3; i++) |
---|
30 | { // run through the cases x = 0, 1, 2 |
---|
31 | |
---|
32 | size_t j; |
---|
33 | for(j = 0; j < 3; j++) |
---|
34 | { // run through the cases y = 0, 1, 2 |
---|
35 | |
---|
36 | CPPAD_TEST_VECTOR< ADdouble > x(1); |
---|
37 | x[0] = double(i); |
---|
38 | Independent(x); |
---|
39 | |
---|
40 | CPPAD_TEST_VECTOR< ADDdouble > y(1); |
---|
41 | y[0] = ADDdouble(j); |
---|
42 | Independent(y); |
---|
43 | |
---|
44 | CPPAD_TEST_VECTOR< ADDdouble > z(2); |
---|
45 | if( j == 0 ) |
---|
46 | z[0] = ADDdouble(0); |
---|
47 | else z[0] = x[0] / y[0]; |
---|
48 | if( i == 0 ) |
---|
49 | z[1] = ADDdouble(0); |
---|
50 | else |
---|
51 | { z[1] = y[0] / x[0]; |
---|
52 | z[1] /= x[0]; |
---|
53 | } |
---|
54 | |
---|
55 | // f(y) = z = { x / y , y / (x * x) } |
---|
56 | ADFun< ADdouble > f(y, z); |
---|
57 | CPPAD_TEST_VECTOR< ADdouble > u( f.Domain() ); |
---|
58 | CPPAD_TEST_VECTOR< ADdouble > v( f.Range() ); |
---|
59 | |
---|
60 | // v = f'(y) |
---|
61 | u[0] = ADdouble(1.); |
---|
62 | v = f.Forward(1, u); |
---|
63 | |
---|
64 | // check derivatives of f |
---|
65 | ADdouble check = - double(i) / double(j * j); |
---|
66 | if( j != 0 ) ok &= NearEqual( |
---|
67 | v[0], check, 1e-10, 1e-10 ); |
---|
68 | |
---|
69 | check = 1. / double(i * i); |
---|
70 | if( i != 0 ) ok &= NearEqual( |
---|
71 | v[1], check, 1e-10, 1e-10); |
---|
72 | |
---|
73 | // g(x) = f'(y) = {-x/y^2 , 1/(x * x)} |
---|
74 | ADFun<double> g(x, v); |
---|
75 | CPPAD_TEST_VECTOR< double > a( g.Domain() ); |
---|
76 | CPPAD_TEST_VECTOR< double > b( g.Range() ); |
---|
77 | |
---|
78 | // b = g'(x) |
---|
79 | a[0] = 1.; |
---|
80 | b = g.Forward(1, a); |
---|
81 | |
---|
82 | // check derivatives of g |
---|
83 | if( j != 0 ) ok &= NearEqual( |
---|
84 | b[0], - 1./double(j*j), 1e-10, 1e-10 ); |
---|
85 | if( i != 0 ) ok &= NearEqual( |
---|
86 | b[1], -2./double(i*i*i), 1e-10, 1e-10); |
---|
87 | |
---|
88 | } |
---|
89 | } |
---|
90 | |
---|
91 | return ok; |
---|
92 | } |
---|