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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
// 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 sparse_jac_for.cpp}
Computing Sparse Jacobian Using Forward Mode: Example and Test
##############################################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end sparse_jac_for.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp>
bool sparse_jac_for(void)
{ bool ok = true;
//
using CppAD::AD;
using CppAD::NearEqual;
using CppAD::sparse_rc;
using CppAD::sparse_rcv;
//
typedef CPPAD_TESTVECTOR(AD<double>) a_vector;
typedef CPPAD_TESTVECTOR(double) d_vector;
typedef CPPAD_TESTVECTOR(size_t) s_vector;
//
// domain space vector
size_t n = 3;
a_vector a_x(n);
for(size_t j = 0; j < n; j++)
a_x[j] = AD<double> (0);
//
// declare independent variables and starting recording
CppAD::Independent(a_x);
//
size_t m = 4;
a_vector a_y(m);
a_y[0] = a_x[0] + a_x[2];
a_y[1] = a_x[0] + a_x[2];
a_y[2] = a_x[1] + a_x[2];
a_y[3] = a_x[1] + a_x[2] * a_x[2] / 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(size_t j = 0; j < n; j++)
x[j] = double(j);
/*
[ 1 0 1 ]
J(x) = [ 1 0 1 ]
[ 0 1 1 ]
[ 0 1 x_2 ]
*/
d_vector check(m * n);
//
// column-major order values of J(x)
size_t nnz = 8;
s_vector check_row(nnz), check_col(nnz);
d_vector check_val(nnz);
for(size_t k = 0; k < nnz; k++)
{ // check_val
if( k < 7 )
check_val[k] = 1.0;
else
check_val[k] = x[2];
//
// check_row and check_col
check_row[k] = k;
if( k < 2 )
check_col[k] = 0;
else if( k < 4 )
check_col[k] = 1;
else
{ check_col[k] = 2;
check_row[k] = k - 4;
}
}
//
// n by n identity matrix sparsity
sparse_rc<s_vector> pattern_in;
pattern_in.resize(n, n, n);
for(size_t k = 0; k < n; k++)
pattern_in.set(k, k, k);
//
// sparsity for J(x)
bool transpose = false;
bool dependency = false;
bool internal_bool = true;
sparse_rc<s_vector> pattern_jac;
f.for_jac_sparsity(
pattern_in, transpose, dependency, internal_bool, pattern_jac
);
//
// compute entire forward mode Jacobian
sparse_rcv<s_vector, d_vector> subset( pattern_jac );
CppAD::sparse_jac_work work;
std::string coloring = "cppad";
size_t group_max = 10;
size_t n_color = f.sparse_jac_for(
group_max, x, subset, pattern_jac, coloring, work
);
ok &= n_color == 2;
//
const s_vector row( subset.row() );
const s_vector col( subset.col() );
const d_vector val( subset.val() );
s_vector col_major = subset.col_major();
ok &= subset.nnz() == nnz;
for(size_t k = 0; k < nnz; k++)
{ ok &= row[ col_major[k] ] == check_row[k];
ok &= col[ col_major[k] ] == check_col[k];
ok &= val[ col_major[k] ] == check_val[k];
}
// compute non-zero in row 3 only
sparse_rc<s_vector> pattern_row3;
pattern_row3.resize(m, n, 2); // nr = m, nc = n, nnz = 2
pattern_row3.set(0, 3, 1); // row[0] = 3, col[0] = 1
pattern_row3.set(1, 3, 2); // row[1] = 3, col[1] = 2
sparse_rcv<s_vector, d_vector> subset_row3( pattern_row3 );
work.clear();
n_color = f.sparse_jac_for(
group_max, x, subset_row3, pattern_jac, coloring, work
);
ok &= n_color == 2;
//
const s_vector row_row3( subset_row3.row() );
const s_vector col_row3( subset_row3.col() );
const d_vector val_row3( subset_row3.val() );
ok &= subset_row3.nnz() == 2;
//
ok &= row_row3[0] == 3;
ok &= col_row3[0] == 1;
ok &= val_row3[0] == 1.0;
//
ok &= row_row3[1] == 3;
ok &= col_row3[1] == 2;
ok &= val_row3[1] == x[2];
//
return ok;
}
// END C++
|
