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