source: trunk/cppad/local/for_sparse_jac.hpp @ 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: 7.9 KB
Line 
1/* $Id: for_sparse_jac.hpp 1532 2009-09-28 11:55:30Z bradbell $ */
2# ifndef CPPAD_FOR_SPARSE_JAC_INCLUDED
3# define CPPAD_FOR_SPARSE_JAC_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Common Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin ForSparseJac$$
18$spell
19        var
20        Jacobian
21        Jac
22        const
23        Bool
24        proportional
25        VecAD
26$$
27
28$section Jacobian Sparsity Pattern: Forward Mode$$
29
30$index ForSparseJac$$
31$index forward, sparsity Jacobian$$
32$index sparsity, forward Jacobian$$
33$index pattern, forward Jacobian$$
34
35$head Syntax$$
36$syntax%%s% = %f%.ForSparseJac(%q%, %r%)%$$
37$pre
38$$
39$syntax%%s% = %f%.ForSparseJac(%q%, %r%, %packed%)%$$
40
41$head Purpose$$
42We use $latex F : B^n \rightarrow B^m$$ to denote the
43$xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$.
44For a fixed $latex n \times q$$ matrix $latex R$$,
45the Jacobian of $latex F[ x + R * u ]$$
46with respect to $latex u$$ at $latex u = 0$$ is
47$latex \[
48        J(x) = F^{(1)} ( x ) * R
49\] $$
50Given a
51$xref/glossary/Sparsity Pattern/sparsity pattern/$$
52for $latex R$$,
53$code ForSparseJac$$ returns a sparsity pattern for the $latex J(x)$$.
54
55$head f$$
56The object $italic f$$ has prototype
57$syntax%
58        ADFun<%Base%> %f%
59%$$
60Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$.
61After this the sparsity pattern
62for each of the variables in the operation sequence
63is stored in the object $italic f$$.
64
65
66$head x$$
67the sparsity pattern is valid for all values of the independent
68variables in $latex x \in B^n$$
69(even if you use $cref/CondExp/$$ or $cref/VecAD/$$).
70
71$head q$$
72The argument $italic q$$ has prototype
73$syntax%
74        size_t %q%
75%$$
76It specifies the number of columns in the Jacobian $latex J(x)$$.
77Note that the memory required for the calculation is proportional
78to $latex q$$ times the total number of variables
79in the AD operation sequence corresponding to $italic f$$
80($xref/SeqProperty/size_var/f.size_var/$$).
81Smaller values for $italic q$$ can be used to
82break the sparsity calculation
83into groups that do not require to much memory.
84
85$head r$$
86The argument $italic r$$ has prototype
87$syntax%
88        const %VectorBool% &%r%
89%$$
90(see $xref/ForSparseJac/VectorBool/VectorBool/$$ below)
91and its size is $latex n * q$$.
92It specifies a
93$xref/glossary/Sparsity Pattern/sparsity pattern/$$
94for the matrix $italic R$$ as follows:
95for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , q-1$$.
96$latex \[
97        R_{i,j} \neq 0 ; \Rightarrow \; r [ i * q + j ] = {\rm true}
98\] $$
99
100$head packed$$
101If $italic packed$$ is true,
102during the sparsity calculation sets of indices are represented
103as vectors of bits that packed into words and operations are done
104on multiple bits at a time (the number of bits in a word is unspecified).
105Otherwise, sets of indices are represents using a sparse structure
106that only includes the non-zero indices and operations are done
107one index at a time.
108$pre
109
110$$
111The default value for $italic packed$$ is true; i.e.,
112the value used if it is not present.
113
114$head s$$
115The return value $italic s$$ has prototype
116$syntax%
117        %VectorBool% %s%
118%$$
119(see $xref/ForSparseJac/VectorBool/VectorBool/$$ below)
120and its size is $latex m * q$$.
121It specifies a
122$xref/glossary/Sparsity Pattern/sparsity pattern/$$
123for the matrix $latex J(x)$$ as follows:
124for $latex x \in B^n$$,
125for $latex i = 0 , \ldots , m-1$$,
126and $latex j = 0 , \ldots , q-1$$
127$latex \[
128        J(x)_{i,j} \neq 0 ; \Rightarrow \; s [ i * q + j ] = {\rm true}
129\] $$
130
131$head VectorBool$$
132The type $italic VectorBool$$ must be a $xref/SimpleVector/$$ class with
133$xref/SimpleVector/Elements of Specified Type/elements of type bool/$$.
134The routine $xref/CheckSimpleVector/$$ will generate an error message
135if this is not the case.
136In order to save memory,
137you may want to use a class that packs multiple elements into one
138storage location; for example,
139$xref/CppAD_vector/vectorBool/vectorBool/$$.
140
141$head Entire Sparsity Pattern$$
142Suppose that $latex q = n$$ and
143$latex R$$ is the $latex n \times n$$ identity matrix,
144If follows that
145$latex \[
146r [ i * q + j ] = \left\{ \begin{array}{ll}
147        {\rm true}  & {\rm if} \; i = j \\
148        {\rm flase} & {\rm otherwise}
149\end{array} \right.
150\] $$
151is an efficient sparsity pattern for $latex R$$;
152i.e., the choice for $italic r$$ has as few true values as possible.
153In this case,
154the corresponding value for $italic s$$ is a
155sparsity pattern for the Jacobian $latex J(x) = F^{(1)} ( x )$$.
156
157$head Example$$
158$children%
159        example/for_sparse_jac.cpp
160%$$
161The file
162$xref/ForSparseJac.cpp/$$
163contains an example and test of this operation.
164It returns true if it succeeds and false otherwise.
165
166$end
167-----------------------------------------------------------------------------
168*/
169
170// BEGIN CppAD namespace
171namespace CppAD {
172
173template <class Base, class VectorBool, class VectorSet> 
174void ForSparseJac(
175        size_t                 total_num_var    ,
176        size_t                 q                , 
177        const VectorBool&      r                ,
178        CppAD::vector<size_t>& dep_taddr        ,
179        CppAD::vector<size_t>& ind_taddr        ,
180        CppAD::player<Base>&   play             ,
181        VectorBool&            s                ,
182        VectorSet&             for_jac_sparsity )
183{
184        // temporary indices
185        size_t i, j;
186
187        // check VectorBool is Simple Vector class with bool elements
188        CheckSimpleVector<bool, VectorBool>();
189
190        // range and domain dimensions for F
191        size_t m = dep_taddr.size();
192        size_t n = ind_taddr.size();
193
194        CPPAD_ASSERT_KNOWN(
195                q > 0,
196                "ForSparseJac: q (first arugment) is not greater than zero"
197        );
198
199        CPPAD_ASSERT_KNOWN(
200                r.size() == n * q,
201                "ForSparseJac: r (second argument) length is not equal to\n"
202                "q (first argument) times domain dimension for ADFun object."
203        );
204
205        // allocate memory for the requested sparsity calculation
206        for_jac_sparsity.resize(total_num_var, q);
207
208        // set values corresponding to independent variables
209        for(i = 0; i < n; i++)
210        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[i] < total_num_var );
211                // ind_taddr[i] is operator taddr for i-th independent variable
212                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[i] ) == InvOp );
213
214                // set bits that are true
215                for(j = 0; j < q; j++) if( r[ i * q + j ] )
216                        for_jac_sparsity.add_element( ind_taddr[i], j);
217        }
218
219        // evaluate the sparsity patterns
220        ForJacSweep(
221                n,
222                total_num_var,
223                &play,
224                for_jac_sparsity
225        );
226
227        // return values corresponding to dependent variables
228        CPPAD_ASSERT_UNKNOWN( s.size() == m * q );
229        for(i = 0; i < m; i++)
230        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
231
232                // set bits
233                for(j = 0; j < q; j++)
234                        s[ i * q + j ] = false;
235                CPPAD_ASSERT_UNKNOWN( for_jac_sparsity.limit() == q );
236                j = for_jac_sparsity.next_element( dep_taddr[i] );
237                while( j < q )
238                {       s[i * q + j ] = true;
239                        j = for_jac_sparsity.next_element( dep_taddr[i] );
240                }
241        }
242}
243
244template <class Base>
245template <class VectorBool>
246VectorBool ADFun<Base>::ForSparseJac(
247        size_t             q      , 
248        const VectorBool&  r      ,
249        bool               packed )
250{       size_t m = dep_taddr_.size();
251        VectorBool s( m * q );
252
253        if( packed )
254        {       // free any results stored in for_jac_sparse_set_       
255                for_jac_sparse_set_.resize(0, 0);
256                // store results in for_jac_sparsity_
257                CppAD::ForSparseJac(
258                        total_num_var_   ,
259                        q                , 
260                        r                ,
261                        dep_taddr_       ,
262                        ind_taddr_       ,
263                        play_            ,
264                        s                ,
265                        for_jac_sparsity_
266                );
267        }
268        else
269        {       // free any results stored in for_jac_sparsity_
270                for_jac_sparsity_.resize(0, 0);
271                // store results in for_jac_sparse_set_
272                CppAD::ForSparseJac(
273                        total_num_var_   ,
274                        q                , 
275                        r                ,
276                        dep_taddr_       ,
277                        ind_taddr_       ,
278                        play_            ,
279                        s                ,
280                        for_jac_sparse_set_
281                );
282        }
283        return s;
284}
285
286} // END CppAD namespace
287
288# endif
Note: See TracBrowser for help on using the repository browser.