source: trunk/test_more/jacobian.cpp @ 2794

Last change on this file since 2794 was 2506, checked in by bradbell, 7 years ago

Change Licenses: CPL-1.0 -> EPL-1.0, GPL-2.0->GPL-3.0

  • Property svn:keywords set to Id
File size: 2.7 KB
Line 
1/* $Id: jacobian.cpp 2506 2012-10-24 19:36:49Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13# include <cppad/cppad.hpp>
14
15namespace { // ---------------------------------------------------------
16
17template <class Vector>
18Vector eval_g(const Vector&  x)
19{       Vector g(2);
20
21        g[0] = x[0] * x[1] * x[2] * x[3];
22        g[1] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
23
24        return g;
25}
26
27template <class Vector>
28Vector eval_jac_g(const Vector&  x)
29{       Vector jac_g(8);
30
31        // g[0] = x[0] * x[1] * x[2] * x[3];
32        jac_g[0] = x[1] * x[2] * x[3];
33        jac_g[1] = x[0] * x[2] * x[3];
34        jac_g[2] = x[0] * x[1] * x[3];
35        jac_g[3] = x[0] * x[1] * x[2];
36
37        // g[1] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
38        jac_g[4+0] = 2. * x[0];
39        jac_g[4+1] = 2. * x[1];
40        jac_g[4+2] = 2. * x[2];
41        jac_g[4+3] = 2. * x[3];
42
43        return jac_g;
44}
45
46
47} // End empty namespace
48
49bool jacobian(void)
50{       bool ok = true;
51        using CppAD::vector;
52        size_t i, j, k;
53
54        size_t n = 4;
55        size_t m = 2;
56        vector< CppAD::AD<double> > ad_x(n);
57        vector< CppAD::AD<double> > ad_g(m);
58
59        vector<double> x(n);
60        x[0] = 1.; x[1] = 5.0; x[2] = 5.0; x[3] = 1.0;
61        for(j = 0; j < n; j++)
62                ad_x[j] = x[j];
63        //
64        CppAD::Independent(ad_x);
65        ad_g = eval_g(ad_x);
66        CppAD::ADFun<double> fun_g(ad_x, ad_g);
67
68        vector<double> check(m * n);
69        check = eval_jac_g(x);
70
71        // regular jacobian
72        vector<double> jac_g = fun_g.Jacobian(x);
73        for(k = 0; k < m *n; k++)
74                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
75
76        // one argument sparse jacobian
77        jac_g = fun_g.SparseJacobian(x);
78        for(k = 0; k < m *n; k++)
79                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
80
81        // sparse jacobian using bool vectors
82        CPPAD_TESTVECTOR(bool) p_b(m * n) , r_b(n * n);
83        for(i = 0; i < n; i++)
84                for(j = 0; j < n; j++)
85                        r_b[i * n + j] = (i == j);
86        p_b = fun_g.ForSparseJac(n, r_b);
87        jac_g = fun_g.SparseJacobian(x, p_b);
88        for(k = 0; k < m *n; k++)
89                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
90
91        // sparse jacobian using vectors of sets
92        std::vector< std::set<size_t> > p_s(m) , r_s(n);
93        for(i = 0; i < n; i++)
94                for(j = 0; j < n; j++)
95                        r_s[i].insert(j);
96        p_s = fun_g.ForSparseJac(n, r_s);
97        jac_g = fun_g.SparseJacobian(x, p_s);
98        for(k = 0; k < m *n; k++)
99                ok &= CppAD::NearEqual(jac_g[k], check[k], 1e-10, 1e-10);
100
101        return ok;
102}
103// END PROGRAM
Note: See TracBrowser for help on using the repository browser.