source: trunk/test_more/sparse_vec_ad.cpp @ 3520

Last change on this file since 3520 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.9 KB
Line 
1/* $Id: sparse_vec_ad.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
14# include <cppad/cppad.hpp>
15
16bool sparse_vec_ad(void)
17{       bool ok = true;
18        using namespace CppAD;
19
20        // dimension of the domain space
21        size_t n = 3; 
22
23        size_t i, j;
24
25        // independent variable vector
26        CPPAD_TESTVECTOR(AD<double>) X(n);
27        for(j = 0; j < n; j++)
28                X[j] = AD<double>(j); 
29        Independent(X);
30
31        // dependent variable vector
32        size_t m = n;
33        CPPAD_TESTVECTOR(AD<double>) Y(m);
34
35        // check results vector
36        CPPAD_TESTVECTOR( bool )  Check(m * n);
37
38        // Create a VecAD so that there are two in the tape and the sparsity
39        // pattern depends on the second one (checks addressing VecAD objects)
40        VecAD<double> W(n);
41
42        // VecAD equal to X
43        VecAD<double> Z(n);
44        AD<double> J;
45        for(j = 0; j < n; j++)
46        {       J    = AD<double>(j);
47                W[J] = X[0];
48                Z[J] = X[j];
49
50                // y[i] depends on x[j] for j <= i
51                // (and is non-linear for j <= 1).
52                if( j == 1 )
53                        Y[j] = Z[J] * Z[J];
54                else    Y[j] = Z[J];
55        }
56
57        // compute dependent variables values
58        AD<double> P = 1;
59        J = AD<double>(0);
60        for(j = 0; j < n; j++)
61        {       for(i = 0; i < m; i++)
62                        Check[ i * m + j ] = (j <= i);
63        }
64       
65        // create function object F : X -> Y
66        ADFun<double> F(X, Y);
67
68        // dependency matrix for the identity function W(x) = x
69        CPPAD_TESTVECTOR( bool ) Identity(n * n);
70        for(i = 0; i < n; i++)
71        {       for(j = 0; j < n; j++)
72                        Identity[ i * n + j ] = false;
73                Identity[ i * n + i ] = true;
74        }
75        // evaluate the dependency matrix for Identity(F(x))
76        CPPAD_TESTVECTOR( bool ) Px(m * n);
77        Px = F.RevSparseJac(n, Identity);
78
79        // check values
80        for(i = 0; i < m; i++)
81        {       for(j = 0; j < n; j++)
82                        ok &= (Px[i * m + j] == Check[i * m + j]);
83        }       
84
85        // evaluate the dependency matrix for F(Identity(x))
86        CPPAD_TESTVECTOR( bool ) Py(m * n);
87        Py = F.ForSparseJac(n, Identity);
88
89        // check values
90        for(i = 0; i < m; i++)
91        {       for(j = 0; j < n; j++)
92                        ok &= (Py[i * m + j] == Check[i * m + j]);
93        }       
94
95        // test sparsity pattern for Hessian of F_2 ( Identity(x) )
96        CPPAD_TESTVECTOR(bool) Hy(m);
97        for(i = 0; i < m; i++)
98                Hy[i] = false;
99        Hy[2] = true;
100        CPPAD_TESTVECTOR(bool) Pxx(n * n);
101        Pxx = F.RevSparseHes(n, Hy);
102        for(i = 0; i < n; i++)
103        {       for(j = 0; j < n; j++)
104                        ok &= (Pxx[i * n + j] == false );
105        }
106
107        // test sparsity pattern for Hessian of F_1 ( Identity(x) )
108        for(i = 0; i < m; i++)
109                Hy[i] = false;
110        Hy[1] = true;
111        Pxx = F.RevSparseHes(n, Hy);
112        for(i = 0; i < n; i++)
113        {       for(j = 0; j < n; j++)
114                        ok &= (Pxx[i * n + j] == ( (i <= 1) & (j <= 1) ) );
115        }
116
117
118        return ok;
119}
Note: See TracBrowser for help on using the repository browser.