source: trunk/test_more/for_sparse_jac.cpp @ 1532

Last change on this file since 1532 was 1532, checked in by bradbell, 11 years ago

trunk: Add option for sparse Jacobians to used vector_set for calculations.

for_sparse_jac.cpp: test forward Jacobian sparsity using vector_set.
rev_sparse_jac.cpp: test reverse Jacobian sparsity using vector_set.
whats_new_09.omh: user's view of the changes.
ad_fun.hpp: add for_jac_sparse_set to private data (move private to front).
for_sparse_jac.hpp: add packed (true or false) option.
rev_sparse_jac.hpp: add packed (true or false) option.

  • Property svn:keywords set to Id
File size: 6.3 KB
Line 
1/* $Id: for_sparse_jac.cpp 1532 2009-09-28 11:55:30Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Common 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# include <vector>
16# include <valarray>
17
18# define CheckOp(Op)                   \
19        Y[index] = X[0] Op 2.;         \
20        Check[index * n + 0] = true;   \
21        Check[index * n + 1] = false;  \
22        Check[index * n + 2] = false;  \
23        index++;                       \
24        Y[index] = X[0] Op X[1];       \
25        Check[index * n + 0] = true;   \
26        Check[index * n + 1] = true;   \
27        Check[index * n + 2] = false;  \
28        index++;                       \
29        Y[index] = 3.   Op X[1];       \
30        Check[index * n + 0] = false;  \
31        Check[index * n + 1] = true;   \
32        Check[index * n + 2] = false;  \
33        index++;
34
35# define CheckUnaryFun(Fun)            \
36        Y[index] = Fun(X[0]);          \
37        Check[index * n + 0] = true;   \
38        Check[index * n + 1] = false;  \
39        Check[index * n + 2] = false;  \
40        index++;                       \
41        Y[index] = Fun(X[0] + X[1]);   \
42        Check[index * n + 0] = true;   \
43        Check[index * n + 1] = true;   \
44        Check[index * n + 2] = false;  \
45        index++;                       \
46        Y[index] = Fun(X[1]);          \
47        Check[index * n + 0] = false;  \
48        Check[index * n + 1] = true;   \
49        Check[index * n + 2] = false;  \
50        index++;
51
52
53# define CheckBinaryFun(Fun)           \
54        Y[index] = Fun( X[0] , 2.);    \
55        Check[index * n + 0] = true;   \
56        Check[index * n + 1] = false;  \
57        Check[index * n + 2] = false;  \
58        index++;                       \
59        Y[index] = Fun( X[0] , X[1]);  \
60        Check[index * n + 0] = true;   \
61        Check[index * n + 1] = true;   \
62        Check[index * n + 2] = false;  \
63        index++;                       \
64        Y[index] = Fun( 3.   , X[1]);  \
65        Check[index * n + 0] = false;  \
66        Check[index * n + 1] = true;   \
67        Check[index * n + 2] = false;  \
68        index++;
69
70
71namespace { // Begin empty namespace
72
73bool case_one(bool packed)
74{       bool ok = true;
75        using namespace CppAD;
76
77        // dimension of the domain space
78        size_t n = 3; 
79
80        // dimension of the range space
81        size_t m = (4 + 11 + 1) * 3 + 4;
82
83        // independent variable vector
84        CPPAD_TEST_VECTOR< AD<double> > X(n);
85        X[0] = .1; 
86        X[1] = .2;
87        X[2] = .3;
88        Independent(X);
89
90        // dependent variable vector
91        CPPAD_TEST_VECTOR< AD<double> > Y(m);
92
93        // check results vector
94        CPPAD_TEST_VECTOR< bool >       Check(m * n);
95
96        // initialize index into Y
97        size_t index = 0;
98
99        // 4 binary operators
100        CheckOp(+);
101        CheckOp(-);
102        CheckOp(*);
103        CheckOp(/);
104
105        // 11 unary functions
106        CheckUnaryFun(abs);
107        CheckUnaryFun(acos);
108        CheckUnaryFun(asin);
109        CheckUnaryFun(atan);
110        CheckUnaryFun(cos);
111        CheckUnaryFun(cosh);
112        CheckUnaryFun(exp);
113        CheckUnaryFun(log);
114        CheckUnaryFun(sin);
115        CheckUnaryFun(sinh);
116        CheckUnaryFun(sqrt);
117
118        // 1 binary function
119        CheckBinaryFun(pow);
120
121        // conditional expression (value of comparision does not matter)
122        Y[index] = CondExpLt(X[0], X[1], X[0], AD<double>(2.));
123        Check[index * n + 0] = true;
124        Check[index * n + 1] = false;
125        Check[index * n + 2] = false;
126        index++;
127        Y[index] = CondExpLt(X[0], X[1], X[0], X[1]);
128        Check[index * n + 0] = true;
129        Check[index * n + 1] = true;
130        Check[index * n + 2] = false;
131        index++;
132        Y[index] = CondExpLt(X[0], X[1], AD<double>(3.), X[1]);
133        Check[index * n + 0] = false;
134        Check[index * n + 1] = true;
135        Check[index * n + 2] = false;
136        index++;
137
138        // non-trival composition
139        Y[index] = X[0] * X[1] + X[1] * X[2];
140        Check[index * n + 0] = true;
141        Check[index * n + 1] = true;
142        Check[index * n + 2] = true;
143        index++;
144
145        // check final index
146        assert( index == m );
147
148        // create function object F : X -> Y
149        ADFun<double> F(X, Y);
150
151        // dependency matrix for the identity function W(x) = x
152        CPPAD_TEST_VECTOR< bool > Px(n * n);
153        size_t i, j;
154        for(i = 0; i < n; i++)
155        {       for(j = 0; j < n; j++)
156                        Px[ i * n + j ] = false;
157                Px[ i * n + i ] = true;
158        }
159
160        // evaluate the dependency matrix for F(X(x))
161        CPPAD_TEST_VECTOR< bool > Py(m * n);
162        Py = F.ForSparseJac(n, Px, packed);
163
164        // check values
165        for(i = 0; i < m; i++)
166        {       for(j = 0; j < n; j++)
167                        ok &= (Py[i * n + j] == Check[i * n + j]);
168        }       
169
170        return ok;
171}
172
173bool case_two(bool packed)
174{       bool ok = true;
175        using namespace CppAD;
176
177        // dimension of the domain space
178        size_t n = 3; 
179
180        // dimension of the range space
181        size_t m = 3;
182
183        // inialize the vector as zero
184        CppAD::VecAD<double> Z(n - 1);
185        size_t k;
186        for(k = 0; k < n-1; k++)
187                Z[k] = 0.;
188
189        // independent variable vector
190        CPPAD_TEST_VECTOR< AD<double> > X(n);
191        X[0] = 0.; 
192        X[1] = 1.;
193        X[2] = 2.;
194        Independent(X);
195
196        // VecAD vector is going to depend on X[1] and X[2]
197        Z[ X[0] ] = X[1];
198        Z[ X[1] ] = X[2]; 
199
200        // dependent variable vector
201        CPPAD_TEST_VECTOR< AD<double> > Y(m);
202
203        // check results vector
204        CPPAD_TEST_VECTOR< bool >       Check(m * n);
205
206        // initialize index into Y
207        size_t index = 0;
208
209        // First component only depends on X[0];
210        Y[index]             = X[0];
211        Check[index * n + 0] = true;
212        Check[index * n + 1] = false;
213        Check[index * n + 2] = false;
214        index++;
215
216        // Second component depends on the vector Z
217        AD<double> zero(0);
218        Y[index]             = Z[zero]; // Load by a parameter
219        Check[index * n + 0] = false;
220        Check[index * n + 1] = true;
221        Check[index * n + 2] = true;
222        index++;
223
224        // Third component depends on the vector Z
225        Y[index]             = Z[ X[0] ]; // Load by a variable
226        Check[index * n + 0] = false;
227        Check[index * n + 1] = true;
228        Check[index * n + 2] = true;
229        index++;
230
231        // check final index
232        assert( index == m );
233
234        // create function object F : X -> Y
235        ADFun<double> F(X, Y);
236
237        // dependency matrix for the identity function W(x) = x
238        CPPAD_TEST_VECTOR< bool > Px(n * n);
239        size_t i, j;
240        for(i = 0; i < n; i++)
241        {       for(j = 0; j < n; j++)
242                        Px[ i * n + j ] = false;
243                Px[ i * n + i ] = true;
244        }
245
246        // evaluate the dependency matrix for F(X(x))
247        CPPAD_TEST_VECTOR< bool > Py(m * n);
248        Py = F.ForSparseJac(n, Px, packed);
249
250        // check values
251        for(i = 0; i < m; i++)
252        {       for(j = 0; j < n; j++)
253                        ok &= (Py[i * n + j] == Check[i * n + j]);
254        }       
255
256        return ok;
257}
258
259} // End empty namespace
260
261bool for_sparse_jac(void)
262{       bool ok = true;
263
264        ok &= case_one(true);
265        ok &= case_two(true);
266
267        ok &= case_one(false);
268        ok &= case_two(false);
269
270        return ok;
271}
Note: See TracBrowser for help on using the repository browser.