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.
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$ */
4
5/* --------------------------------------------------------------------------
7
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.
14-------------------------------------------------------------------------- */
15
16/*
17$begin RevSparseJac$$18spell 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 31index RevSparseJac$$
32$index reverse, sparse Jacobian$$33index sparse, reverse Jacobian$$ 34$index pattern, reverse Jacobian$$35 36head Syntax$$
37$syntax%%r% = %F%.RevSparseJac(%p%, %s%)%$$38pre 39$$ 40$syntax%%r% = %F%.RevSparseJac(%p%, %s%, %packed%)%$$41 42 43head Purpose$$
44We use $latex F : B^n \rightarrow B^m$$to denote the 45xref/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 49latex $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%
61%$$62 63head 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 68head p$$
69The argument $italic p$$has prototype 70syntax% 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 82head s$$
83The argument $italic s$$has prototype 84syntax% 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 97head 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. 105pre 106 107$$
108The default value for $italic packed$$is true; i.e., 109the value used if it is not present. 110 111head r$$ 112The return value$italic r$$has prototype 113syntax% 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$$124latex $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 130xref/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, 136xref/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 156head Example$$
157$children% 158 example/rev_sparse_jac.cpp 159%$$160The file 161xref/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
171
172template <class Base, class VectorBool, class VectorSet>
173void ForSparseJac(
174        size_t                 total_num_var    ,
175        size_t                 p                ,
176        const VectorBool&      s                ,
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
191
193                p > 0,
194                "RevSparseJac: p (first argument) is not greater than zero"
195        );
196
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++)
210
211                for(j = 0; j < p; j++) if( s[ i * m + 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++)
227
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>
247        size_t              p      ,
248        const VectorBool&   s      ,
249        bool                packed )
251        VectorBool r( p * n );
252
253        if( packed )
255                        total_num_var_   ,
256                        p                ,
257                        s                ,
260                        play_            ,
261                        r
262                );
263        }
264        else
266                        total_num_var_   ,
267                        p                ,
268                        s                ,