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