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
|
// 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 subgraph_jac_rev.cpp}
Computing Sparse Jacobian Using Reverse Mode: Example and Test
##############################################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end subgraph_jac_rev.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp>
bool subgraph_jac_rev(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;
typedef CPPAD_TESTVECTOR(bool) b_vector;
//
// domain space vector
size_t n = 4;
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 = 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);
ok &= f.size_random() == 0;
//
// new value for the independent variable vector
d_vector x(n);
for(size_t j = 0; j < n; j++)
x[j] = double(j);
/*
[ 1 1 0 0 ]
J(x) = [ 0 0 1 1 ]
[ 1 1 1 x_3]
*/
//
// row-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[3];
//
// check_row and check_col
check_col[k] = k;
if( k < 2 )
check_row[k] = 0;
else if( k < 4 )
check_row[k] = 1;
else
{ check_row[k] = 2;
check_col[k] = k - 4;
}
}
//
// select all range components of domain and range
b_vector select_domain(n), select_range(m);
for(size_t j = 0; j < n; ++j)
select_domain[j] = true;
for(size_t i = 0; i < m; ++i)
select_range[i] = true;
// -----------------------------------------------------------------------
// Compute Jacobian using f.subgraph_jac_rev(x, subset)
// -----------------------------------------------------------------------
//
// get sparsity pattern
bool transpose = false;
sparse_rc<s_vector> pattern_jac;
f.subgraph_sparsity(
select_domain, select_range, transpose, pattern_jac
);
// f.subgraph_jac_rev(x, subset)
sparse_rcv<s_vector, d_vector> subset( pattern_jac );
f.subgraph_jac_rev(x, subset);
//
// check result
ok &= subset.nnz() == nnz;
s_vector row_major = subset.row_major();
for(size_t k = 0; k < nnz; k++)
{ ok &= subset.row()[ row_major[k] ] == check_row[k];
ok &= subset.col()[ row_major[k] ] == check_col[k];
ok &= subset.val()[ row_major[k] ] == check_val[k];
}
// -----------------------------------------------------------------------
// f.subgraph_jac_rev(select_domain, select_range, x, matrix_out)
// -----------------------------------------------------------------------
sparse_rcv<s_vector, d_vector> matrix_out;
f.subgraph_jac_rev(select_domain, select_range, x, matrix_out);
//
// check result
ok &= matrix_out.nnz() == nnz;
row_major = matrix_out.row_major();
for(size_t k = 0; k < nnz; k++)
{ ok &= matrix_out.row()[ row_major[k] ] == check_row[k];
ok &= matrix_out.col()[ row_major[k] ] == check_col[k];
ok &= matrix_out.val()[ row_major[k] ] == check_val[k];
}
//
ok &= f.size_random() > 0;
f.clear_subgraph();
ok &= f.size_random() == 0;
return ok;
}
// END C++
|
