1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 Bradley M. Bell
3
4 CppAD is distributed under the terms of the
5 Eclipse Public License Version 2.0.
6
7 This Source Code may also be made available under the following
8 Secondary License when the conditions for such availability set forth
9 in the Eclipse Public License, Version 2.0 are satisfied:
10 GNU General Public License, Version 2.0 or later.
11 ---------------------------------------------------------------------------- */
12
13 /*
14 $begin subgraph_jac_rev.cpp$$
15 $spell
16 Cpp
17 Jacobian
18 $$
19
20 $section Computing Sparse Jacobian Using Reverse Mode: Example and Test$$
21
22 $srcthisfile%0%// BEGIN C++%// END C++%1%$$
23
24 $end
25 */
26 // BEGIN C++
27 # include <cppad/cppad.hpp>
subgraph_jac_rev(void)28 bool subgraph_jac_rev(void)
29 { bool ok = true;
30 //
31 using CppAD::AD;
32 using CppAD::NearEqual;
33 using CppAD::sparse_rc;
34 using CppAD::sparse_rcv;
35 //
36 typedef CPPAD_TESTVECTOR(AD<double>) a_vector;
37 typedef CPPAD_TESTVECTOR(double) d_vector;
38 typedef CPPAD_TESTVECTOR(size_t) s_vector;
39 typedef CPPAD_TESTVECTOR(bool) b_vector;
40 //
41 // domain space vector
42 size_t n = 4;
43 a_vector a_x(n);
44 for(size_t j = 0; j < n; j++)
45 a_x[j] = AD<double> (0);
46 //
47 // declare independent variables and starting recording
48 CppAD::Independent(a_x);
49 //
50 size_t m = 3;
51 a_vector a_y(m);
52 a_y[0] = a_x[0] + a_x[1];
53 a_y[1] = a_x[2] + a_x[3];
54 a_y[2] = a_x[0] + a_x[1] + a_x[2] + a_x[3] * a_x[3] / 2.;
55 //
56 // create f: x -> y and stop tape recording
57 CppAD::ADFun<double> f(a_x, a_y);
58 ok &= f.size_random() == 0;
59 //
60 // new value for the independent variable vector
61 d_vector x(n);
62 for(size_t j = 0; j < n; j++)
63 x[j] = double(j);
64 /*
65 [ 1 1 0 0 ]
66 J(x) = [ 0 0 1 1 ]
67 [ 1 1 1 x_3]
68 */
69 //
70 // row-major order values of J(x)
71 size_t nnz = 8;
72 s_vector check_row(nnz), check_col(nnz);
73 d_vector check_val(nnz);
74 for(size_t k = 0; k < nnz; k++)
75 { // check_val
76 if( k < 7 )
77 check_val[k] = 1.0;
78 else
79 check_val[k] = x[3];
80 //
81 // check_row and check_col
82 check_col[k] = k;
83 if( k < 2 )
84 check_row[k] = 0;
85 else if( k < 4 )
86 check_row[k] = 1;
87 else
88 { check_row[k] = 2;
89 check_col[k] = k - 4;
90 }
91 }
92 //
93 // select all range components of domain and range
94 b_vector select_domain(n), select_range(m);
95 for(size_t j = 0; j < n; ++j)
96 select_domain[j] = true;
97 for(size_t i = 0; i < m; ++i)
98 select_range[i] = true;
99 // -----------------------------------------------------------------------
100 // Compute Jacobian using f.subgraph_jac_rev(x, subset)
101 // -----------------------------------------------------------------------
102 //
103 // get sparsity pattern
104 bool transpose = false;
105 sparse_rc<s_vector> pattern_jac;
106 f.subgraph_sparsity(
107 select_domain, select_range, transpose, pattern_jac
108 );
109 // f.subgraph_jac_rev(x, subset)
110 sparse_rcv<s_vector, d_vector> subset( pattern_jac );
111 f.subgraph_jac_rev(x, subset);
112 //
113 // check result
114 ok &= subset.nnz() == nnz;
115 s_vector row_major = subset.row_major();
116 for(size_t k = 0; k < nnz; k++)
117 { ok &= subset.row()[ row_major[k] ] == check_row[k];
118 ok &= subset.col()[ row_major[k] ] == check_col[k];
119 ok &= subset.val()[ row_major[k] ] == check_val[k];
120 }
121 // -----------------------------------------------------------------------
122 // f.subgraph_jac_rev(select_domain, select_range, x, matrix_out)
123 // -----------------------------------------------------------------------
124 sparse_rcv<s_vector, d_vector> matrix_out;
125 f.subgraph_jac_rev(select_domain, select_range, x, matrix_out);
126 //
127 // check result
128 ok &= matrix_out.nnz() == nnz;
129 row_major = matrix_out.row_major();
130 for(size_t k = 0; k < nnz; k++)
131 { ok &= matrix_out.row()[ row_major[k] ] == check_row[k];
132 ok &= matrix_out.col()[ row_major[k] ] == check_col[k];
133 ok &= matrix_out.val()[ row_major[k] ] == check_val[k];
134 }
135 //
136 ok &= f.size_random() > 0;
137 f.clear_subgraph();
138 ok &= f.size_random() == 0;
139 return ok;
140 }
141 // END C++
142