source: trunk/cppad/local/for_sparse_jac.hpp @ 1531

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

trunk: Change BaseVector? to VectorBase?, BoolVector? to VectorBool?.

makefile.in: update corresponding version of makefile.am.

  • Property svn:keywords set to Id
File size: 6.1 KB
Line 
1/* $Id: for_sparse_jac.hpp 1531 2009-09-25 16:23:54Z 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
38
39$head Purpose$$
40We use $latex F : B^n \rightarrow B^m$$ to denote the
41$xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$.
42For a fixed $latex n \times q$$ matrix $latex R$$,
43the Jacobian of $latex F[ x + R * u ]$$
44with respect to $latex u$$ at $latex u = 0$$ is
45$latex \[
46        J(x) = F^{(1)} ( x ) * R
47\] $$
48Given a
49$xref/glossary/Sparsity Pattern/sparsity pattern/$$
50for $latex R$$,
51$code ForSparseJac$$ returns a sparsity pattern for the $latex J(x)$$.
52
53$head f$$
54The object $italic f$$ has prototype
55$syntax%
56        ADFun<%Base%> %f%
57%$$
58Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$.
59After this the sparsity pattern
60for each of the variables in the operation sequence
61is stored in the object $italic f$$.
62
63
64$head x$$
65the sparsity pattern is valid for all values of the independent
66variables in $latex x \in B^n$$
67(even if you use $cref/CondExp/$$ or $cref/VecAD/$$).
68
69$head q$$
70The argument $italic q$$ has prototype
71$syntax%
72        size_t %q%
73%$$
74It specifies the number of columns in the Jacobian $latex J(x)$$.
75Note that the memory required for the calculation is proportional
76to $latex q$$ times the total number of variables
77in the AD operation sequence corresponding to $italic f$$
78($xref/SeqProperty/size_var/f.size_var/$$).
79Smaller values for $italic q$$ can be used to
80break the sparsity calculation
81into groups that do not require to much memory.
82
83$head r$$
84The argument $italic r$$ has prototype
85$syntax%
86        const %VectorBool% &%r%
87%$$
88(see $xref/ForSparseJac/VectorBool/VectorBool/$$ below)
89and its size is $latex n * q$$.
90It specifies a
91$xref/glossary/Sparsity Pattern/sparsity pattern/$$
92for the matrix $italic R$$ as follows:
93for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , q-1$$.
94$latex \[
95        R_{i,j} \neq 0 ; \Rightarrow \; r [ i * q + j ] = {\rm true}
96\] $$
97
98$head s$$
99The return value $italic s$$ has prototype
100$syntax%
101        %VectorBool% %s%
102%$$
103(see $xref/ForSparseJac/VectorBool/VectorBool/$$ below)
104and its size is $latex m * q$$.
105It specifies a
106$xref/glossary/Sparsity Pattern/sparsity pattern/$$
107for the matrix $latex J(x)$$ as follows:
108for $latex x \in B^n$$,
109for $latex i = 0 , \ldots , m-1$$,
110and $latex j = 0 , \ldots , q-1$$
111$latex \[
112        J(x)_{i,j} \neq 0 ; \Rightarrow \; s [ i * q + j ] = {\rm true}
113\] $$
114
115$head VectorBool$$
116The type $italic VectorBool$$ must be a $xref/SimpleVector/$$ class with
117$xref/SimpleVector/Elements of Specified Type/elements of type bool/$$.
118The routine $xref/CheckSimpleVector/$$ will generate an error message
119if this is not the case.
120In order to save memory,
121you may want to use a class that packs multiple elements into one
122storage location; for example,
123$xref/CppAD_vector/vectorBool/vectorBool/$$.
124
125$head Entire Sparsity Pattern$$
126Suppose that $latex q = n$$ and
127$latex R$$ is the $latex n \times n$$ identity matrix,
128If follows that
129$latex \[
130r [ i * q + j ] = \left\{ \begin{array}{ll}
131        {\rm true}  & {\rm if} \; i = j \\
132        {\rm flase} & {\rm otherwise}
133\end{array} \right.
134\] $$
135is an efficient sparsity pattern for $latex R$$;
136i.e., the choice for $italic r$$ has as few true values as possible.
137In this case,
138the corresponding value for $italic s$$ is a
139sparsity pattern for the Jacobian $latex J(x) = F^{(1)} ( x )$$.
140
141$head Example$$
142$children%
143        example/for_sparse_jac.cpp
144%$$
145The file
146$xref/ForSparseJac.cpp/$$
147contains an example and test of this operation.
148It returns true if it succeeds and false otherwise.
149
150$end
151-----------------------------------------------------------------------------
152*/
153
154// BEGIN CppAD namespace
155namespace CppAD {
156
157template <class Base>
158template <class VectorBool>
159VectorBool ADFun<Base>::ForSparseJac(size_t q, const VectorBool &r)
160{
161        // temporary indices
162        size_t i, j;
163
164        // check VectorBool is Simple Vector class with bool elements
165        CheckSimpleVector<bool, VectorBool>();
166
167        // range and domain dimensions for F
168        size_t m = dep_taddr_.size();
169        size_t n = ind_taddr_.size();
170
171        CPPAD_ASSERT_KNOWN(
172                q > 0,
173                "ForSparseJac: q (first arugment) is not greater than zero"
174        );
175
176        CPPAD_ASSERT_KNOWN(
177                r.size() == n * q,
178                "ForSparseJac: r (second argument) length is not equal to\n"
179                "q (first argument) times domain dimension for ADFun object."
180        );
181
182        // allocate memory for the requested sparsity calculation
183        for_jac_sparsity_.resize(total_num_var_, q);
184
185        // set values corresponding to independent variables
186        for(i = 0; i < n; i++)
187        {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[i] < total_num_var_ );
188                // ind_taddr_[i] is operator taddr for i-th independent variable
189                CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[i] ) == InvOp );
190
191                // set bits that are true
192                for(j = 0; j < q; j++) if( r[ i * q + j ] )
193                        for_jac_sparsity_.add_element( ind_taddr_[i], j);
194        }
195
196        // evaluate the sparsity patterns
197        ForJacSweep(
198                n,
199                total_num_var_,
200                &play_,
201                for_jac_sparsity_
202        );
203
204        // return values corresponding to dependent variables
205        VectorBool s(m * q);
206        for(i = 0; i < m; i++)
207        {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
208
209                // set bits
210                for(j = 0; j < q; j++)
211                        s[ i * q + j ] = false;
212                CPPAD_ASSERT_UNKNOWN( for_jac_sparsity_.limit() == q );
213                j = for_jac_sparsity_.next_element( dep_taddr_[i] );
214                while( j < q )
215                {       s[i * q + j ] = true;
216                        j = for_jac_sparsity_.next_element( dep_taddr_[i] );
217                }
218        }
219
220        return s;
221}
222
223} // END CppAD namespace
224
225# endif
Note: See TracBrowser for help on using the repository browser.