CppAD: A C++ Algorithmic Differentiation Package  20130102
user_atomic.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_USER_ATOMIC_INCLUDED
00003 # define CPPAD_USER_ATOMIC_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin user_atomic$$
00018 $spell
00019      hes
00020      std
00021      Jacobian
00022      jac
00023      Tvector
00024      afun
00025      vx
00026      vy
00027      bool
00028      namespace
00029      CppAD
00030      const
00031      Taylor
00032      tx
00033      ty
00034      px
00035      py
00036 $$
00037 
00038 $section User Defined Atomic AD Functions$$
00039 $index CPPAD_USER_ATOMIC$$
00040 $index atomic, user function$$
00041 $index user, atomic function$$
00042 $index operation, user atomic$$
00043 $index function, user atomic$$
00044 
00045 $head Syntax$$
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 
00085 $head Partial Implementation$$
00086 The routines 
00087 $cref/forward/user_atomic/forward/$$,
00088 $cref/reverse/user_atomic/reverse/$$,
00089 $cref/for_jac_sparse/user_atomic/for_jac_sparse/$$,
00090 $cref/rev_jac_sparse/user_atomic/rev_jac_sparse/$$, and
00091 $cref/rev_hes_sparse/user_atomic/rev_hes_sparse/$$,
00092 must be defined by the user.
00093 The $icode forward$$ the routine, 
00094 for the case $icode%k% = 0%$$, must be implemented.
00095 Functions with the correct prototype,
00096 that just return $code false$$, 
00097 can be used for the other cases 
00098 (unless they are required by your calculations). 
00099 For example, you need not implement
00100 $icode forward$$ for the case $icode%k% == 2%$$ until you require
00101 forward mode calculation of second derivatives.
00102 
00103 $head CPPAD_USER_ATOMIC$$
00104 $index CPPAD_USER_ATOMIC$$
00105 The macro 
00106 $codei%
00107 CPPAD_USER_ATOMIC(%afun%, %Tvector%, %Base%, 
00108      %forward%, %reverse%, %for_jac_sparse%, %rev_jac_sparse%, %rev_hes_sparse%
00109 )
00110 %$$ 
00111 defines the $codei%AD<%Base%>%$$ routine $icode afun$$.
00112 This macro can be placed within a namespace 
00113 (not the $code CppAD$$ namespace) 
00114 but must be outside of any routine.
00115 
00116 $subhead Tvector$$
00117 The macro argument $icode Tvector$$ must be a
00118 $cref/simple vector template class/SimpleVector/$$.
00119 It determines the type of vectors used as arguments to the routine 
00120 $icode afun$$.
00121 
00122 $subhead Base$$
00123 The macro argument $icode Base$$ specifies the 
00124 $cref/base type/base_require/$$
00125 corresponding to $codei%AD<%Base>%$$ operation sequences.
00126 Calling the routine $icode afun$$ will add the operator defined
00127 by this macro to an $codei%AD<%Base>%$$ operation sequence.
00128 
00129 $head ok$$
00130 For all routines documented below,
00131 the return value $icode ok$$ has prototype
00132 $codei%
00133      bool %ok%
00134 %$$
00135 If it is $code true$$, the corresponding evaluation succeeded,
00136 otherwise it failed.
00137 
00138 $head id$$
00139 For all routines documented below,
00140 the argument $icode id$$ has prototype
00141 $codei%
00142      size_t %id%
00143 %$$
00144 Its value in all other calls is the same as in the corresponding
00145 call to $icode afun$$.
00146 It can be used to store and retrieve extra information about
00147 a specific call to $icode afun$$. 
00148 
00149 $head k$$
00150 For all routines documented below, the argument $icode k$$ has prototype
00151 $codei%
00152      size_t %k%
00153 %$$
00154 The value $icode%k%$$ is the order of the Taylor coefficient that
00155 we are evaluating ($cref/forward/user_atomic/forward/$$)
00156 or taking the derivative of ($cref/reverse/user_atomic/reverse/$$).
00157 
00158 $head n$$
00159 For all routines documented below, 
00160 the argument $icode n$$ has prototype
00161 $codei%
00162      size_t %n%
00163 %$$
00164 It is the size of the vector $icode ax$$ in the corresponding call to
00165 $icode%afun%(%id%, %ax%, %ay%)%$$; i.e.,
00166 the dimension of the domain space for $latex y = f(x)$$.
00167 
00168 $head m$$
00169 For all routines documented below, the argument $icode m$$ has prototype
00170 $codei%
00171      size_t %m%
00172 %$$
00173 It is the size of the vector $icode ay$$ in the corresponding call to
00174 $icode%afun%(%id%, %ax%, %ay%)%$$; i.e.,
00175 the dimension of the range space for $latex y = f(x)$$.
00176 
00177 $head tx$$
00178 For all routines documented below, 
00179 the argument $icode tx$$ has prototype
00180 $codei%
00181      const CppAD::vector<%Base%>& %tx%
00182 %$$
00183 and $icode%tx%.size() >= (%k% + 1) * %n%$$.
00184 For $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , k$$,
00185 we use the Taylor coefficient notation
00186 $latex \[
00187 \begin{array}{rcl}
00188      x_j^\ell & = & tx [ j * ( k + 1 ) + \ell ]
00189      \\
00190      X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^k t^k
00191 \end{array}
00192 \] $$
00193 If $icode%tx%.size() > (%k% + 1) * %n%$$,
00194 the other components of $icode tx$$ are not specified and should not be used.
00195 Note that superscripts represent an index for $latex x_j^\ell$$
00196 and an exponent for $latex t^\ell$$.
00197 Also note that the Taylor coefficients for $latex X(t)$$ correspond
00198 to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
00199 $latex \[
00200      x_j^\ell = \frac{1}{ \ell ! } X_j^{(\ell)} (0)
00201 \] $$
00202 
00203 $head ty$$
00204 In calls to $cref/forward/user_atomic/forward/$$, 
00205 the argument $icode ty$$ has prototype
00206 $codei%
00207      CppAD::vector<%Base%>& %ty%
00208 %$$
00209 while in calls to $cref/reverse/user_atomic/reverse/$$ it has prototype
00210 $codei%
00211      const CppAD::vector<%Base%>& %ty%
00212 %$$
00213 For all calls, $icode%tx%.size() >= (%k% + 1) * %m%$$.
00214 For $latex i = 0 , \ldots , m-1$$ and $latex \ell = 0 , \ldots , k$$,
00215 we use the Taylor coefficient notation
00216 $latex \[
00217 \begin{array}{rcl}
00218      y_i^\ell & = & ty [ i * ( k + 1 ) + \ell ]
00219      \\
00220      Y_i (t)  & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^k t^k + o ( t^k )
00221 \end{array}
00222 \] $$
00223 where $latex o( t^k ) / t^k \rightarrow 0$$ as $latex t \rightarrow 0$$.
00224 If $icode%ty%.size() > (%k% + 1) * %m%$$,
00225 the other components of $icode ty$$ are not specified and should not be used.
00226 Note that superscripts represent an index for $latex y_j^\ell$$
00227 and an exponent for $latex t^\ell$$.
00228 Also note that the Taylor coefficients for $latex Y(t)$$ correspond
00229 to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
00230 $latex \[
00231      y_j^\ell = \frac{1}{ \ell ! } Y_j^{(\ell)} (0)
00232 \] $$
00233 
00234 $head afun$$
00235 The macro argument $icode afun$$,
00236 is the name of the AD function corresponding to this atomic
00237 operation (as it is used in the source code).
00238 CppAD uses the other functions,
00239 where the arguments are vectors with elements of type $icode Base$$,
00240 to implement the function 
00241 $codei%
00242      %afun%(%id%, %ax%, %ay%)
00243 %$$
00244 where the argument are vectors with elements of type $codei%AD<%Base%>%$$.
00245 
00246 $subhead ax$$
00247 The $icode afun$$ argument $icode ax$$ has prototype
00248 $codei%
00249      const %Tvector%< AD<%Base%> >& %ax%
00250 %$$
00251 It is the argument vector $latex x \in B^n$$ 
00252 at which the $codei%AD<%Base%>%$$ version of 
00253 $latex y = f(x)$$ is to be evaluated.
00254 The dimension of the domain space for $latex y = f (x)$$
00255 is specified by $cref/n/user_atomic/n/$$ $codei%= %ax%.size()%$$,
00256 which must be greater than zero.
00257 
00258 $subhead ay$$
00259 The $icode afun$$ result $icode ay$$ has prototype
00260 $codei%
00261      %Tvector%< AD<%Base%> >& %ay%
00262 %$$
00263 The input values of its elements do not matter.
00264 Upon return, it is the $codei%AD<%Base%>%$$ version of the 
00265 result vector $latex y = f(x)$$.
00266 The dimension of the range space for $latex y = f (x)$$
00267 is specified by $cref/m/user_atomic/m/$$ $codei%= %ay%.size()%$$,
00268 which must be greater than zero.
00269 
00270 $subhead Parallel Mode$$
00271 $index parallel, user_atomic$$
00272 $index user_atomic, parallel$$
00273 The first call to 
00274 $codei%
00275      %afun%(%id%, %ax%, %ay%)
00276 %$$
00277 must not be in $cref/parallel/ta_in_parallel/$$ mode.
00278 In addition, the
00279 $cref/user_atomic clear/user_atomic/clear/$$
00280 routine cannot be called while in parallel mode.
00281 
00282 $head forward$$
00283 The macro argument $icode forward$$ is a
00284 user defined function
00285 $codei%
00286      %ok% = %forward%(%id%, %k%, %n%, %m%, %vx%, %vy%, %tx%, %ty%)
00287 %$$
00288 that computes results during a $cref/forward/Forward/$$ mode sweep.
00289 For this call, we are given the Taylor coefficients in $icode tx$$ 
00290 form order zero through $icode k$$,
00291 and the Taylor coefficients in  $icode ty$$ with order less than $icode k$$.
00292 The $icode forward$$ routine computes the 
00293 $icode k$$ order Taylor coefficients for $latex y$$ using the definition
00294 $latex Y(t) = f[ X(t) ]$$. 
00295 For example, for $latex i = 0 , \ldots , m-1$$,
00296 $latex \[
00297 \begin{array}{rcl}
00298 y_i^0 & = & Y(0) 
00299         = f_i ( x^0 )
00300 \\
00301 y_i^1 & = & Y^{(1)} ( 0 ) 
00302         = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 ) 
00303         = f_i^{(1)} ( x^0 ) x^1 
00304 \\
00305 y_i^2 
00306 & = & \frac{1}{2 !} Y^{(2)} (0)
00307 \\
00308 & = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 )
00309   +   \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 )
00310 \\
00311 & = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1 
00312   +    f_i^{(1)} ( x^0 ) x^2
00313 \end{array}
00314 \] $$
00315 Then, for $latex i = 0 , \ldots , m-1$$, it sets
00316 $latex \[
00317      ty [ i * (k + 1) + k ] = y_i^k
00318 \] $$
00319 The other components of $icode ty$$ must be left unchanged.
00320 
00321 $subhead Usage$$
00322 This routine is used,
00323 with $icode%vx%.size() > 0%$$ and $icode%k% == 0%$$,
00324 by calls to $icode afun$$.
00325 It is used,
00326 with $icode%vx%.size() = 0%$$ and
00327 $icode k$$ equal to the order of the derivative begin computed,
00328 by calls to $cref/forward/ForwardAny/$$.
00329 
00330 $subhead vx$$
00331 The $icode forward$$ argument $icode vx$$ has prototype
00332 $codei%
00333      const CppAD::vector<bool>& %vx%
00334 %$$
00335 The case $icode%vx%.size() > 0%$$ occurs 
00336 once for each call to $icode afun$$,
00337 during the call,
00338 and before any of the other callbacks corresponding to that call.
00339 Hence such a call can be used to cache information attached to 
00340 the corresponding $icode id$$
00341 (such as the elements of $icode vx$$).
00342 If $icode%vx%.size() > 0%$$ then
00343 $icode%k% == 0%$$, 
00344 $icode%vx%.size() >= %n%$$, and
00345 for $latex j = 0 , \ldots , n-1$$,
00346 $icode%vx%[%j%]%$$ is true if and only if
00347 $icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$.
00348 $pre
00349 
00350 $$
00351 If $icode%vx%.size() == 0%$$, 
00352 then $icode%vy%.size() == 0%$$ and neither of these vectors
00353 should be used.
00354 
00355 $subhead vy$$
00356 The $icode forward$$ argument $icode vy$$ has prototype
00357 $codei%
00358      CppAD::vector<bool>& %vy%
00359 %$$
00360 If $icode%vy%.size() == 0%$$, it should not be used.
00361 Otherwise, 
00362 $icode%k% == 0%$$ and $icode%vy%.size() >= %m%$$.
00363 The input values of the elements of $icode vy$$ do not matter.
00364 Upon return, for $latex j = 0 , \ldots , m-1$$,
00365 $icode%vy%[%i%]%$$ is true if and only if
00366 $icode%ay%[%j%]%$$ is a variable.
00367 (CppAD uses $icode vy$$ to reduce the necessary computations.)
00368 
00369 $head reverse$$
00370 The macro argument $icode reverse$$
00371 is a user defined function
00372 $codei%
00373      %ok% = %reverse%(%id%, %k%, %n%, %m%, %tx%, %ty%, %px%, %py%)
00374 %$$
00375 that computes results during a $cref/reverse/Reverse/$$ mode sweep. 
00376 The input value of the vectors $icode tx$$ and $icode ty$$
00377 contain Taylor coefficient, up to order $icode k$$,
00378 for $latex X(t)$$ and $latex Y(t)$$ respectively.
00379 We use the $latex \{ x_j^\ell \}$$ and $latex \{ y_i^\ell \}$$
00380 to denote these Taylor coefficients where the implicit range indices are
00381 $latex i = 0 , \ldots , m-1$$,
00382 $latex j = 0 , \ldots , n-1$$,
00383 $latex \ell = 0 , \ldots , k$$.
00384 Using the calculations done by $cref/forward/user_atomic/forward/$$,
00385 the Taylor coefficients $latex \{ y_i^\ell \}$$ are a function of the Taylor
00386 coefficients for $latex \{ x_j^\ell \}$$; i.e., given $latex y = f(x)$$
00387 we define the function
00388 $latex F : B^{n \times (k+1)} \rightarrow B^{m \times (k+1)}$$ by
00389 $latex \[
00390 y_i^\ell =  F_i^\ell ( \{ x_j^\ell \} )
00391 \] $$
00392 We use $latex G : B^{m \times (k+1)} \rightarrow B$$
00393 to denote an arbitrary scalar valued function of the Taylor coefficients for 
00394 $latex Y(t)$$ and write  $latex z = G( \{ y_i^\ell \} )$$.
00395 The $code reverse$$ routine
00396 is given the derivative of $latex z$$ with respect to
00397 $latex \{ y_i^\ell \}$$ and computes its derivative with respect
00398 to $latex \{ x_j^\ell \}$$.
00399 
00400 $subhead Usage$$
00401 This routine is used,
00402 with $icode%k% + 1%$$ equal to the order of the derivative being calculated,
00403 by calls to $cref/reverse/reverse_any/$$.
00404 
00405 $subhead py$$
00406 The $icode reverse$$ argument $icode py$$ has prototype
00407 $codei%
00408      const CppAD::vector<%Base%>& %py%
00409 %$$
00410 and $icode%py%.size() >= (%k% + 1) * %m%$$.
00411 For $latex i = 0 , \ldots , m-1$$ and $latex \ell = 0 , \ldots , k$$,
00412 $latex \[
00413      py[ i * (k + 1 ) + \ell ] = \partial G / \partial y_i^\ell
00414 \] $$
00415 If $icode%py%.size() > (%k% + 1) * %m%$$,
00416 the other components of $icode py$$ are not specified and should not be used.
00417 
00418 $subhead px$$
00419 We define the function 
00420 $latex \[
00421 H ( \{ x_j^\ell \} ) = G[ F( \{ x_j^\ell \} ) ]
00422 \] $$
00423 The $icode reverse$$ argument $icode px$$ has prototype
00424 $codei%
00425      CppAD::vector<%Base%>& %px%
00426 %$$
00427 and $icode%px%.size() >= (%k% + 1) * %n%$$.
00428 The input values of the elements of $icode px$$ do not matter.
00429 Upon return,
00430 for $latex j = 0 , \ldots , n-1$$ and $latex p = 0 , \ldots , k$$,
00431 $latex \[
00432 \begin{array}{rcl}
00433 px [ j * (k + 1) + p ] & = & \partial H / \partial x_j^p
00434 \\
00435 & = & 
00436 ( \partial G / \partial \{ y_i^\ell \} ) 
00437      ( \partial \{ y_i^\ell \} / \partial x_j^p )
00438 \\
00439 & = & 
00440 \sum_{i=0}^{m-1} \sum_{\ell=0}^k
00441 ( \partial G / \partial y_i^\ell ) ( \partial y_i^\ell / \partial x_j^p )
00442 \\
00443 & = & 
00444 \sum_{i=0}^{m-1} \sum_{\ell=p}^k
00445 py[ i * (k + 1 ) + \ell ] ( \partial F_i^\ell / \partial x_j^p )
00446 \end{array}
00447 \] $$
00448 Note that we have used the fact that for $latex \ell < p$$,
00449 $latex \partial F_i^\ell / \partial x_j^p = 0$$.
00450 If $icode%px%.size() > (%k% + 1) * %n%$$,
00451 the other components of $icode px$$ are not specified and should not be used.
00452 
00453 $head for_jac_sparse$$
00454 The macro argument $icode for_jac_sparse$$
00455 is a user defined function
00456 $codei%
00457      %ok% = %for_jac_sparse%(%id%, %n%, %m%, %q%, %r%, %s%)
00458 %$$
00459 that is used to compute results during a forward Jacobian sparsity sweep.
00460 For a fixed $latex n \times q$$ matrix $latex R$$,
00461 the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is
00462 $latex \[
00463      S(x) = f^{(1)} (x) * R
00464 \] $$
00465 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
00466 $icode for_jac_sparse$$ computes a sparsity pattern for $latex S(x)$$.
00467 
00468 $subhead Usage$$
00469 This routine is used by calls to $cref ForSparseJac$$.
00470 
00471 $subhead q$$
00472 The $icode for_jac_sparse$$ argument $icode q$$ has prototype
00473 $codei%
00474      size_t %q%
00475 %$$
00476 It specifies the number of columns in 
00477 $latex R \in B^{n \times q}$$ and the Jacobian 
00478 $latex S(x) \in B^{m \times q}$$. 
00479 
00480 $subhead r$$
00481 The $icode for_jac_sparse$$ argument $icode r$$ has prototype
00482 $codei%
00483      const CppAD::vector< std::set<size_t> >& %r%
00484 %$$
00485 and $icode%r%.size() >= %n%$$.
00486 For $latex j = 0 , \ldots , n-1$$,
00487 all the elements of $icode%r%[%j%]%$$ are between
00488 zero and $icode%q%-1%$$ inclusive.
00489 This specifies a sparsity pattern for the matrix $latex R$$.
00490 
00491 $subhead s$$
00492 The $icode for_jac_sparse$$ return value $icode s$$ has prototype
00493 $codei%
00494      CppAD::vector< std::set<size_t> >& %s%
00495 %$$
00496 and $icode%s%.size() >= %m%%$$.
00497 The input values of its sets do not matter. Upon return 
00498 for $latex i = 0 , \ldots , m-1$$,
00499 all the elements of $icode%s%[%i%]%$$ are between
00500 zero and $icode%q%-1%$$ inclusive.
00501 This represents a sparsity pattern for the matrix $latex S(x)$$.
00502 
00503 $head rev_jac_sparse$$
00504 The macro argument $icode rev_jac_sparse$$
00505 is a user defined function
00506 $codei%
00507      %ok% = %rev_jac_sparse%(%id%, %n%, %m%, %q%, %r%, %s%)
00508 %$$
00509 that is used to compute results during a reverse Jacobian sparsity sweep.
00510 For a fixed $latex q \times m$$ matrix $latex S$$,
00511 the Jacobian of $latex S * f( x )$$ with respect to $latex x \in B^n$$ is
00512 $latex \[
00513      R(x) = S * f^{(1)} (x)
00514 \] $$
00515 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex S$$,
00516 $icode rev_jac_sparse$$ computes a sparsity pattern for $latex R(x)$$.
00517 
00518 $subhead Usage$$
00519 This routine is used by calls to $cref RevSparseJac$$
00520 and to $cref optimize$$.
00521 
00522 
00523 $subhead q$$
00524 The $icode rev_jac_sparse$$ argument $icode q$$ has prototype
00525 $codei%
00526      size_t %q%
00527 %$$
00528 It specifies the number of rows in 
00529 $latex S \in B^{q \times m}$$ and the Jacobian 
00530 $latex R(x) \in B^{q \times n}$$. 
00531 
00532 $subhead s$$
00533 The $icode rev_jac_sparse$$ argument $icode s$$ has prototype
00534 $codei%
00535      const CppAD::vector< std::set<size_t> >& %s%
00536 %$$
00537 and $icode%s%.size() >= %m%$$.
00538 For $latex i = 0 , \ldots , m-1$$, 
00539 all the elements of $icode%s%[%i%]%$$
00540 are between zero and $icode%q%-1%$$ inclusive.
00541 This specifies a sparsity pattern for the matrix $latex S^\R{T}$$.
00542 
00543 $subhead r$$
00544 The $icode rev_jac_sparse$$ return value $icode r$$ has prototype
00545 $codei%
00546      CppAD::vector< std::set<size_t> >& %r%
00547 %$$
00548 and $icode%r%.size() >= %n%$$.
00549 The input values of its sets do not matter.
00550 Upon return for $latex j = 0 , \ldots , n-1$$,
00551 all the elements of $icode%r%[%j%]%$$
00552 are between zero and $icode%q%-1%$$ inclusive.
00553 This represents a sparsity pattern for the matrix $latex R(x)^\R{T}$$.
00554 
00555 $head rev_hes_sparse$$
00556 The macro argument $icode rev_hes_sparse$$
00557 is a user defined function
00558 $codei%
00559      %ok% = %rev_hes_sparse%(%id%, %n%, %m%, %q%, %r%, %s%, %t%, %u%, %v%)
00560 %$$
00561 There is an unspecified scalar valued function 
00562 $latex g : B^m \rightarrow B$$.
00563 Given a sparsity pattern for $latex R$$
00564 and information about the function $latex z = g(y)$$,
00565 this routine computes the sparsity pattern for
00566 $latex \[
00567      V(x) = (g \circ f)^{(2)}( x ) R
00568 \] $$
00569 
00570 $subhead Usage$$
00571 This routine is used by calls to $cref RevSparseHes$$.
00572 
00573 $subhead q$$
00574 The $icode rev_hes_sparse$$ argument $icode q$$ has prototype
00575 $codei%
00576      size_t %q%
00577 %$$
00578 It specifies the number of columns in the sparsity patterns.
00579 
00580 $subhead r$$
00581 The $icode rev_hes_sparse$$ argument $icode r$$ has prototype
00582 $codei%
00583      const CppAD::vector< std::set<size_t> >& %r%
00584 %$$
00585 and $icode%r%.size() >= %n%$$.
00586 For $latex j = 0 , \ldots , n-1$$,
00587 all the elements of $icode%r%[%j%]%$$ are between
00588 zero and $icode%q%-1%$$ inclusive.
00589 This specifies a sparsity pattern for the matrix $latex R \in B^{n \times q}$$.
00590 
00591 $subhead s$$
00592 The $icode rev_hes_sparse$$ argument $icode s$$ has prototype
00593 $codei%
00594      const CppAD::vector<bool>& %s%
00595 %$$
00596 and $icode%s%.size() >= %m%$$.
00597 This specifies a sparsity pattern for the matrix 
00598 $latex S(x) = g^{(1)} (y) \in B^{1 \times m}$$.
00599 
00600 $subhead t$$
00601 The $icode rev_hes_sparse$$ argument $icode t$$ has prototype
00602 $codei%
00603      CppAD::vector<bool>& %t%
00604 %$$
00605 and $icode%t%.size() >= %n%$$.
00606 The input values of its elements do not matter.
00607 Upon return it represents a sparsity pattern for the matrix 
00608 $latex T(x) \in B^{1 \times n}$$ defined by
00609 $latex \[
00610 T(x)  =  (g \circ f)^{(1)} (x) =  S(x) * f^{(1)} (x)
00611 \] $$
00612 
00613 $subhead u$$
00614 The $icode rev_hes_sparse$$ argument $icode u$$ has prototype
00615 $codei%
00616      const CppAD::vector< std::set<size_t> >& %u%
00617 %$$
00618 and $icode%u%.size() >= %m%$$.
00619 For $latex i = 0 , \ldots , m-1$$,
00620 all the elements of $icode%u%[%i%]%$$
00621 are between zero and $icode%q%-1%$$ inclusive.
00622 This specifies a sparsity pattern
00623 for the matrix $latex U(x) \in B^{m \times q}$$ defined by
00624 $latex \[
00625 \begin{array}{rcl}
00626 U(x) 
00627 & = & 
00628 \partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{u=0}
00629 \\
00630 & = & 
00631 \partial_u \{ g^{(1)} [ y + f^{(1)} (x) R u ] \}_{u=0}
00632 \\
00633 & = &
00634 g^{(2)} (y) f^{(1)} (x) R
00635 \end{array}
00636 \] $$
00637 
00638 $subhead v$$
00639 The $icode rev_hes_sparse$$ argument $icode v$$ has prototype
00640 $codei%
00641      CppAD::vector< std::set<size_t> >& %v%
00642 %$$
00643 and $icode%v%.size() >= %n%$$.
00644 The input values of its elements do not matter.
00645 Upon return, for $latex j = 0, \ldots , n-1$$,
00646 all the elements of $icode%v%[%j%]%$$
00647 are between zero and $icode%q%-1%$$ inclusive.
00648 This represents a sparsity pattern for the matrix 
00649 $latex V(x) \in B^{n \times q}$$ defined by
00650 $latex \[
00651 \begin{array}{rcl}
00652 V(x) 
00653 & = & 
00654 \partial_u [ \partial_x (g \circ f) ( x + R u )  ]_{u=0}
00655 \\
00656 & = &
00657 \partial_u [ (g \circ f)^{(1)}( x + R u )  ]_{u=0}
00658 \\
00659 & = &
00660 (g \circ f)^{(2)}( x ) R
00661 \\
00662 & = &
00663 f^{(1)} (x)^\R{T} g^{(2)} ( y ) f^{(1)} (x)  R
00664 +
00665 \sum_{i=1}^m [ g^{(1)} (y) ]_i \; f_i^{(2)} (x) R
00666 \\
00667 & = &
00668 f^{(1)} (x)^\R{T} U(x)
00669 +
00670 \sum_{i=1}^m S(x)_i \; f_i^{(2)} (x) R
00671 \end{array}
00672 \] $$
00673 
00674 $head clear$$
00675 User atomic functions hold onto static work space in order to
00676 increase speed by avoiding system memory allocation calls.
00677 The function call $codei%
00678      user_atomic<%Base%>::clear()
00679 %$$ 
00680 makes to work space $cref/available/ta_available/$$ to
00681 for other uses by the same thread.
00682 This should be called when you are done using the 
00683 user atomic functions for a specific value of $icode Base$$.
00684 
00685 $subhead Restriction$$
00686 The user atomic $code clear$$ routine cannot be called
00687 while in $cref/parallel/ta_in_parallel/$$ execution mode.
00688 
00689 $children%
00690      example/reciprocal.cpp%
00691      example/user_tan.cpp%
00692      example/mat_mul.cpp
00693 %$$
00694 $head Example$$
00695 
00696 $subhead Reciprocal$$
00697 The file $cref reciprocal.cpp$$ contains the simplest example and test
00698 of a user atomic operation.
00699 
00700 $subhead Tangent Function$$
00701 The file $cref user_tan.cpp$$ contains an example and test
00702 implementation of the tangent function as a user atomic operation.
00703 
00704 $subhead Matrix Multiplication$$
00705 The file  $cref mat_mul.cpp$$ contains an example and test
00706 implementation of matrix multiplication a a user atomic operation.
00707 
00708 $end
00709 ------------------------------------------------------------------------------
00710 */
00711 # include <set>
00712 # include <cppad/local/cppad_assert.hpp>
00713 
00714 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
00715 # include <cppad/thread_alloc.hpp>
00716 
00717 CPPAD_BEGIN_NAMESPACE
00718 /*!
00719 \defgroup user_atomic_hpp user_atomic.hpp
00720 \{
00721 \file user_atomic.hpp
00722 user defined atomic operations.
00723 */
00724 
00725 /*!
00726 \def CPPAD_USER_ATOMIC(afun, Tvector, 
00727      forward, reverse, for_jac_sparse, rev_jac_sparse, rev_hes_sparse 
00728 )
00729 Defines the function <tt>afun(id, ax, ay)</tt>  
00730 where \c id is \c ax and \c ay are vectors with <tt>AD<Base></tt> elements.
00731 
00732 \par Tvector
00733 the Simple Vector template class for this function.
00734 
00735 \par Base
00736 the base type for the atomic operation.
00737 
00738 \par afun 
00739 name of the CppAD defined function that corresponding to this operation.
00740 Note that \c afun, preceeded by a pound sign, 
00741 is a version of \c afun with quotes arround it.
00742 
00743 \par forward
00744 name of the user defined function that computes corresponding
00745 results during forward mode.
00746 
00747 \par reverse
00748 name of the user defined function that computes corresponding
00749 results during reverse mode.
00750 
00751 \par for_jac_sparse
00752 name of the user defined routine that computes corresponding
00753 results during forward mode jacobian sparsity sweeps.
00754 
00755 \par rev_jac_sparse
00756 name of the user defined routine that computes corresponding
00757 results during reverse mode jacobian sparsity sweeps.
00758 
00759 \par rev_hes_sparse
00760 name of the user defined routine that computes corresponding
00761 results during reverse mode Hessian sparsity sweeps.
00762 
00763 \par memory allocation
00764 Note that user_atomic is used as a static object, so its objects
00765 do note get deallocated until the program terminates. 
00766 */
00767 
00768 # define CPPAD_USER_ATOMIC(                                           \
00769      afun            ,                                                \
00770      Tvector         ,                                                \
00771      Base            ,                                                \
00772      forward         ,                                                \
00773      reverse         ,                                                \
00774      for_jac_sparse  ,                                                \
00775      rev_jac_sparse  ,                                                \
00776      rev_hes_sparse                                                   \
00777 )                                                                     \
00778 inline void afun (                                                    \
00779      size_t                               id ,                        \
00780      const Tvector< CppAD::AD<Base> >&    ax ,                        \
00781      Tvector< CppAD::AD<Base> >&          ay                          \
00782 )                                                                     \
00783 {    CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;                            \
00784      static CppAD::user_atomic<Base> fun(                             \
00785           #afun          ,                                            \
00786           forward        ,                                            \
00787           reverse        ,                                            \
00788           for_jac_sparse ,                                            \
00789           rev_jac_sparse ,                                            \
00790           rev_hes_sparse                                              \
00791      );                                                               \
00792      fun.ad(id, ax, ay);                                              \
00793 }
00794 
00795 /*!
00796 Class that actually implements the <tt>afun(id, ax, ay)</tt> calls.
00797 
00798 A new user_atomic object is generated each time the user invokes
00799 the CPPAD_USER_ATOMIC macro; see static object in that macro.
00800 */
00801 template <class Base>
00802 class user_atomic {
00803      /// type for user routine that computes forward mode results
00804      typedef bool (*F) (
00805           size_t                  id ,
00806           size_t                   k , 
00807           size_t                   n ,
00808           size_t                   m ,
00809           const vector<bool>&     vx ,
00810           vector<bool>&           vy ,
00811           const vector<Base>&     tx , 
00812           vector<Base>&           ty
00813      );
00814      /// type for user routine that computes reverse mode results
00815      typedef bool (*R) (
00816           size_t                  id ,
00817           size_t                   k , 
00818           size_t                   n , 
00819           size_t                   m , 
00820           const vector<Base>&     tx , 
00821           const vector<Base>&     ty ,
00822           vector<Base>&           px ,
00823           const vector<Base>&     py
00824      );
00825      /// type for user routine that computes forward mode Jacobian sparsity
00826      typedef bool (*FJS) (
00827           size_t                           id ,
00828           size_t                            n ,
00829           size_t                            m ,
00830           size_t                            q ,
00831           const vector< std::set<size_t> >& r ,
00832           vector< std::set<size_t>  >&      s
00833      );
00834      /// type for user routine that computes reverse mode Jacobian sparsity
00835      typedef bool (*RJS) (
00836           size_t                           id ,
00837           size_t                            n ,
00838           size_t                            m ,
00839           size_t                            q ,
00840           vector< std::set<size_t> >&       r ,
00841           const vector< std::set<size_t> >& s
00842      );
00843      /// type for user routine that computes reverse mode Hessian sparsity
00844      typedef bool (*RHS) (
00845           size_t                           id ,
00846           size_t                            n ,
00847           size_t                            m ,
00848           size_t                            q ,
00849           const vector< std::set<size_t> >& r ,
00850           const vector<bool>&               s ,
00851           vector<bool>&                     t ,
00852           const vector< std::set<size_t> >& u ,
00853           vector< std::set<size_t> >&       v 
00854      );
00855 private:
00856      /// users name for the AD version of this atomic operation
00857      const std::string     name_;
00858      /// user's implementation of forward mode
00859      const F                  f_;
00860      /// user's implementation of reverse mode
00861      const R                  r_;
00862      /// user's implementation of forward jacobian sparsity calculations
00863      const FJS              fjs_;
00864      /// user's implementation of reverse jacobian sparsity calculations
00865      const RJS              rjs_;
00866      /// user's implementation of reverse Hessian sparsity calculations
00867      const RHS              rhs_;
00868      /// index of this object in the vector of all objects in this class
00869      const size_t         index_;
00870 
00871      /// temporary work space used to avoid memory allocation/deallocation
00872      /// extra information to be passed to the functions
00873      vector<bool>            vx_[CPPAD_MAX_NUM_THREADS];
00874      vector<bool>            vy_[CPPAD_MAX_NUM_THREADS];
00875      vector<Base>             x_[CPPAD_MAX_NUM_THREADS];
00876      vector<Base>             y_[CPPAD_MAX_NUM_THREADS];
00877 
00878      /// List of all objects in this class.
00879      static std::vector<user_atomic *>& List(void)
00880      {    CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
00881           static std::vector<user_atomic *> list;
00882           return list;
00883      }
00884 public:
00885      /*!
00886      Constructor called for each invocation of CPPAD_USER_ATOMIC.
00887 
00888      Put this object in the list of all objects for this class and set
00889      the constant private data name_, f_, r_, and index_.
00890 
00891      \param afun
00892      is the user's name for the AD version of this atomic operation.
00893 
00894      \param f
00895      user routine that does forward mode calculations for this operation.
00896 
00897      \param r
00898      user routine that does reverse mode calculations for this operation.
00899 
00900      \param fjs
00901      user routine that does forward Jacobian sparsity calculations.
00902 
00903      \param rjs
00904      user routine that does reverse Jacobian sparsity calculations.
00905 
00906      \param rhs
00907      user routine that does reverse Hessian sparsity calculations.
00908 
00909      \par
00910      This constructor can not be used in parallel mode because
00911      the object it constructs is static. In addition, it changes the
00912      static object \c List.
00913      */
00914      user_atomic(const char* afun, F f, R r, FJS fjs, RJS rjs, RHS rhs) : 
00915      name_(afun)
00916      , f_(f)
00917      , r_(r)
00918      , fjs_(fjs)
00919      , rjs_(rjs)
00920      , rhs_(rhs)
00921      , index_( List().size() )
00922      {    CPPAD_ASSERT_KNOWN(
00923                ! thread_alloc::in_parallel() ,
00924                "First call to the function *afun is in parallel mode."
00925           );
00926           List().push_back(this);
00927      }
00928 
00929      /*!
00930      Implement the user call to <tt>afun(id, ax, ay)</tt>.
00931 
00932      \tparam ADVector
00933      A simple vector class with elements of type <code>AD<Base></code>.
00934 
00935      \param id
00936      extra information vector that is just passed through by CppAD,
00937      and possibly used by user's routines.
00938 
00939      \param ax
00940      is the argument vector for this call,
00941      <tt>ax.size()</tt> determines the number of arguments.
00942 
00943      \param ay
00944      is the result vector for this call,
00945      <tt>ay.size()</tt> determines the number of results.
00946 
00947      This routine is not \c const because it may modify the works
00948      space vectors \c x_ and \c y_.
00949      */
00950      template <class ADVector>
00951      void ad(size_t id, const ADVector& ax, ADVector& ay)
00952      {    size_t i, j, k;
00953           size_t n = ax.size();
00954           size_t m = ay.size();
00955           size_t thread = thread_alloc::thread_num();
00956 # ifndef NDEBUG
00957           bool ok;
00958           std::string msg = "user_atomoc: ";
00959 # endif
00960           vector <Base>& x  = x_[thread];
00961           vector <Base>& y  = y_[thread];
00962           vector <bool>& vx = vx_[thread];
00963           vector <bool>& vy = vy_[thread];
00964           //
00965           if( x.size() < n )
00966           {    x.resize(n);
00967                vx.resize(n);
00968           }
00969           if( y.size() < m )
00970           {    y.resize(m);
00971                vy.resize(m);
00972           }
00973           // 
00974           // Determine if we are going to have to tape this operation
00975           tape_id_t tape_id     = 0;
00976           ADTape<Base>* tape = CPPAD_NULL;
00977           for(j = 0; j < n; j++)
00978           {    x[j]   = ax[j].value_;
00979                vx[j]  = Variable( ax[j] );
00980                if ( (tape == CPPAD_NULL) & vx[j] )
00981                {    tape    = ax[j].tape_this();
00982                     tape_id = ax[j].tape_id_;
00983                }
00984 # ifndef NDEBUG
00985                ok  = (tape_id == 0) | Parameter(ax[j]);
00986                ok |= (tape_id == ax[j].tape_id_);
00987                if( ! ok )
00988                {    msg = msg + name_ + 
00989                     ": ax contains variables from different threads.";
00990                     CPPAD_ASSERT_KNOWN(false, msg.c_str());
00991                }
00992 # endif
00993           }
00994           // Use zero order forward mode to compute values
00995           k  = 0;
00996 # if NDEBUG
00997           f_(id, k, n, m, vx, vy, x, y);  
00998 # else
00999           ok = f_(id, k, n, m, vx, vy, x, y);  
01000           if( ! ok )
01001           {    msg = msg + name_ + ": ok returned false from "
01002                     "zero order forward mode calculation.";
01003                CPPAD_ASSERT_KNOWN(false, msg.c_str());
01004           }
01005 # endif
01006           // pass back values
01007           for(i = 0; i < m; i++)
01008           {    ay[i].value_ = y[i];
01009 
01010                // initialize entire vector as a constant (not on tape)
01011                ay[i].tape_id_ = 0;
01012                ay[i].taddr_   = 0;
01013           }
01014           // if tape is not null, ay is on the tape
01015           if( tape != CPPAD_NULL )
01016           {
01017                // Note the actual number of results is m
01018                CPPAD_ASSERT_UNKNOWN( NumRes(UserOp) == 0 );
01019                CPPAD_ASSERT_UNKNOWN( NumArg(UserOp) == 4 );
01020 
01021                // Begin operators corresponding to one user_atomic operation.
01022                // Put function index, domain size, range size, and id in tape
01023                tape->Rec_.PutArg(index_, id, n, m);
01024                tape->Rec_.PutOp(UserOp);
01025                // n + m operators follow for this one atomic operation
01026 
01027                // Now put the information for the argument vector in the tape
01028                CPPAD_ASSERT_UNKNOWN( NumRes(UsravOp) == 0 );
01029                CPPAD_ASSERT_UNKNOWN( NumRes(UsrapOp) == 0 );
01030                CPPAD_ASSERT_UNKNOWN( NumArg(UsravOp) == 1 );
01031                CPPAD_ASSERT_UNKNOWN( NumArg(UsrapOp) == 1 );
01032                for(j = 0; j < n; j++)
01033                {    if( vx[j] )
01034                     {    // information for an argument that is a variable
01035                          tape->Rec_.PutArg(ax[j].taddr_);
01036                          tape->Rec_.PutOp(UsravOp);
01037                     }
01038                     else
01039                     {    // information for an arugment that is parameter
01040                          addr_t p = tape->Rec_.PutPar(ax[j].value_);
01041                          tape->Rec_.PutArg(p);
01042                          tape->Rec_.PutOp(UsrapOp);
01043                     }
01044                }
01045 
01046                // Now put the information for the results in the tape
01047                CPPAD_ASSERT_UNKNOWN( NumArg(UsrrpOp) == 1 );
01048                CPPAD_ASSERT_UNKNOWN( NumRes(UsrrpOp) == 0 );
01049                CPPAD_ASSERT_UNKNOWN( NumArg(UsrrvOp) == 0 );
01050                CPPAD_ASSERT_UNKNOWN( NumRes(UsrrvOp) == 1 );
01051                for(i = 0; i < m; i++)
01052                {    if( vy[i] )
01053                     {    ay[i].taddr_    = tape->Rec_.PutOp(UsrrvOp);
01054                          ay[i].tape_id_  = tape_id;
01055                     }
01056                     else
01057                     {    addr_t p = tape->Rec_.PutPar(ay[i].value_);
01058                          tape->Rec_.PutArg(p);
01059                          tape->Rec_.PutOp(UsrrpOp);
01060                     }
01061                }
01062 
01063                // Put a duplicate UserOp at end of UserOp sequence
01064                tape->Rec_.PutArg(index_, id, n, m);
01065                tape->Rec_.PutOp(UserOp);
01066           } 
01067           return;
01068      }
01069 
01070      /// Name corresponding to a user_atomic object
01071      static const char* name(size_t index)
01072      {    return List()[index]->name_.c_str(); }
01073      /*!
01074      Link from forward mode sweep to users routine.
01075 
01076      \param index
01077      index for this function in the list of all user_atomic objects
01078 
01079      \param id
01080      extra information vector that is just passed through by CppAD,
01081      and possibly used by user's routines.
01082 
01083      \param k
01084      order for this forward mode calculation.
01085 
01086      \param n
01087      domain space size for this calcualtion.
01088 
01089      \param m
01090      range space size for this calculation.
01091 
01092      \param tx
01093      Taylor coefficients corresponding to \c x for this calculation.
01094 
01095      \param ty
01096      Taylor coefficient corresponding to \c y for this calculation
01097 
01098      See the forward mode in user's documentation for user_atomic 
01099      */
01100      static void forward(
01101           size_t                index , 
01102           size_t                   id ,
01103           size_t                    k ,
01104           size_t                    n , 
01105           size_t                    m , 
01106           const vector<Base>&      tx ,
01107           vector<Base>&            ty )
01108      {    
01109 # ifdef _OPENMP
01110           vector<bool> empty(0);
01111 # else
01112           static vector<bool> empty(0);
01113 # endif
01114           
01115           CPPAD_ASSERT_UNKNOWN( tx.size() >= n * k );
01116           CPPAD_ASSERT_UNKNOWN( ty.size() >= m * k );
01117 
01118           CPPAD_ASSERT_UNKNOWN(index < List().size() );
01119           user_atomic* op = List()[index];
01120 
01121           bool ok = op->f_(id, k, n, m, empty, empty, tx, ty);
01122           if( ! ok )
01123           {    std::stringstream ss;
01124                ss << k;
01125                std::string msg = "user_atomic: ";
01126                msg = msg + op->name_ + ": ok returned false from " + ss.str()
01127                    + " order forward mode calculation";
01128                CPPAD_ASSERT_KNOWN(false, msg.c_str());
01129           }
01130      }
01131 
01132 
01133      /*!
01134      Link from reverse mode sweep to users routine.
01135 
01136      \param index
01137      index in the list of all user_atomic objects
01138      corresponding to this function.
01139 
01140 
01141      \param id
01142      extra information vector that is just passed through by CppAD,
01143      and possibly used by user's routines.
01144 
01145      \param k
01146      order for this forward mode calculation.
01147 
01148      \param n
01149      domain space size for this calcualtion.
01150 
01151      \param m
01152      range space size for this calculation.
01153 
01154      \param tx
01155      Taylor coefficients corresponding to \c x for this calculation.
01156 
01157      \param ty
01158      Taylor coefficient corresponding to \c y for this calculation
01159 
01160      \param px
01161      Partials w.r.t. the \c x Taylor coefficients.
01162 
01163      \param py
01164      Partials w.r.t. the \c y Taylor coefficients.
01165 
01166      See reverse mode documentation for user_atomic 
01167      */
01168      static void reverse(
01169           size_t               index , 
01170           size_t                  id ,
01171           size_t                   k ,
01172           size_t                   n , 
01173           size_t                   m , 
01174           const vector<Base>&     tx ,
01175           const vector<Base>&     ty ,
01176           vector<Base>&           px ,
01177           const vector<Base>&     py )
01178      {
01179           CPPAD_ASSERT_UNKNOWN(index < List().size() );
01180           CPPAD_ASSERT_UNKNOWN( tx.size() >= n * k );
01181           CPPAD_ASSERT_UNKNOWN( px.size() >= n * k );
01182           CPPAD_ASSERT_UNKNOWN( ty.size() >= m * k );
01183           CPPAD_ASSERT_UNKNOWN( py.size() >= m * k );
01184           user_atomic* op = List()[index];
01185 
01186           bool ok = op->r_(id, k, n, m, tx, ty, px, py);
01187           if( ! ok )
01188           {    std::stringstream ss;
01189                ss << k;
01190                std::string msg = "user_atomic: ";
01191                msg = op->name_ + ": ok returned false from " + ss.str() 
01192                    + " order reverse mode calculation";
01193                CPPAD_ASSERT_KNOWN(false, msg.c_str());
01194           }
01195      }
01196 
01197      /*!
01198      Link from forward Jacobian sparsity sweep to users routine.
01199 
01200      \param index
01201      index in the list of all user_atomic objects
01202      corresponding to this function.
01203 
01204      \param id
01205      extra information vector that is just passed through by CppAD,
01206      and possibly used by user's routines.
01207 
01208      \param n
01209      domain space size for this calcualtion.
01210 
01211      \param m
01212      range space size for this calculation.
01213 
01214      \param q
01215      is the column dimension for the Jacobian sparsity partterns.
01216 
01217      \param r
01218      is the Jacobian sparsity pattern for the argument vector x
01219 
01220      \param s
01221      is the Jacobian sparsity pattern for the result vector y
01222      */
01223      static void for_jac_sparse(
01224           size_t                            index ,
01225           size_t                               id ,
01226           size_t                                n , 
01227           size_t                                m , 
01228           size_t                                q ,
01229           const vector< std::set<size_t> >&     r ,
01230           vector< std::set<size_t> >&           s )
01231      {
01232           CPPAD_ASSERT_UNKNOWN(index < List().size() );
01233           CPPAD_ASSERT_UNKNOWN( r.size() >= n );
01234           CPPAD_ASSERT_UNKNOWN( s.size() >= m );
01235           user_atomic* op = List()[index];
01236 
01237           bool ok = op->fjs_(id, n, m, q, r, s);
01238           if( ! ok )
01239           {    std::string msg = "user_atomic: ";
01240                msg = msg + op->name_ 
01241                    + ": ok returned false from for_jac_sparse calculation";
01242                CPPAD_ASSERT_KNOWN(false, msg.c_str());
01243           }
01244      }
01245 
01246      /*!
01247      Link from reverse Jacobian sparsity sweep to users routine.
01248 
01249      \param index
01250      index in the list of all user_atomic objects
01251      corresponding to this function.
01252 
01253      \param id
01254      extra information vector that is just passed through by CppAD,
01255      and possibly used by user's routines.
01256 
01257      \param n
01258      domain space size for this calcualtion.
01259 
01260      \param m
01261      range space size for this calculation.
01262 
01263      \param q
01264      is the row dimension for the Jacobian sparsity partterns.
01265 
01266      \param r
01267      is the Jacobian sparsity pattern for the argument vector x
01268 
01269      \param s
01270      is the Jacobian sparsity pattern for the result vector y
01271      */
01272      static void rev_jac_sparse(
01273           size_t                            index ,
01274           size_t                               id ,
01275           size_t                                n , 
01276           size_t                                m , 
01277           size_t                                q ,
01278           vector< std::set<size_t> >&           r ,
01279           const vector< std::set<size_t> >&     s )
01280      {
01281           CPPAD_ASSERT_UNKNOWN(index < List().size() );
01282           CPPAD_ASSERT_UNKNOWN( r.size() >= n );
01283           CPPAD_ASSERT_UNKNOWN( s.size() >= m );
01284           user_atomic* op = List()[index];
01285 
01286           bool ok = op->rjs_(id, n, m, q, r, s);
01287           if( ! ok )
01288           {    std::string msg = "user_atomic: ";
01289                msg = msg + op->name_ 
01290                    + ": ok returned false from rev_jac_sparse calculation";
01291                CPPAD_ASSERT_KNOWN(false, msg.c_str());
01292           }
01293      }
01294 
01295      /*!
01296      Link from reverse Hessian sparsity sweep to users routine.
01297 
01298      \param index
01299      index in the list of all user_atomic objects
01300      corresponding to this function.
01301 
01302      \param id
01303      extra information vector that is just passed through by CppAD,
01304      and possibly used by user's routines.
01305 
01306      \param n
01307      domain space size for this calcualtion.
01308 
01309      \param m
01310      range space size for this calculation.
01311 
01312      \param q
01313      is the column dimension for the sparsity partterns.
01314 
01315      \param r
01316      is the forward Jacobian sparsity pattern w.r.t the argument vector x
01317 
01318      \param s
01319      is the reverse Jacobian sparsity pattern w.r.t the result vector y.
01320 
01321      \param t
01322      is the reverse Jacobian sparsity pattern w.r.t the argument vector x.
01323 
01324      \param u
01325      is the Hessian sparsity pattern w.r.t the result vector y.
01326 
01327      \param v
01328      is the Hessian sparsity pattern w.r.t the argument vector x.
01329      */
01330      static void rev_hes_sparse(
01331           size_t                            index ,
01332           size_t                               id ,
01333           size_t                                n , 
01334           size_t                                m , 
01335           size_t                                q ,
01336           vector< std::set<size_t> >&           r ,
01337           const vector<bool>&                   s ,
01338           vector<bool>&                         t ,
01339           const vector< std::set<size_t> >&     u ,
01340           vector< std::set<size_t> >&           v )
01341      {
01342           CPPAD_ASSERT_UNKNOWN(index < List().size() );
01343           CPPAD_ASSERT_UNKNOWN( r.size() >= n );
01344           CPPAD_ASSERT_UNKNOWN( s.size() >= m );
01345           CPPAD_ASSERT_UNKNOWN( t.size() >= n );
01346           CPPAD_ASSERT_UNKNOWN( u.size() >= m );
01347           CPPAD_ASSERT_UNKNOWN( v.size() >= n );
01348           user_atomic* op = List()[index];
01349 
01350           bool ok = op->rhs_(id, n, m, q, r, s, t, u, v);
01351           if( ! ok )
01352           {    std::string msg = "user_atomic: ";
01353                msg = msg + op->name_ 
01354                    + ": ok returned false from rev_jac_sparse calculation";
01355                CPPAD_ASSERT_KNOWN(false, msg.c_str());
01356           }
01357      }
01358 
01359      /// Free static CppAD::vector memory used by this class (work space)
01360      static void clear(void)
01361      {    CPPAD_ASSERT_KNOWN(
01362                ! thread_alloc::in_parallel() ,
01363                "cannot use user_atomic clear during parallel execution"
01364           );
01365           size_t i = List().size();
01366           while(i--)
01367           {    size_t thread = CPPAD_MAX_NUM_THREADS;
01368                while(thread--)
01369                {
01370                     user_atomic* op = List()[i];
01371                     op->vx_[thread].clear();
01372                     op->vy_[thread].clear();
01373                     op->x_[thread].clear();
01374                     op->y_[thread].clear();
01375                }
01376           }
01377           return;
01378      }
01379 };
01380 
01381 /*! \} */
01382 CPPAD_END_NAMESPACE
01383 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines