source: trunk/example/sparse/colpack_hes.cpp @ 3875

Last change on this file since 3875 was 3875, checked in by bradbell, 3 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 9b561e6163a715c60405323b27831e920e887464
end hash code: c696698d38ecd5c03c12e6b87cce6e80c26386bf

commit c696698d38ecd5c03c12e6b87cce6e80c26386bf
Author: Brad Bell <bradbell@…>
Date: Mon Jan 30 10:11:28 2017 -0700

Automatic changes generated by auto-tools.

commit 6e283b7b0c690233d8f429248a0d2341f47a1301
Author: Brad Bell <bradbell@…>
Date: Mon Jan 30 10:10:16 2017 -0700

Move sparse examples to example/sparse directory.

File size: 3.7 KB
Line 
1/* --------------------------------------------------------------------------
2CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
3
4CppAD is distributed under multiple licenses. This distribution is under
5the terms of the
6                    Eclipse Public License Version 1.0.
7
8A copy of this license is included in the COPYING file of this distribution.
9Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
10-------------------------------------------------------------------------- */
11
12/*
13$begin colpack_hes.cpp$$
14$spell
15        colpack_hes
16        jacobian
17$$
18
19$section Using ColPack: Example and Test$$
20$mindex colpack hessian sparse$$
21
22
23$code
24$srcfile%example/sparse/colpack_hes.cpp%0%// BEGIN C++%// END C++%1%$$
25$$
26
27$end
28*/
29// BEGIN C++
30
31# include <cppad/cppad.hpp>
32bool colpack_hes(void)
33{       bool ok = true;
34        using CppAD::AD;
35        using CppAD::NearEqual;
36        typedef CPPAD_TESTVECTOR(AD<double>) a_vector;
37        typedef CPPAD_TESTVECTOR(double)     d_vector;
38        typedef CppAD::vector<size_t>        i_vector;
39        size_t i, j, k, ell;
40        double eps = 10. * CppAD::numeric_limits<double>::epsilon();
41
42        // domain space vector
43        size_t n = 5;
44        a_vector  a_x(n);
45        for(j = 0; j < n; j++)
46                a_x[j] = AD<double> (0);
47
48        // declare independent variables and starting recording
49        CppAD::Independent(a_x);
50
51        // colpack example case where hessian is a spear head
52        // i.e, H(i, j) non zero implies i = 0, j = 0, or i = j
53        AD<double> sum = 0.0;
54        // partial_0 partial_j = x[j]
55        // partial_j partial_j = x[0]
56        for(j = 1; j < n; j++)
57                sum += a_x[0] * a_x[j] * a_x[j] / 2.0;
58        //
59        // partial_i partial_i = 2 * x[i]
60        for(i = 0; i < n; i++)
61                sum += a_x[i] * a_x[i] * a_x[i] / 3.0;
62
63        // declare dependent variables
64        size_t m = 1;
65        a_vector  a_y(m);
66        a_y[0] = sum;
67
68        // create f: x -> y and stop tape recording
69        CppAD::ADFun<double> f(a_x, a_y);
70
71        // new value for the independent variable vector
72        d_vector x(n);
73        for(j = 0; j < n; j++)
74                x[j] = double(j + 1);
75
76        /*
77              [ 2  2  3  4  5 ]
78        hes = [ 2  5  0  0  0 ]
79              [ 3  0  7  0  0 ]
80              [ 4  0  0  9  0 ]
81              [ 5  0  0  0 11 ]
82        */
83        d_vector check(n * n);
84        for(i = 0; i < n; i++)
85        {       for(j = 0; j < n; j++)
86                {       size_t index = i * n + j;
87                        check[index] = 0.0;
88                        if( i == 0 && 1 <= j )
89                                check[index] += x[j];
90                        if( 1 <= i && j == 0 )
91                                check[index] += x[i];
92                        if( i == j )
93                        {       check[index] += 2.0 * x[i];
94                                if( i != 0 )
95                                        check[index] += x[0];
96                        }
97                }
98        }
99        // Normally one would use f.RevSparseHes to compute
100        // sparsity pattern, but for this example we extract it from check.
101        std::vector< std::set<size_t> >  p(n);
102        i_vector row, col;
103        for(i = 0; i < n; i++)
104        {       for(j = 0; j < n; j++)
105                {       ell = i * n + j;
106                        if( check[ell] != 0. )
107                        {       // insert this non-zero entry in sparsity pattern
108                                p[i].insert(j);
109
110                                // the Hessian is symmetric, so only upper lower triangle
111                                if( j <= i )
112                                {       row.push_back(i);
113                                        col.push_back(j);
114                                }
115                        }
116                }
117        }
118        size_t K = row.size();
119        d_vector hes(K);
120
121        // contrast and check results using both cppad and colpack
122        CppAD::sparse_hessian_work work;
123        for(size_t i_method = 0; i_method < 3; i_method++)
124        {       // empty work structure
125                ok &= work.color_method == "cppad.symmetric";
126                if( i_method == 2 )
127                        work.color_method = "colpack.star";
128
129                // compute Hessian
130                d_vector w(m);
131                w[0] = 1.0;
132                size_t n_sweep = f.SparseHessian(x, w, p, row, col, hes, work);
133
134                // check result
135                for(k = 0; k < K; k++)
136                {       ell = row[k] * n + col[k];
137                        ok &= NearEqual(check[ell], hes[k], eps, eps);
138                }
139                if( work.color_method != "cppad.general" )
140                        ok &= n_sweep == 2;
141                else
142                        ok &= n_sweep == 5;
143                //
144                // check that clear resets color_method to cppad.symmetric
145                work.clear();
146                ok &= work.color_method == "cppad.symmetric";
147        }
148
149        return ok;
150}
151// END C++
Note: See TracBrowser for help on using the repository browser.