CppAD: A C++ Algorithmic Differentiation Package
20130102
|
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