CppAD: A C++ Algorithmic Differentiation Package  20130918
old_atomic.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines