1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------
# include <cppad/cppad.hpp>
namespace { // ---------------------------------------------------------
template <class Vector>
Vector eval_g(const Vector& x)
{ Vector g(2);
g[0] = x[0] * x[1] * x[2] * x[3];
g[1] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
return g;
}
template <class Vector>
Vector eval_jac_g(const Vector& x)
{ Vector jac_g(8);
// g[0] = x[0] * x[1] * x[2] * x[3];
jac_g[0] = x[1] * x[2] * x[3];
jac_g[1] = x[0] * x[2] * x[3];
jac_g[2] = x[0] * x[1] * x[3];
jac_g[3] = x[0] * x[1] * x[2];
// g[1] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
jac_g[4+0] = 2. * x[0];
jac_g[4+1] = 2. * x[1];
jac_g[4+2] = 2. * x[2];
jac_g[4+3] = 2. * x[3];
return jac_g;
}
} // End empty namespace
bool jacobian(void)
{ bool ok = true;
using CppAD::vector;
size_t i, j, k;
using CppAD::NearEqual;
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
size_t n = 4;
size_t m = 2;
vector< CppAD::AD<double> > ad_x(n);
vector< CppAD::AD<double> > ad_g(m);
vector<double> x(n);
x[0] = 1.; x[1] = 5.0; x[2] = 5.0; x[3] = 1.0;
for(j = 0; j < n; j++)
ad_x[j] = x[j];
//
CppAD::Independent(ad_x);
ad_g = eval_g(ad_x);
CppAD::ADFun<double> fun_g(ad_x, ad_g);
vector<double> check(m * n);
check = eval_jac_g(x);
// regular jacobian
vector<double> jac_g = fun_g.Jacobian(x);
for(k = 0; k < m *n; k++)
ok &= NearEqual(jac_g[k], check[k], eps99, eps99);
// one argument sparse jacobian
jac_g = fun_g.SparseJacobian(x);
for(k = 0; k < m *n; k++)
ok &= NearEqual(jac_g[k], check[k], eps99, eps99);
// sparse jacobian using bool vectors
CPPAD_TESTVECTOR(bool) p_b(m * n) , r_b(n * n);
for(i = 0; i < n; i++)
for(j = 0; j < n; j++)
r_b[i * n + j] = (i == j);
p_b = fun_g.ForSparseJac(n, r_b);
jac_g = fun_g.SparseJacobian(x, p_b);
for(k = 0; k < m *n; k++)
ok &= NearEqual(jac_g[k], check[k], eps99, eps99);
// sparse jacobian using vectors of sets
std::vector< std::set<size_t> > p_s(m) , r_s(n);
for(i = 0; i < n; i++)
for(j = 0; j < n; j++)
r_s[i].insert(j);
p_s = fun_g.ForSparseJac(n, r_s);
jac_g = fun_g.SparseJacobian(x, p_s);
for(k = 0; k < m *n; k++)
ok &= NearEqual(jac_g[k], check[k], eps99, eps99);
return ok;
}
// END PROGRAM
|
