00001 # ifndef CPPAD_REV_SPARSE_JAC_INCLUDED 00002 # define CPPAD_REV_SPARSE_JAC_INCLUDED 00003 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell 00006 00007 CppAD is distributed under multiple licenses. This distribution is under 00008 the terms of the 00009 Common Public License Version 1.0. 00010 00011 A copy of this license is included in the COPYING file of this distribution. 00012 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00013 -------------------------------------------------------------------------- */ 00014 00015 /* 00016 $begin RevSparseJac$$ 00017 $spell 00018 VecAD 00019 var 00020 Jacobian 00021 Jac 00022 const 00023 Bool 00024 Dep 00025 proportional 00026 $$ 00027 00028 $section Jacobian Sparsity Pattern: Reverse Mode$$ 00029 00030 $index RevSparseJac$$ 00031 $index reverse, sparse Jacobian$$ 00032 $index sparse, reverse Jacobian$$ 00033 $index pattern, reverse Jacobian$$ 00034 00035 $head Syntax$$ 00036 $syntax%%r% = %F%.RevSparseJac(%p%, %s%)%$$ 00037 $pre 00038 $$ 00039 00040 00041 $head Purpose$$ 00042 We use $latex F : B^n \rightarrow B^m$$ to denote the 00043 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$. 00044 For a fixed $latex p \times m$$ matrix $latex S$$, 00045 the Jacobian of $latex S * F( x )$$ 00046 with respect to $latex x$$ is 00047 $latex \[ 00048 J(x) = S * F^{(1)} ( x ) 00049 \] $$ 00050 Given a 00051 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00052 for $latex S$$, 00053 $code RevSparseJac$$ returns a sparsity pattern for the $latex J(x)$$. 00054 00055 $head f$$ 00056 The object $italic f$$ has prototype 00057 $syntax% 00058 const ADFun<%Base%> %f% 00059 %$$ 00060 00061 $head x$$ 00062 If no $xref/VecAD/$$ objects are used by the 00063 AD of $italic Base$$ 00064 $xref/glossary/Operation/Sequence/operation sequence/1/$$ 00065 stored in $italic f$$, 00066 the sparsity pattern is valid for all values $latex x \in B^n$$. 00067 $pre 00068 00069 $$ 00070 If $xref/SeqProperty/use_VecAD/f.use_VecAD/$$ is true, 00071 the sparsity patter is only valid for the value of $italic x$$ 00072 in the previous $xref/ForwardZero//zero order forward mode/$$ call 00073 $syntax% 00074 %f%.Forward(0, %x%) 00075 %$$ 00076 If there is no previous zero order forward mode call using $italic f$$, 00077 the value of the $xref/Independent//independent/$$ variables 00078 during the recording of the AD sequence of operations is used 00079 for $italic x$$. 00080 00081 $head p$$ 00082 The argument $italic p$$ has prototype 00083 $syntax% 00084 size_t %p% 00085 %$$ 00086 It specifies the number of rows in the Jacobian $latex J(x)$$. 00087 Note that the memory required for the calculation is proportional 00088 to $latex p$$ times the total number of variables 00089 in the AD operation sequence corresponding to $italic f$$ 00090 ($xref/SeqProperty/size_var/f.size_var/$$). 00091 Smaller values for $italic p$$ can be used to 00092 break the sparsity calculation 00093 into groups that do not require to much memory. 00094 00095 $head s$$ 00096 The argument $italic s$$ has prototype 00097 $syntax% 00098 const %Vector% &%s% 00099 %$$ 00100 (see $xref/RevSparseJac/Vector/Vector/$$ below) 00101 and its size is $latex p * m$$. 00102 It specifies a 00103 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00104 for the matrix $italic S$$ as follows: 00105 for $latex i = 0 , \ldots , p-1$$ and $latex j = 0 , \ldots , m-1$$. 00106 $latex \[ 00107 S_{i,j} \neq 0 ; \Rightarrow \; s [ i * m + j ] = {\rm true} 00108 \] $$ 00109 00110 $head r$$ 00111 The return value $italic r$$ has prototype 00112 $syntax% 00113 %Vector% %r% 00114 %$$ 00115 (see $xref/RevSparseJac/Vector/Vector/$$ below) 00116 and its size is $latex p * n$$. 00117 It specifies a 00118 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00119 for the matrix $latex J(x)$$ as follows: 00120 for $latex x \in B^n$$, 00121 for $latex i = 0 , \ldots , p-1$$, 00122 and $latex j = 0 , \ldots , n-1$$ 00123 $latex \[ 00124 J(x)_{i,j} \neq 0 ; \Rightarrow \; r [ i * n + j ] = {\rm true} 00125 \] $$ 00126 00127 $head Vector$$ 00128 The type $italic Vector$$ must be a $xref/SimpleVector/$$ class with 00129 $xref/SimpleVector/Elements of Specified Type/elements of type bool/$$. 00130 The routine $xref/CheckSimpleVector/$$ will generate an error message 00131 if this is not the case. 00132 In order to save memory, 00133 you may want to use a class that packs multiple elements into one 00134 storage location; for example, 00135 $xref/CppAD_vector/vectorBool/vectorBool/$$. 00136 00137 00138 $head Entire Sparsity Pattern$$ 00139 Suppose that $latex p = m$$ and 00140 $latex S$$ is the $latex m \times m$$ identity matrix, 00141 If follows that 00142 $latex \[ 00143 s [ i * q + j ] = \left\{ \begin{array}{ll} 00144 {\rm true} & {\rm if} \; i = j \\ 00145 {\rm flase} & {\rm otherwise} 00146 \end{array} \right. 00147 \] $$ 00148 is an efficient sparsity pattern for the Jacobian of $latex S$$; 00149 i.e., the choice for $italic s$$ has as few true values as possible. 00150 In this case, 00151 the corresponding value for $italic r$$ is a 00152 sparsity pattern for the Jacobian $latex J(x) = F^{(1)} ( x )$$. 00153 00154 00155 $head Example$$ 00156 $children% 00157 example/rev_sparse_jac.cpp 00158 %$$ 00159 The file 00160 $xref/RevSparseJac.cpp/$$ 00161 contains an example and test of this operation. 00162 It returns true if it succeeds and false otherwise. 00163 00164 $end 00165 ----------------------------------------------------------------------------- 00166 */ 00167 00168 // BEGIN CppAD namespace 00169 namespace CppAD { 00170 00171 template <class Base> 00172 template <class Vector> 00173 Vector ADFun<Base>::RevSparseJac(size_t p, const Vector &s) const 00174 { 00175 // type used to pack bits (must support standard bit operations) 00176 typedef size_t Pack; 00177 00178 // temporary indices 00179 size_t i, j, k, q; 00180 00181 // check Vector is Simple Vector class with bool elements 00182 CheckSimpleVector<bool, Vector>(); 00183 00184 // range and domain dimensions for F 00185 size_t m = dep_taddr.size(); 00186 size_t n = ind_taddr.size(); 00187 00188 CPPAD_ASSERT_KNOWN( 00189 p > 0, 00190 "RevSparseJac: p (first argument) is not greater than zero" 00191 ); 00192 00193 CPPAD_ASSERT_KNOWN( 00194 s.size() == p * m, 00195 "RevSparseJac: s (second argument) length is not equal to\n" 00196 "p (first argument) times range dimension for ADFun object." 00197 ); 00198 00199 00200 // number of packed values per variable on the tape 00201 size_t npv = 1 + (p - 1) / sizeof(Pack); 00202 00203 // array that will hold packed values 00204 Pack *RevJac = CPPAD_NULL; 00205 RevJac = CPPAD_TRACK_NEW_VEC(totalNumVar * npv, RevJac); 00206 00207 // update maximum memory requirement 00208 // memoryMax = std::max( memoryMax, 00209 // Memory() + totalNumVar * npv * sizeof(Pack) 00210 // ); 00211 00212 // initialize entire RevJac matrix as false 00213 for(i = 0; i < totalNumVar; i++) 00214 for(k = 0; k < npv; k++) 00215 RevJac[ i * npv + k ] = 0; 00216 00217 // set values corresponding to dependent variables 00218 Pack mask; 00219 for(j = 0; j < m; j++) 00220 { CPPAD_ASSERT_UNKNOWN( dep_taddr[j] < totalNumVar ); 00221 00222 // set bits that are true 00223 for(i = 0; i < p; i++) 00224 { k = i / sizeof(Pack); 00225 q = i - k * sizeof(Pack); 00226 mask = Pack(1) << q; 00227 if( s[ i * m + j ] ) 00228 RevJac[ dep_taddr[j] * npv + k ] |= mask; 00229 } 00230 } 00231 00232 // evaluate the sparsity patterns 00233 RevJacSweep(npv, totalNumVar, &Rec, TaylorColDim, Taylor, RevJac); 00234 00235 // return values corresponding to dependent variables 00236 Vector r(p * n); 00237 for(j = 0; j < n; j++) 00238 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < totalNumVar ); 00239 // ind_taddr[j] is operator taddr for j-th independent variable 00240 CPPAD_ASSERT_UNKNOWN( Rec.GetOp( ind_taddr[j] ) == InvOp ); 00241 00242 // set bits 00243 for(i = 0; i < p; i++) 00244 { k = i / sizeof(Pack); 00245 q = i - k * sizeof(Pack); 00246 mask = Pack(1) << q; 00247 mask &= RevJac[ ind_taddr[j] * npv + k ]; 00248 r[ i * n + j ] = (mask != 0); 00249 } 00250 } 00251 00252 // done with buffer 00253 CPPAD_TRACK_DEL_VEC(RevJac); 00254 00255 return r; 00256 } 00257 00258 } // END CppAD namespace 00259 00260 00261 # endif