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 sparse_hes.cpp$$
15 $spell
16 Cpp
17 Hessian
18 $$
19
20 $section Computing Sparse Hessian: Example and Test$$
21
22
23 $srcthisfile%0%// BEGIN C++%// END C++%1%$$
24
25 $end
26 */
27 // BEGIN C++
28 # include <cppad/cppad.hpp>
sparse_hes(void)29 bool sparse_hes(void)
30 { bool ok = true;
31 using CppAD::AD;
32 using CppAD::NearEqual;
33 //
34 typedef CPPAD_TESTVECTOR(AD<double>) a_vector;
35 typedef CPPAD_TESTVECTOR(double) d_vector;
36 typedef CPPAD_TESTVECTOR(size_t) s_vector;
37 typedef CPPAD_TESTVECTOR(bool) b_vector;
38 //
39 // domain space vector
40 size_t n = 12; // must be greater than or equal 3; see n_sweep below
41 a_vector a_x(n);
42 for(size_t j = 0; j < n; j++)
43 a_x[j] = AD<double> (0);
44 //
45 // declare independent variables and starting recording
46 CppAD::Independent(a_x);
47 //
48 // range space vector
49 size_t m = 1;
50 a_vector a_y(m);
51 a_y[0] = a_x[0] * a_x[1];
52 for(size_t j = 0; j < n; j++)
53 a_y[0] += a_x[j] * a_x[j] * a_x[j];
54 //
55 // create f: x -> y and stop tape recording
56 // (without executing zero order forward calculation)
57 CppAD::ADFun<double> f;
58 f.Dependent(a_x, a_y);
59 //
60 // new value for the independent variable vector, and weighting vector
61 d_vector w(m), x(n);
62 for(size_t j = 0; j < n; j++)
63 x[j] = double(j);
64 w[0] = 1.0;
65 //
66 // vector used to check the value of the hessian
67 d_vector check(n * n);
68 size_t ij = 0 * n + 1;
69 for(ij = 0; ij < n * n; ij++)
70 check[ij] = 0.0;
71 ij = 0 * n + 1;
72 check[ij] = 1.0;
73 ij = 1 * n + 0;
74 check[ij] = 1.0 ;
75 for(size_t j = 0; j < n; j++)
76 { ij = j * n + j;
77 check[ij] = 6.0 * x[j];
78 }
79 //
80 // compute Hessian sparsity pattern
81 b_vector select_domain(n), select_range(m);
82 for(size_t j = 0; j < n; j++)
83 select_domain[j] = true;
84 select_range[0] = true;
85 //
86 CppAD::sparse_rc<s_vector> hes_pattern;
87 bool internal_bool = false;
88 f.for_hes_sparsity(
89 select_domain, select_range, internal_bool, hes_pattern
90 );
91 //
92 // compute entire sparse Hessian (really only need lower triangle)
93 CppAD::sparse_rcv<s_vector, d_vector> subset( hes_pattern );
94 CppAD::sparse_hes_work work;
95 std::string coloring = "cppad.symmetric";
96 size_t n_sweep = f.sparse_hes(x, w, subset, hes_pattern, coloring, work);
97 ok &= n_sweep == 2;
98 //
99 const s_vector row( subset.row() );
100 const s_vector col( subset.col() );
101 const d_vector val( subset.val() );
102 size_t nnz = subset.nnz();
103 ok &= nnz == n + 2;
104 for(size_t k = 0; k < nnz; k++)
105 { ij = row[k] * n + col[k];
106 ok &= val[k] == check[ij];
107 }
108 return ok;
109 }
110 // END C++
111