source: trunk/example/sparse/colpack_jac.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.0 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_jac.cpp$$
14$spell
15        colpack_jac
16        jacobian
17$$
18
19$section Using ColPack: Example and Test$$
20$mindex colpack jacobian sparse$$
21
22
23$code
24$srcfile%example/sparse/colpack_jac.cpp%0%// BEGIN C++%// END C++%1%$$
25$$
26
27$end
28*/
29// BEGIN C++
30
31# include <cppad/cppad.hpp>
32bool colpack_jac(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 = 4;
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        size_t m = 3;
52        a_vector  a_y(m);
53        a_y[0] = a_x[0] + a_x[1];
54        a_y[1] = a_x[2] + a_x[3];
55        a_y[2] = a_x[0] + a_x[1] + a_x[2] + a_x[3] * a_x[3] / 2.;
56
57        // create f: x -> y and stop tape recording
58        CppAD::ADFun<double> f(a_x, a_y);
59
60        // new value for the independent variable vector
61        d_vector x(n);
62        for(j = 0; j < n; j++)
63                x[j] = double(j);
64
65        /*
66              [ 1 1 0 0  ]
67        jac = [ 0 0 1 1  ]
68              [ 1 1 1 x_3]
69        */
70        d_vector check(m * n);
71        check[0] = 1.; check[1] = 1.; check[2]  = 0.; check[3]  = 0.;
72        check[4] = 0.; check[5] = 0.; check[6]  = 1.; check[7]  = 1.;
73        check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3];
74
75        // Normally one would use f.ForSparseJac or f.RevSparseJac to compute
76        // sparsity pattern, but for this example we extract it from check.
77        std::vector< std::set<size_t> >  p(m);
78
79        // using row and column indices to compute non-zero in rows 1 and 2
80        i_vector row, col;
81        for(i = 0; i < m; i++)
82        {       for(j = 0; j < n; j++)
83                {       ell = i * n + j;
84                        if( check[ell] != 0. )
85                        {       row.push_back(i);
86                                col.push_back(j);
87                                p[i].insert(j);
88                        }
89                }
90        }
91        size_t K = row.size();
92        d_vector jac(K);
93
94        // empty work structure
95        CppAD::sparse_jacobian_work work;
96        ok &= work.color_method == "cppad";
97
98        // choose to use ColPack
99        work.color_method = "colpack";
100
101        // forward mode
102        size_t n_sweep = f.SparseJacobianForward(x, p, row, col, jac, work);
103        for(k = 0; k < K; k++)
104        {       ell = row[k] * n + col[k];
105                ok &= NearEqual(check[ell], jac[k], eps, eps);
106        }
107        ok &= n_sweep == 4;
108
109        // reverse mode
110        work.clear();
111        work.color_method = "colpack";
112        n_sweep = f.SparseJacobianReverse(x, p, row, col, jac, work);
113        for(k = 0; k < K; k++)
114        {       ell = row[k] * n + col[k];
115                ok &= NearEqual(check[ell], jac[k], eps, eps);
116        }
117        ok &= n_sweep == 2;
118
119        return ok;
120}
121// END C++
Note: See TracBrowser for help on using the repository browser.