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