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