CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_OLD_ATOMIC_INCLUDED 00003 # define CPPAD_OLD_ATOMIC_INCLUDED 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell 00006 00007 CppAD is distributed under multiple licenses. This distribution is under 00008 the terms of the 00009 Eclipse 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 $begin old_atomic$$ 00016 $spell 00017 hes 00018 std 00019 Jacobian 00020 jac 00021 Tvector 00022 afun 00023 vx 00024 vy 00025 bool 00026 namespace 00027 CppAD 00028 const 00029 Taylor 00030 tx 00031 ty 00032 px 00033 py 00034 $$ 00035 00036 $section User Defined Atomic AD Functions$$ 00037 $index CPPAD_USER_ATOMIC$$ 00038 $index atomic, user function$$ 00039 $index user, atomic function$$ 00040 $index operation, user atomic$$ 00041 $index function, user atomic$$ 00042 00043 $head Syntax$$ 00044 Using $code CPPAD_USER_ATOMIC$$ has been deprecated. 00045 Use $cref atomic_base$$ instead. 00046 00047 $subhead Define Function$$ 00048 $codei%CPPAD_USER_ATOMIC(%afun%, %Tvector%, %Base%, 00049 %forward%, %reverse%, %for_jac_sparse%, %rev_jac_sparse%, %rev_hes_sparse% 00050 ) 00051 %$$ 00052 00053 $subhead Use Function$$ 00054 $icode%afun%(%id%, %ax%, %ay%) 00055 %$$ 00056 00057 $subhead Callback Routines$$ 00058 $icode%ok% = %forward%(%id%, %k%, %n%, %m%, %vx%, %vy%, %tx%, %ty%) 00059 %$$ 00060 $icode%ok% = %reverse%(%id%, %k%, %n%, %m%, %tx%, %ty%, %px%, %py%) 00061 %$$ 00062 $icode%ok% = %for_jac_sparse%(%id%, %n%, %m%, %q%, %r%, %s%) 00063 %$$ 00064 $icode%ok% = %rev_jac_sparse%(%id%, %n%, %m%, %q%, %r%, %s%) 00065 %$$ 00066 $icode%ok% = %rev_hes_sparse%(%id%, %n%, %m%, %q%, %r%, %s%, %t%, %u%, %v%) 00067 %$$ 00068 00069 $subhead Free Static Memory$$ 00070 $codei%user_atomic<%Base%>::clear()%$$ 00071 00072 $head Purpose$$ 00073 In some cases, the user knows how to compute the derivative 00074 of a function 00075 $latex \[ 00076 y = f(x) \; {\rm where} \; f : B^n \rightarrow B^m 00077 \] $$ 00078 more efficiently than by coding it using $codei%AD<%Base%>%$$ 00079 $cref/atomic/glossary/Operation/Atomic/$$ operations 00080 and letting CppAD do the rest. 00081 In this case, $code CPPAD_USER_ATOMIC$$ can be used 00082 add the user code for $latex f(x)$$, and its derivatives, 00083 to the set of $codei%AD<%Base%>%$$ atomic operations. 00084 $pre 00085 00086 $$ 00087 Another possible purpose is to reduce the size of the tape; 00088 see $cref/use AD/old_atomic/Example/Use AD/$$ 00089 00090 $head Partial Implementation$$ 00091 The routines 00092 $cref/forward/old_atomic/forward/$$, 00093 $cref/reverse/old_atomic/reverse/$$, 00094 $cref/for_jac_sparse/old_atomic/for_jac_sparse/$$, 00095 $cref/rev_jac_sparse/old_atomic/rev_jac_sparse/$$, and 00096 $cref/rev_hes_sparse/old_atomic/rev_hes_sparse/$$, 00097 must be defined by the user. 00098 The $icode forward$$ the routine, 00099 for the case $icode%k% = 0%$$, must be implemented. 00100 Functions with the correct prototype, 00101 that just return $code false$$, 00102 can be used for the other cases 00103 (unless they are required by your calculations). 00104 For example, you need not implement 00105 $icode forward$$ for the case $icode%k% == 2%$$ until you require 00106 forward mode calculation of second derivatives. 00107 00108 $head CPPAD_USER_ATOMIC$$ 00109 $index CPPAD_USER_ATOMIC$$ 00110 The macro 00111 $codei% 00112 CPPAD_USER_ATOMIC(%afun%, %Tvector%, %Base%, 00113 %forward%, %reverse%, %for_jac_sparse%, %rev_jac_sparse%, %rev_hes_sparse% 00114 ) 00115 %$$ 00116 defines the $codei%AD<%Base%>%$$ routine $icode afun$$. 00117 This macro can be placed within a namespace 00118 (not the $code CppAD$$ namespace) 00119 but must be outside of any routine. 00120 00121 $subhead Tvector$$ 00122 The macro argument $icode Tvector$$ must be a 00123 $cref/simple vector template class/SimpleVector/$$. 00124 It determines the type of vectors used as arguments to the routine 00125 $icode afun$$. 00126 00127 $subhead Base$$ 00128 The macro argument $icode Base$$ specifies the 00129 $cref/base type/base_require/$$ 00130 corresponding to $codei%AD<%Base>%$$ operation sequences. 00131 Calling the routine $icode afun$$ will add the operator defined 00132 by this macro to an $codei%AD<%Base>%$$ operation sequence. 00133 00134 $head ok$$ 00135 For all routines documented below, 00136 the return value $icode ok$$ has prototype 00137 $codei% 00138 bool %ok% 00139 %$$ 00140 If it is $code true$$, the corresponding evaluation succeeded, 00141 otherwise it failed. 00142 00143 $head id$$ 00144 For all routines documented below, 00145 the argument $icode id$$ has prototype 00146 $codei% 00147 size_t %id% 00148 %$$ 00149 Its value in all other calls is the same as in the corresponding 00150 call to $icode afun$$. 00151 It can be used to store and retrieve extra information about 00152 a specific call to $icode afun$$. 00153 00154 $head k$$ 00155 For all routines documented below, the argument $icode k$$ has prototype 00156 $codei% 00157 size_t %k% 00158 %$$ 00159 The value $icode%k%$$ is the order of the Taylor coefficient that 00160 we are evaluating ($cref/forward/old_atomic/forward/$$) 00161 or taking the derivative of ($cref/reverse/old_atomic/reverse/$$). 00162 00163 $head n$$ 00164 For all routines documented below, 00165 the argument $icode n$$ has prototype 00166 $codei% 00167 size_t %n% 00168 %$$ 00169 It is the size of the vector $icode ax$$ in the corresponding call to 00170 $icode%afun%(%id%, %ax%, %ay%)%$$; i.e., 00171 the dimension of the domain space for $latex y = f(x)$$. 00172 00173 $head m$$ 00174 For all routines documented below, the argument $icode m$$ has prototype 00175 $codei% 00176 size_t %m% 00177 %$$ 00178 It is the size of the vector $icode ay$$ in the corresponding call to 00179 $icode%afun%(%id%, %ax%, %ay%)%$$; i.e., 00180 the dimension of the range space for $latex y = f(x)$$. 00181 00182 $head tx$$ 00183 For all routines documented below, 00184 the argument $icode tx$$ has prototype 00185 $codei% 00186 const CppAD::vector<%Base%>& %tx% 00187 %$$ 00188 and $icode%tx%.size() >= (%k% + 1) * %n%$$. 00189 For $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , k$$, 00190 we use the Taylor coefficient notation 00191 $latex \[ 00192 \begin{array}{rcl} 00193 x_j^\ell & = & tx [ j * ( k + 1 ) + \ell ] 00194 \\ 00195 X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^k t^k 00196 \end{array} 00197 \] $$ 00198 If $icode%tx%.size() > (%k% + 1) * %n%$$, 00199 the other components of $icode tx$$ are not specified and should not be used. 00200 Note that superscripts represent an index for $latex x_j^\ell$$ 00201 and an exponent for $latex t^\ell$$. 00202 Also note that the Taylor coefficients for $latex X(t)$$ correspond 00203 to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way: 00204 $latex \[ 00205 x_j^\ell = \frac{1}{ \ell ! } X_j^{(\ell)} (0) 00206 \] $$ 00207 00208 $head ty$$ 00209 In calls to $cref/forward/old_atomic/forward/$$, 00210 the argument $icode ty$$ has prototype 00211 $codei% 00212 CppAD::vector<%Base%>& %ty% 00213 %$$ 00214 while in calls to $cref/reverse/old_atomic/reverse/$$ it has prototype 00215 $codei% 00216 const CppAD::vector<%Base%>& %ty% 00217 %$$ 00218 For all calls, $icode%tx%.size() >= (%k% + 1) * %m%$$. 00219 For $latex i = 0 , \ldots , m-1$$ and $latex \ell = 0 , \ldots , k$$, 00220 we use the Taylor coefficient notation 00221 $latex \[ 00222 \begin{array}{rcl} 00223 y_i^\ell & = & ty [ i * ( k + 1 ) + \ell ] 00224 \\ 00225 Y_i (t) & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^k t^k + o ( t^k ) 00226 \end{array} 00227 \] $$ 00228 where $latex o( t^k ) / t^k \rightarrow 0$$ as $latex t \rightarrow 0$$. 00229 If $icode%ty%.size() > (%k% + 1) * %m%$$, 00230 the other components of $icode ty$$ are not specified and should not be used. 00231 Note that superscripts represent an index for $latex y_j^\ell$$ 00232 and an exponent for $latex t^\ell$$. 00233 Also note that the Taylor coefficients for $latex Y(t)$$ correspond 00234 to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way: 00235 $latex \[ 00236 y_j^\ell = \frac{1}{ \ell ! } Y_j^{(\ell)} (0) 00237 \] $$ 00238 00239 $subhead forward$$ 00240 In the case of $icode forward$$, 00241 for $latex i = 0 , \ldots , m-1$$, $latex ty[ i *( k + 1) + k ]$$ is an output 00242 and all the other components of $icode ty$$ are inputs. 00243 00244 $subhead reverse$$ 00245 In the case of $icode reverse$$, 00246 all the components of $icode ty$$ are inputs. 00247 00248 $head afun$$ 00249 The macro argument $icode afun$$, 00250 is the name of the AD function corresponding to this atomic 00251 operation (as it is used in the source code). 00252 CppAD uses the other functions, 00253 where the arguments are vectors with elements of type $icode Base$$, 00254 to implement the function 00255 $codei% 00256 %afun%(%id%, %ax%, %ay%) 00257 %$$ 00258 where the argument are vectors with elements of type $codei%AD<%Base%>%$$. 00259 00260 $subhead ax$$ 00261 The $icode afun$$ argument $icode ax$$ has prototype 00262 $codei% 00263 const %Tvector%< AD<%Base%> >& %ax% 00264 %$$ 00265 It is the argument vector $latex x \in B^n$$ 00266 at which the $codei%AD<%Base%>%$$ version of 00267 $latex y = f(x)$$ is to be evaluated. 00268 The dimension of the domain space for $latex y = f (x)$$ 00269 is specified by $cref/n/old_atomic/n/$$ $codei%= %ax%.size()%$$, 00270 which must be greater than zero. 00271 00272 $subhead ay$$ 00273 The $icode afun$$ result $icode ay$$ has prototype 00274 $codei% 00275 %Tvector%< AD<%Base%> >& %ay% 00276 %$$ 00277 The input values of its elements 00278 are not specified (must not matter). 00279 Upon return, it is the $codei%AD<%Base%>%$$ version of the 00280 result vector $latex y = f(x)$$. 00281 The dimension of the range space for $latex y = f (x)$$ 00282 is specified by $cref/m/old_atomic/m/$$ $codei%= %ay%.size()%$$, 00283 which must be greater than zero. 00284 00285 $subhead Parallel Mode$$ 00286 $index parallel, old_atomic$$ 00287 $index old_atomic, parallel$$ 00288 The first call to 00289 $codei% 00290 %afun%(%id%, %ax%, %ay%) 00291 %$$ 00292 must not be in $cref/parallel/ta_in_parallel/$$ mode. 00293 In addition, the 00294 $cref/old_atomic clear/old_atomic/clear/$$ 00295 routine cannot be called while in parallel mode. 00296 00297 $head forward$$ 00298 The macro argument $icode forward$$ is a 00299 user defined function 00300 $codei% 00301 %ok% = %forward%(%id%, %k%, %n%, %m%, %vx%, %vy%, %tx%, %ty%) 00302 %$$ 00303 that computes results during a $cref/forward/Forward/$$ mode sweep. 00304 For this call, we are given the Taylor coefficients in $icode tx$$ 00305 form order zero through $icode k$$, 00306 and the Taylor coefficients in $icode ty$$ with order less than $icode k$$. 00307 The $icode forward$$ routine computes the 00308 $icode k$$ order Taylor coefficients for $latex y$$ using the definition 00309 $latex Y(t) = f[ X(t) ]$$. 00310 For example, for $latex i = 0 , \ldots , m-1$$, 00311 $latex \[ 00312 \begin{array}{rcl} 00313 y_i^0 & = & Y(0) 00314 = f_i ( x^0 ) 00315 \\ 00316 y_i^1 & = & Y^{(1)} ( 0 ) 00317 = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 ) 00318 = f_i^{(1)} ( x^0 ) x^1 00319 \\ 00320 y_i^2 00321 & = & \frac{1}{2 !} Y^{(2)} (0) 00322 \\ 00323 & = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 ) 00324 + \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 ) 00325 \\ 00326 & = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1 00327 + f_i^{(1)} ( x^0 ) x^2 00328 \end{array} 00329 \] $$ 00330 Then, for $latex i = 0 , \ldots , m-1$$, it sets 00331 $latex \[ 00332 ty [ i * (k + 1) + k ] = y_i^k 00333 \] $$ 00334 The other components of $icode ty$$ must be left unchanged. 00335 00336 $subhead Usage$$ 00337 This routine is used, 00338 with $icode%vx%.size() > 0%$$ and $icode%k% == 0%$$, 00339 by calls to $icode afun$$. 00340 It is used, 00341 with $icode%vx%.size() = 0%$$ and 00342 $icode k$$ equal to the order of the derivative begin computed, 00343 by calls to $cref/forward/forward_order/$$. 00344 00345 $subhead vx$$ 00346 The $icode forward$$ argument $icode vx$$ has prototype 00347 $codei% 00348 const CppAD::vector<bool>& %vx% 00349 %$$ 00350 The case $icode%vx%.size() > 0%$$ occurs 00351 once for each call to $icode afun$$, 00352 during the call, 00353 and before any of the other callbacks corresponding to that call. 00354 Hence such a call can be used to cache information attached to 00355 the corresponding $icode id$$ 00356 (such as the elements of $icode vx$$). 00357 If $icode%vx%.size() > 0%$$ then 00358 $icode%k% == 0%$$, 00359 $icode%vx%.size() >= %n%$$, and 00360 for $latex j = 0 , \ldots , n-1$$, 00361 $icode%vx%[%j%]%$$ is true if and only if 00362 $icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$. 00363 $pre 00364 00365 $$ 00366 If $icode%vx%.size() == 0%$$, 00367 then $icode%vy%.size() == 0%$$ and neither of these vectors 00368 should be used. 00369 00370 $subhead vy$$ 00371 The $icode forward$$ argument $icode vy$$ has prototype 00372 $codei% 00373 CppAD::vector<bool>& %vy% 00374 %$$ 00375 If $icode%vy%.size() == 0%$$, it should not be used. 00376 Otherwise, 00377 $icode%k% == 0%$$ and $icode%vy%.size() >= %m%$$. 00378 The input values of the elements of $icode vy$$ 00379 are not specified (must not matter). 00380 Upon return, for $latex j = 0 , \ldots , m-1$$, 00381 $icode%vy%[%i%]%$$ is true if and only if 00382 $icode%ay%[%j%]%$$ is a variable. 00383 (CppAD uses $icode vy$$ to reduce the necessary computations.) 00384 00385 $head reverse$$ 00386 The macro argument $icode reverse$$ 00387 is a user defined function 00388 $codei% 00389 %ok% = %reverse%(%id%, %k%, %n%, %m%, %tx%, %ty%, %px%, %py%) 00390 %$$ 00391 that computes results during a $cref/reverse/Reverse/$$ mode sweep. 00392 The input value of the vectors $icode tx$$ and $icode ty$$ 00393 contain Taylor coefficient, up to order $icode k$$, 00394 for $latex X(t)$$ and $latex Y(t)$$ respectively. 00395 We use the $latex \{ x_j^\ell \}$$ and $latex \{ y_i^\ell \}$$ 00396 to denote these Taylor coefficients where the implicit range indices are 00397 $latex i = 0 , \ldots , m-1$$, 00398 $latex j = 0 , \ldots , n-1$$, 00399 $latex \ell = 0 , \ldots , k$$. 00400 Using the calculations done by $cref/forward/old_atomic/forward/$$, 00401 the Taylor coefficients $latex \{ y_i^\ell \}$$ are a function of the Taylor 00402 coefficients for $latex \{ x_j^\ell \}$$; i.e., given $latex y = f(x)$$ 00403 we define the function 00404 $latex F : B^{n \times (k+1)} \rightarrow B^{m \times (k+1)}$$ by 00405 $latex \[ 00406 y_i^\ell = F_i^\ell ( \{ x_j^\ell \} ) 00407 \] $$ 00408 We use $latex G : B^{m \times (k+1)} \rightarrow B$$ 00409 to denote an arbitrary scalar valued function of the Taylor coefficients for 00410 $latex Y(t)$$ and write $latex z = G( \{ y_i^\ell \} )$$. 00411 The $code reverse$$ routine 00412 is given the derivative of $latex z$$ with respect to 00413 $latex \{ y_i^\ell \}$$ and computes its derivative with respect 00414 to $latex \{ x_j^\ell \}$$. 00415 00416 $subhead Usage$$ 00417 This routine is used, 00418 with $icode%k% + 1%$$ equal to the order of the derivative being calculated, 00419 by calls to $cref/reverse/reverse_any/$$. 00420 00421 $subhead py$$ 00422 The $icode reverse$$ argument $icode py$$ has prototype 00423 $codei% 00424 const CppAD::vector<%Base%>& %py% 00425 %$$ 00426 and $icode%py%.size() >= (%k% + 1) * %m%$$. 00427 For $latex i = 0 , \ldots , m-1$$ and $latex \ell = 0 , \ldots , k$$, 00428 $latex \[ 00429 py[ i * (k + 1 ) + \ell ] = \partial G / \partial y_i^\ell 00430 \] $$ 00431 If $icode%py%.size() > (%k% + 1) * %m%$$, 00432 the other components of $icode py$$ are not specified and should not be used. 00433 00434 $subhead px$$ 00435 We define the function 00436 $latex \[ 00437 H ( \{ x_j^\ell \} ) = G[ F( \{ x_j^\ell \} ) ] 00438 \] $$ 00439 The $icode reverse$$ argument $icode px$$ has prototype 00440 $codei% 00441 CppAD::vector<%Base%>& %px% 00442 %$$ 00443 and $icode%px%.size() >= (%k% + 1) * %n%$$. 00444 The input values of the elements of $icode px$$ 00445 are not specified (must not matter). 00446 Upon return, 00447 for $latex j = 0 , \ldots , n-1$$ and $latex p = 0 , \ldots , k$$, 00448 $latex \[ 00449 \begin{array}{rcl} 00450 px [ j * (k + 1) + p ] & = & \partial H / \partial x_j^p 00451 \\ 00452 & = & 00453 ( \partial G / \partial \{ y_i^\ell \} ) 00454 ( \partial \{ y_i^\ell \} / \partial x_j^p ) 00455 \\ 00456 & = & 00457 \sum_{i=0}^{m-1} \sum_{\ell=0}^k 00458 ( \partial G / \partial y_i^\ell ) ( \partial y_i^\ell / \partial x_j^p ) 00459 \\ 00460 & = & 00461 \sum_{i=0}^{m-1} \sum_{\ell=p}^k 00462 py[ i * (k + 1 ) + \ell ] ( \partial F_i^\ell / \partial x_j^p ) 00463 \end{array} 00464 \] $$ 00465 Note that we have used the fact that for $latex \ell < p$$, 00466 $latex \partial F_i^\ell / \partial x_j^p = 0$$. 00467 If $icode%px%.size() > (%k% + 1) * %n%$$, 00468 the other components of $icode px$$ are not specified and should not be used. 00469 00470 $head for_jac_sparse$$ 00471 The macro argument $icode for_jac_sparse$$ 00472 is a user defined function 00473 $codei% 00474 %ok% = %for_jac_sparse%(%id%, %n%, %m%, %q%, %r%, %s%) 00475 %$$ 00476 that is used to compute results during a forward Jacobian sparsity sweep. 00477 For a fixed $latex n \times q$$ matrix $latex R$$, 00478 the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is 00479 $latex \[ 00480 S(x) = f^{(1)} (x) * R 00481 \] $$ 00482 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$, 00483 $icode for_jac_sparse$$ computes a sparsity pattern for $latex S(x)$$. 00484 00485 $subhead Usage$$ 00486 This routine is used by calls to $cref ForSparseJac$$. 00487 00488 $subhead q$$ 00489 The $icode for_jac_sparse$$ argument $icode q$$ has prototype 00490 $codei% 00491 size_t %q% 00492 %$$ 00493 It specifies the number of columns in 00494 $latex R \in B^{n \times q}$$ and the Jacobian 00495 $latex S(x) \in B^{m \times q}$$. 00496 00497 $subhead r$$ 00498 The $icode for_jac_sparse$$ argument $icode r$$ has prototype 00499 $codei% 00500 const CppAD::vector< std::set<size_t> >& %r% 00501 %$$ 00502 and $icode%r%.size() >= %n%$$. 00503 For $latex j = 0 , \ldots , n-1$$, 00504 all the elements of $icode%r%[%j%]%$$ are between 00505 zero and $icode%q%-1%$$ inclusive. 00506 This specifies a sparsity pattern for the matrix $latex R$$. 00507 00508 $subhead s$$ 00509 The $icode for_jac_sparse$$ return value $icode s$$ has prototype 00510 $codei% 00511 CppAD::vector< std::set<size_t> >& %s% 00512 %$$ 00513 and $icode%s%.size() >= %m%%$$. 00514 The input values of its sets 00515 are not specified (must not matter). Upon return 00516 for $latex i = 0 , \ldots , m-1$$, 00517 all the elements of $icode%s%[%i%]%$$ are between 00518 zero and $icode%q%-1%$$ inclusive. 00519 This represents a sparsity pattern for the matrix $latex S(x)$$. 00520 00521 $head rev_jac_sparse$$ 00522 The macro argument $icode rev_jac_sparse$$ 00523 is a user defined function 00524 $codei% 00525 %ok% = %rev_jac_sparse%(%id%, %n%, %m%, %q%, %r%, %s%) 00526 %$$ 00527 that is used to compute results during a reverse Jacobian sparsity sweep. 00528 For a fixed $latex q \times m$$ matrix $latex S$$, 00529 the Jacobian of $latex S * f( x )$$ with respect to $latex x \in B^n$$ is 00530 $latex \[ 00531 R(x) = S * f^{(1)} (x) 00532 \] $$ 00533 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex S$$, 00534 $icode rev_jac_sparse$$ computes a sparsity pattern for $latex R(x)$$. 00535 00536 $subhead Usage$$ 00537 This routine is used by calls to $cref RevSparseJac$$ 00538 and to $cref optimize$$. 00539 00540 00541 $subhead q$$ 00542 The $icode rev_jac_sparse$$ argument $icode q$$ has prototype 00543 $codei% 00544 size_t %q% 00545 %$$ 00546 It specifies the number of rows in 00547 $latex S \in B^{q \times m}$$ and the Jacobian 00548 $latex R(x) \in B^{q \times n}$$. 00549 00550 $subhead s$$ 00551 The $icode rev_jac_sparse$$ argument $icode s$$ has prototype 00552 $codei% 00553 const CppAD::vector< std::set<size_t> >& %s% 00554 %$$ 00555 and $icode%s%.size() >= %m%$$. 00556 For $latex i = 0 , \ldots , m-1$$, 00557 all the elements of $icode%s%[%i%]%$$ 00558 are between zero and $icode%q%-1%$$ inclusive. 00559 This specifies a sparsity pattern for the matrix $latex S^\R{T}$$. 00560 00561 $subhead r$$ 00562 The $icode rev_jac_sparse$$ return value $icode r$$ has prototype 00563 $codei% 00564 CppAD::vector< std::set<size_t> >& %r% 00565 %$$ 00566 and $icode%r%.size() >= %n%$$. 00567 The input values of its sets 00568 are not specified (must not matter). 00569 Upon return for $latex j = 0 , \ldots , n-1$$, 00570 all the elements of $icode%r%[%j%]%$$ 00571 are between zero and $icode%q%-1%$$ inclusive. 00572 This represents a sparsity pattern for the matrix $latex R(x)^\R{T}$$. 00573 00574 $head rev_hes_sparse$$ 00575 The macro argument $icode rev_hes_sparse$$ 00576 is a user defined function 00577 $codei% 00578 %ok% = %rev_hes_sparse%(%id%, %n%, %m%, %q%, %r%, %s%, %t%, %u%, %v%) 00579 %$$ 00580 There is an unspecified scalar valued function 00581 $latex g : B^m \rightarrow B$$. 00582 Given a sparsity pattern for $latex R$$ 00583 and information about the function $latex z = g(y)$$, 00584 this routine computes the sparsity pattern for 00585 $latex \[ 00586 V(x) = (g \circ f)^{(2)}( x ) R 00587 \] $$ 00588 00589 $subhead Usage$$ 00590 This routine is used by calls to $cref RevSparseHes$$. 00591 00592 $subhead q$$ 00593 The $icode rev_hes_sparse$$ argument $icode q$$ has prototype 00594 $codei% 00595 size_t %q% 00596 %$$ 00597 It specifies the number of columns in the sparsity patterns. 00598 00599 $subhead r$$ 00600 The $icode rev_hes_sparse$$ argument $icode r$$ has prototype 00601 $codei% 00602 const CppAD::vector< std::set<size_t> >& %r% 00603 %$$ 00604 and $icode%r%.size() >= %n%$$. 00605 For $latex j = 0 , \ldots , n-1$$, 00606 all the elements of $icode%r%[%j%]%$$ are between 00607 zero and $icode%q%-1%$$ inclusive. 00608 This specifies a sparsity pattern for the matrix $latex R \in B^{n \times q}$$. 00609 00610 $subhead s$$ 00611 The $icode rev_hes_sparse$$ argument $icode s$$ has prototype 00612 $codei% 00613 const CppAD::vector<bool>& %s% 00614 %$$ 00615 and $icode%s%.size() >= %m%$$. 00616 This specifies a sparsity pattern for the matrix 00617 $latex S(x) = g^{(1)} (y) \in B^{1 \times m}$$. 00618 00619 $subhead t$$ 00620 The $icode rev_hes_sparse$$ argument $icode t$$ has prototype 00621 $codei% 00622 CppAD::vector<bool>& %t% 00623 %$$ 00624 and $icode%t%.size() >= %n%$$. 00625 The input values of its elements 00626 are not specified (must not matter). 00627 Upon return it represents a sparsity pattern for the matrix 00628 $latex T(x) \in B^{1 \times n}$$ defined by 00629 $latex \[ 00630 T(x) = (g \circ f)^{(1)} (x) = S(x) * f^{(1)} (x) 00631 \] $$ 00632 00633 $subhead u$$ 00634 The $icode rev_hes_sparse$$ argument $icode u$$ has prototype 00635 $codei% 00636 const CppAD::vector< std::set<size_t> >& %u% 00637 %$$ 00638 and $icode%u%.size() >= %m%$$. 00639 For $latex i = 0 , \ldots , m-1$$, 00640 all the elements of $icode%u%[%i%]%$$ 00641 are between zero and $icode%q%-1%$$ inclusive. 00642 This specifies a sparsity pattern 00643 for the matrix $latex U(x) \in B^{m \times q}$$ defined by 00644 $latex \[ 00645 \begin{array}{rcl} 00646 U(x) 00647 & = & 00648 \partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{u=0} 00649 \\ 00650 & = & 00651 \partial_u \{ g^{(1)} [ y + f^{(1)} (x) R u ] \}_{u=0} 00652 \\ 00653 & = & 00654 g^{(2)} (y) f^{(1)} (x) R 00655 \end{array} 00656 \] $$ 00657 00658 $subhead v$$ 00659 The $icode rev_hes_sparse$$ argument $icode v$$ has prototype 00660 $codei% 00661 CppAD::vector< std::set<size_t> >& %v% 00662 %$$ 00663 and $icode%v%.size() >= %n%$$. 00664 The input values of its elements 00665 are not specified (must not matter). 00666 Upon return, for $latex j = 0, \ldots , n-1$$, 00667 all the elements of $icode%v%[%j%]%$$ 00668 are between zero and $icode%q%-1%$$ inclusive. 00669 This represents a sparsity pattern for the matrix 00670 $latex V(x) \in B^{n \times q}$$ defined by 00671 $latex \[ 00672 \begin{array}{rcl} 00673 V(x) 00674 & = & 00675 \partial_u [ \partial_x (g \circ f) ( x + R u ) ]_{u=0} 00676 \\ 00677 & = & 00678 \partial_u [ (g \circ f)^{(1)}( x + R u ) ]_{u=0} 00679 \\ 00680 & = & 00681 (g \circ f)^{(2)}( x ) R 00682 \\ 00683 & = & 00684 f^{(1)} (x)^\R{T} g^{(2)} ( y ) f^{(1)} (x) R 00685 + 00686 \sum_{i=1}^m [ g^{(1)} (y) ]_i \; f_i^{(2)} (x) R 00687 \\ 00688 & = & 00689 f^{(1)} (x)^\R{T} U(x) 00690 + 00691 \sum_{i=1}^m S(x)_i \; f_i^{(2)} (x) R 00692 \end{array} 00693 \] $$ 00694 00695 $head clear$$ 00696 User atomic functions hold onto static work space in order to 00697 increase speed by avoiding system memory allocation calls. 00698 The function call $codei% 00699 user_atomic<%Base%>::clear() 00700 %$$ 00701 makes to work space $cref/available/ta_available/$$ to 00702 for other uses by the same thread. 00703 This should be called when you are done using the 00704 user atomic functions for a specific value of $icode Base$$. 00705 00706 $subhead Restriction$$ 00707 The user atomic $code clear$$ routine cannot be called 00708 while in $cref/parallel/ta_in_parallel/$$ execution mode. 00709 00710 $children% 00711 example/atomic/old_reciprocal.cpp% 00712 example/atomic/old_usead_1.cpp% 00713 example/atomic/old_usead_2.cpp% 00714 example/atomic/old_tan.cpp% 00715 example/atomic/old_mat_mul.cpp 00716 %$$ 00717 $head Example$$ 00718 00719 $subhead Simple$$ 00720 The file $cref old_reciprocal.cpp$$ contains the simplest example and test 00721 of a user atomic operation. 00722 00723 $subhead Use AD$$ 00724 The examples 00725 $cref old_usead_1.cpp$$ and $cref old_usead_2.cpp$$ 00726 use AD to compute the derivatives 00727 inside a user defined atomic function. 00728 This may have the advantage of reducing the size of the tape, because 00729 a repeated section of code would only be taped once. 00730 00731 $subhead Tangent Function$$ 00732 The file $cref old_tan.cpp$$ contains an example and test 00733 implementation of the tangent function as a user atomic operation. 00734 00735 $subhead Matrix Multiplication$$ 00736 The file $cref old_mat_mul.cpp$$ contains an example and test 00737 implementation of matrix multiplication a a user atomic operation. 00738 00739 $end 00740 ------------------------------------------------------------------------------ 00741 */ 00742 # include <set> 00743 # include <cppad/local/cppad_assert.hpp> 00744 00745 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL 00746 # include <cppad/thread_alloc.hpp> 00747 00748 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00749 /*! 00750 \file old_atomic.hpp 00751 user defined atomic operations. 00752 */ 00753 00754 /*! 00755 \def CPPAD_USER_ATOMIC(afun, Tvector, 00756 forward, reverse, for_jac_sparse, rev_jac_sparse, rev_hes_sparse 00757 ) 00758 Defines the function <tt>afun(id, ax, ay)</tt> 00759 where \c id is \c ax and \c ay are vectors with <tt>AD<Base></tt> elements. 00760 00761 \par Tvector 00762 the Simple Vector template class for this function. 00763 00764 \par Base 00765 the base type for the atomic operation. 00766 00767 \par afun 00768 name of the CppAD defined function that corresponding to this operation. 00769 Note that \c afun, preceeded by a pound sign, 00770 is a version of \c afun with quotes arround it. 00771 00772 \par forward 00773 name of the user defined function that computes corresponding 00774 results during forward mode. 00775 00776 \par reverse 00777 name of the user defined function that computes corresponding 00778 results during reverse mode. 00779 00780 \par for_jac_sparse 00781 name of the user defined routine that computes corresponding 00782 results during forward mode jacobian sparsity sweeps. 00783 00784 \par rev_jac_sparse 00785 name of the user defined routine that computes corresponding 00786 results during reverse mode jacobian sparsity sweeps. 00787 00788 \par rev_hes_sparse 00789 name of the user defined routine that computes corresponding 00790 results during reverse mode Hessian sparsity sweeps. 00791 00792 \par memory allocation 00793 Note that old_atomic is used as a static object, so its objects 00794 do note get deallocated until the program terminates. 00795 */ 00796 # define CPPAD_USER_ATOMIC( \ 00797 afun , \ 00798 Tvector , \ 00799 Base , \ 00800 forward , \ 00801 reverse , \ 00802 for_jac_sparse , \ 00803 rev_jac_sparse , \ 00804 rev_hes_sparse \ 00805 ) \ 00806 inline void afun ( \ 00807 size_t id , \ 00808 const Tvector< CppAD::AD<Base> >& ax , \ 00809 Tvector< CppAD::AD<Base> >& ay \ 00810 ) \ 00811 { CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL; \ 00812 static CppAD::old_atomic<Base> fun( \ 00813 #afun , \ 00814 forward , \ 00815 reverse , \ 00816 for_jac_sparse , \ 00817 rev_jac_sparse , \ 00818 rev_hes_sparse \ 00819 ); \ 00820 fun(id, ax, ay); \ 00821 } 00822 00823 /// link so that user_atomic<Base>::clear() still works 00824 template <class Base> class user_atomic : public atomic_base<Base> { 00825 }; 00826 00827 /*! 00828 Class that actually implements the <tt>afun(id, ax, ay)</tt> calls. 00829 00830 A new old_atomic object is generated each time the user invokes 00831 the CPPAD_USER_ATOMIC macro; see static object in that macro. 00832 */ 00833 template <class Base> 00834 class old_atomic : public atomic_base<Base> { 00835 public: 00836 /// disable old_atomic<Base>::clear(void) 00837 static void clear(void) 00838 { CPPAD_ASSERT_KNOWN( 00839 false, 00840 "Depreacted API uses user_atomic<Base>::clear()" 00841 ); 00842 } 00843 /// type for user routine that computes forward mode results 00844 typedef bool (*F) ( 00845 size_t id , 00846 size_t k , 00847 size_t n , 00848 size_t m , 00849 const vector<bool>& vx , 00850 vector<bool>& vy , 00851 const vector<Base>& tx , 00852 vector<Base>& ty 00853 ); 00854 /// type for user routine that computes reverse mode results 00855 typedef bool (*R) ( 00856 size_t id , 00857 size_t k , 00858 size_t n , 00859 size_t m , 00860 const vector<Base>& tx , 00861 const vector<Base>& ty , 00862 vector<Base>& px , 00863 const vector<Base>& py 00864 ); 00865 /// type for user routine that computes forward mode Jacobian sparsity 00866 typedef bool (*FJS) ( 00867 size_t id , 00868 size_t n , 00869 size_t m , 00870 size_t q , 00871 const vector< std::set<size_t> >& r , 00872 vector< std::set<size_t> >& s 00873 ); 00874 /// type for user routine that computes reverse mode Jacobian sparsity 00875 typedef bool (*RJS) ( 00876 size_t id , 00877 size_t n , 00878 size_t m , 00879 size_t q , 00880 vector< std::set<size_t> >& r , 00881 const vector< std::set<size_t> >& s 00882 ); 00883 /// type for user routine that computes reverse mode Hessian sparsity 00884 typedef bool (*RHS) ( 00885 size_t id , 00886 size_t n , 00887 size_t m , 00888 size_t q , 00889 const vector< std::set<size_t> >& r , 00890 const vector<bool>& s , 00891 vector<bool>& t , 00892 const vector< std::set<size_t> >& u , 00893 vector< std::set<size_t> >& v 00894 ); 00895 private: 00896 /// id value corresponding to next virtual callback 00897 size_t id_; 00898 /// user's implementation of forward mode 00899 const F f_; 00900 /// user's implementation of reverse mode 00901 const R r_; 00902 /// user's implementation of forward jacobian sparsity calculations 00903 const FJS fjs_; 00904 /// user's implementation of reverse jacobian sparsity calculations 00905 const RJS rjs_; 00906 /// user's implementation of reverse Hessian sparsity calculations 00907 const RHS rhs_; 00908 00909 public: 00910 /*! 00911 Constructor called for each invocation of CPPAD_USER_ATOMIC. 00912 00913 Put this object in the list of all objects for this class and set 00914 the constant private data f_, r_, fjs_, rjs_, rhs_. 00915 00916 \param afun 00917 is the user's name for the AD version of this atomic operation. 00918 00919 \param f 00920 user routine that does forward mode calculations for this operation. 00921 00922 \param r 00923 user routine that does reverse mode calculations for this operation. 00924 00925 \param fjs 00926 user routine that does forward Jacobian sparsity calculations. 00927 00928 \param rjs 00929 user routine that does reverse Jacobian sparsity calculations. 00930 00931 \param rhs 00932 user routine that does reverse Hessian sparsity calculations. 00933 00934 \par 00935 This constructor can not be used in parallel mode because 00936 atomic_base has this restriction. 00937 */ 00938 old_atomic(const char* afun, F f, R r, FJS fjs, RJS rjs, RHS rhs) : 00939 atomic_base<Base>(afun) // name = afun 00940 , f_(f) 00941 , r_(r) 00942 , fjs_(fjs) 00943 , rjs_(rjs) 00944 , rhs_(rhs) 00945 { this->option( atomic_base<Base>::set_sparsity_enum ); 00946 } 00947 /*! 00948 Implement the user call to <tt>afun(id, ax, ay)</tt>. 00949 00950 \tparam ADVector 00951 A simple vector class with elements of type <code>AD<Base></code>. 00952 00953 \param id 00954 extra information vector that is just passed through by CppAD, 00955 and possibly used by user's routines. 00956 00957 \param ax 00958 is the argument vector for this call, 00959 <tt>ax.size()</tt> determines the number of arguments. 00960 00961 \param ay 00962 is the result vector for this call, 00963 <tt>ay.size()</tt> determines the number of results. 00964 */ 00965 template <class ADVector> 00966 void operator()(size_t id, const ADVector& ax, ADVector& ay) 00967 { // call atomic_base function object 00968 this->atomic_base<Base>::operator()(ax, ay, id); 00969 return; 00970 } 00971 /*! 00972 Store id for next virtual function callback 00973 00974 \param id 00975 id value corresponding to next virtual callback 00976 */ 00977 virtual void set_id(size_t id) 00978 { id_ = id; } 00979 /*! 00980 Link from old_atomic to forward mode 00981 00982 \copydetails atomic_base::forward 00983 */ 00984 virtual bool forward( 00985 size_t p , 00986 size_t q , 00987 const vector<bool>& vx , 00988 vector<bool>& vy , 00989 const vector<Base>& tx , 00990 vector<Base>& ty ) 00991 { CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 ); 00992 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 ); 00993 size_t n = tx.size() / (q+1); 00994 size_t m = ty.size() / (q+1); 00995 size_t i, j, k, ell; 00996 00997 vector<Base> x(n * (q+1)); 00998 vector<Base> y(m * (q+1)); 00999 vector<bool> empty; 01000 01001 // old_atomic interface can only handel one order at a time 01002 // so must just throuh hoops to get multiple orders at one time. 01003 bool ok = true; 01004 for(k = p; k <= q; k++) 01005 { for(j = 0; j < n; j++) 01006 for(ell = 0; ell <= k; ell++) 01007 x[ j * (k+1) + ell ] = tx[ j * (q+1) + ell ]; 01008 for(i = 0; i < m; i++) 01009 for(ell = 0; ell < k; ell++) 01010 y[ i * (k+1) + ell ] = ty[ i * (q+1) + ell ]; 01011 if( k == 0 ) 01012 ok &= f_(id_, k, n, m, vx, vy, x, y); 01013 else 01014 ok &= f_(id_, k, n, m, empty, empty, x, y); 01015 for(i = 0; i < m; i++) 01016 ty[ i * (q+1) + k ] = y[ i * (k+1) + k]; 01017 } 01018 return ok; 01019 } 01020 /*! 01021 Link from old_atomic to reverse mode 01022 01023 \copydetails atomic_base::reverse 01024 */ 01025 virtual bool reverse( 01026 size_t q , 01027 const vector<Base>& tx , 01028 const vector<Base>& ty , 01029 vector<Base>& px , 01030 const vector<Base>& py ) 01031 { CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 ); 01032 CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 ); 01033 size_t n = tx.size() / (q+1); 01034 size_t m = ty.size() / (q+1); 01035 bool ok = r_(id_, q, n, m, tx, ty, px, py); 01036 return ok; 01037 } 01038 /*! 01039 Link from forward Jacobian sparsity sweep to old_atomic 01040 01041 \copydetails atomic_base::for_sparse_jac 01042 */ 01043 virtual bool for_sparse_jac( 01044 size_t q , 01045 const vector< std::set<size_t> >& r , 01046 vector< std::set<size_t> >& s ) 01047 { size_t n = r.size(); 01048 size_t m = s.size(); 01049 bool ok = fjs_(id_, n, m, q, r, s); 01050 return ok; 01051 } 01052 01053 /*! 01054 Link from reverse Jacobian sparsity sweep to old_atomic. 01055 01056 \copydetails atomic_base::rev_sparse_jac 01057 */ 01058 virtual bool rev_sparse_jac( 01059 size_t q , 01060 const vector< std::set<size_t> >& rt , 01061 vector< std::set<size_t> >& st ) 01062 { size_t n = st.size(); 01063 size_t m = rt.size(); 01064 bool ok = rjs_(id_, n, m, q, st, rt); 01065 return ok; 01066 } 01067 /*! 01068 Link from reverse Hessian sparsity sweep to old_atomic 01069 01070 \copydetails atomic_base::rev_sparse_hes 01071 */ 01072 virtual bool rev_sparse_hes( 01073 const vector<bool>& vx, 01074 const vector<bool>& s , 01075 vector<bool>& t , 01076 size_t q , 01077 const vector< std::set<size_t> >& r , 01078 const vector< std::set<size_t> >& u , 01079 vector< std::set<size_t> >& v ) 01080 { size_t m = u.size(); 01081 size_t n = v.size(); 01082 CPPAD_ASSERT_UNKNOWN( r.size() == n ); 01083 CPPAD_ASSERT_UNKNOWN( s.size() == m ); 01084 CPPAD_ASSERT_UNKNOWN( t.size() == n ); 01085 // 01086 // old interface used id instead of vx 01087 bool ok = rhs_(id_, n, m, q, r, s, t, u, v); 01088 return ok; 01089 } 01090 }; 01091 01092 } // END_CPPAD_NAMESPACE 01093 # endif