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
99
100
101
102
103
104
105
106
107
108
109
110
111
|
// 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
// ----------------------------------------------------------------------------
/*
{xrst_begin colpack_jacobian.cpp}
ColPack: Sparse Jacobian Example and Test
#########################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end colpack_jacobian.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp>
bool colpack_jacobian(void)
{ bool ok = true;
using CppAD::AD;
using CppAD::NearEqual;
typedef CPPAD_TESTVECTOR(AD<double>) a_vector;
typedef CPPAD_TESTVECTOR(double) d_vector;
typedef CppAD::vector<size_t> i_vector;
size_t i, j, k, ell;
double eps = 10. * CppAD::numeric_limits<double>::epsilon();
// domain space vector
size_t n = 4;
a_vector a_x(n);
for(j = 0; j < n; j++)
a_x[j] = AD<double> (0);
// declare independent variables and starting recording
CppAD::Independent(a_x);
size_t m = 3;
a_vector a_y(m);
a_y[0] = a_x[0] + a_x[1];
a_y[1] = a_x[2] + a_x[3];
a_y[2] = a_x[0] + a_x[1] + a_x[2] + a_x[3] * a_x[3] / 2.;
// create f: x -> y and stop tape recording
CppAD::ADFun<double> f(a_x, a_y);
// new value for the independent variable vector
d_vector x(n);
for(j = 0; j < n; j++)
x[j] = double(j);
/*
[ 1 1 0 0 ]
jac = [ 0 0 1 1 ]
[ 1 1 1 x_3]
*/
d_vector check(m * n);
check[0] = 1.; check[1] = 1.; check[2] = 0.; check[3] = 0.;
check[4] = 0.; check[5] = 0.; check[6] = 1.; check[7] = 1.;
check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3];
// Normally one would use f.ForSparseJac or f.RevSparseJac to compute
// sparsity pattern, but for this example we extract it from check.
std::vector< std::set<size_t> > p(m);
// using row and column indices to compute non-zero in rows 1 and 2
i_vector row, col;
for(i = 0; i < m; i++)
{ for(j = 0; j < n; j++)
{ ell = i * n + j;
if( check[ell] != 0. )
{ row.push_back(i);
col.push_back(j);
p[i].insert(j);
}
}
}
size_t K = row.size();
d_vector jac(K);
// empty work structure
CppAD::sparse_jacobian_work work;
ok &= work.color_method == "cppad";
// choose to use ColPack
work.color_method = "colpack";
// forward mode
size_t n_sweep = f.SparseJacobianForward(x, p, row, col, jac, work);
for(k = 0; k < K; k++)
{ ell = row[k] * n + col[k];
ok &= NearEqual(check[ell], jac[k], eps, eps);
}
ok &= n_sweep == 4;
// reverse mode
work.clear();
work.color_method = "colpack";
n_sweep = f.SparseJacobianReverse(x, p, row, col, jac, work);
for(k = 0; k < K; k++)
{ ell = row[k] * n + col[k];
ok &= NearEqual(check[ell], jac[k], eps, eps);
}
ok &= n_sweep == 2;
return ok;
}
// END C++
|
