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.9 KB
Line
1/* $Id: for_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 ForSparseJac$$18spell 19 var 20 Jacobian 21 Jac 22 const 23 Bool 24 proportional 25 VecAD 26$$ 27 28$section Jacobian Sparsity Pattern: Forward Mode$$29 30index ForSparseJac$$
31$index forward, sparsity Jacobian$$32index sparsity, forward Jacobian$$ 33$index pattern, forward Jacobian$$34 35head Syntax$$
36$syntax%%s% = %f%.ForSparseJac(%q%, %r%)%$$37pre 38$$ 39$syntax%%s% = %f%.ForSparseJac(%q%, %r%, %packed%)%$$40 41head Purpose$$
42We use $latex F : B^n \rightarrow B^m$$to denote the 43xref/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 51xref/glossary/Sparsity Pattern/sparsity pattern/$$ 52for$latex R$$, 53code ForSparseJac$$ returns a sparsity pattern for the $latex J(x)$$. 54 55head f$$ 56The object$italic f$$has prototype 57syntax% 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 93xref/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$$. 96latex $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 122xref/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 131head 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 141head Entire Sparsity Pattern$$
142Suppose that $latex q = n$$and 143latex R$$ is the$latex n \times n$$identity matrix, 144If follows that 145latex $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$$158children% 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
172
173template <class Base, class VectorBool, class VectorSet>
174void ForSparseJac(
175        size_t                 total_num_var    ,
176        size_t                 q                ,
177        const VectorBool&      r                ,
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
193
195                q > 0,
196                "ForSparseJac: q (first arugment) is not greater than zero"
197        );
198
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++)
213
214                // set bits that are true
215                for(j = 0; j < q; j++) if( r[ i * q + 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++)
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>
247        size_t             q      ,
248        const VectorBool&  r      ,
249        bool               packed )
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_
258                        total_num_var_   ,
259                        q                ,
260                        r                ,
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_
273                        total_num_var_   ,
274                        q                ,
275                        r                ,