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 code_gen_fun_sparse_jac_as_fun.cpp$$
15 $spell
16 jacobian
17 $$
18
19 $section Pass Sparse Jacobian as Code Gen Function: Example and Test$$
20
21 $srcthisfile%0%// BEGIN C++%// END C++%1%$$
22
23 $end
24 */
25 // BEGIN C++
26 # include <cppad/example/code_gen_fun.hpp>
27
sparse_jac_as_fun(void)28 bool sparse_jac_as_fun(void)
29 { bool ok = true;
30 //
31 typedef CppAD::cg::CG<double> c_double;
32 typedef CppAD::AD<c_double> ac_double;
33 //
34 typedef CppAD::vector<size_t> s_vector;
35 typedef CppAD::vector<double> d_vector;
36 typedef CppAD::vector<ac_double> ac_vector;
37 //
38 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
39
40 // domain space vector
41 size_t n = 2;
42 ac_vector ac_x(n);
43 for(size_t j = 0; j < n; ++j)
44 ac_x[j] = 1.0 / double(j + 1);
45
46 // declare independent variables and start tape recording
47 CppAD::Independent(ac_x);
48
49 // range space vector
50 size_t m = 3;
51 ac_vector ac_y(m);
52 for(size_t i = 0; i < m; ++i)
53 ac_y[i] = double(i + 1) * sin( ac_x[i % n] );
54
55 // create f: x -> y and stop tape recording
56 CppAD::ADFun<c_double> c_f(ac_x, ac_y);
57
58 // create a version of f that evalutes using ac_double
59 CppAD::ADFun<ac_double, c_double> ac_f = c_f.base2ad();
60
61 // Independent varialbes while evaluating Jacobian
62 CppAD::Independent(ac_x);
63
64 // ----------------------------------------------------------------
65 // Sparse Jacobian evaluation using any CppAD method
66 // (there are lots of choices here)
67 //
68 // pattern_eye (pattern for identity matrix)
69 CppAD::sparse_rc<s_vector> pattern_eye(n, n, n);
70 for(size_t k = 0; k < n; ++k)
71 pattern_eye.set(k, k, k);
72 //
73 // pattern_jac
74 bool transpose = false;
75 bool dependency = false;
76 bool internal_bool = true;
77 CppAD::sparse_rc<s_vector> pattern_jac;
78 // note that c_f and ac_f have the same sparsity pattern
79 c_f.for_jac_sparsity(
80 pattern_eye, transpose, dependency, internal_bool, pattern_jac
81 );
82 //
83 // ac_Jrcv
84 CppAD::sparse_rcv<s_vector, ac_vector> ac_Jrcv( pattern_jac );
85 CppAD::sparse_jac_work work;
86 std::string coloring = "cppad";
87 size_t group_max = n;
88 ac_f.sparse_jac_for(
89 group_max, ac_x, ac_Jrcv, pattern_jac, coloring, work
90 );
91 //
92 // create g: x -> non-zero elements of Jacobian
93 CppAD::ADFun<c_double> c_g(ac_x, ac_Jrcv.val());
94
95 // create compiled version of c_g
96 std::string file_name = "example_lib";
97 code_gen_fun g(file_name, c_g);
98
99 // evaluate the compiled jacobian
100 d_vector x(n);
101 for(size_t j = 0; j < n; ++j)
102 x[j] = 1.0 / double(j + 2);
103 d_vector val = g(x);
104
105 // check Jaociban values
106 size_t nnz = pattern_jac.nnz();
107 const s_vector& row(pattern_jac.row());
108 const s_vector& col(pattern_jac.col());
109 s_vector row_major = pattern_jac.row_major();
110 //
111 size_t k = 0;
112 ok &= val.size() == nnz;
113 for(size_t i = 0; i < m; ++i)
114 { for(size_t j = 0; j < n; ++j)
115 { if( j == i % n )
116 { size_t ell = row_major[k];
117 double check = double(i + 1) * cos( x[i % n] );
118 ok &= i == row[ell];
119 ok &= j == col[ell];
120 ok &= CppAD::NearEqual(val[k], check, eps99, eps99);
121 ++k;
122 }
123 }
124 }
125 ok &= k == nnz;
126 return ok;
127 }
128 // END C++
129