CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_ATOMIC_BASE_INCLUDED 00003 # define CPPAD_ATOMIC_BASE_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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 # include <set> 00017 # include <cppad/local/cppad_assert.hpp> 00018 // needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL 00019 # include <cppad/thread_alloc.hpp> 00020 00021 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00022 /*! 00023 \file atomic_base.hpp 00024 Base class for atomic user operations. 00025 */ 00026 00027 template <class Base> 00028 class atomic_base { 00029 // =================================================================== 00030 public: 00031 enum option_enum { bool_sparsity_enum, set_sparsity_enum}; 00032 private: 00033 // ------------------------------------------------------ 00034 // constants 00035 // 00036 /// index of this object in class_object 00037 const size_t index_; 00038 00039 // ----------------------------------------------------- 00040 // variables 00041 // 00042 /// sparsity pattern this object is currently using 00043 /// (set by constructor and option member functions) 00044 option_enum sparsity_; 00045 00046 /// temporary work space used afun, declared here to avoid memory 00047 /// allocation/deallocation for each call to afun 00048 vector<bool> afun_vx_[CPPAD_MAX_NUM_THREADS]; 00049 vector<bool> afun_vy_[CPPAD_MAX_NUM_THREADS]; 00050 vector<Base> afun_tx_[CPPAD_MAX_NUM_THREADS]; 00051 vector<Base> afun_ty_[CPPAD_MAX_NUM_THREADS]; 00052 00053 // ----------------------------------------------------- 00054 // static member functions 00055 // 00056 /// List of all the object in this class 00057 static std::vector<atomic_base *>& class_object(void) 00058 { CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL; 00059 static std::vector<atomic_base *> list_; 00060 return list_; 00061 } 00062 /// List of names for each object in this class 00063 static std::vector<std::string>& class_name(void) 00064 { CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL; 00065 static std::vector<std::string> list_; 00066 return list_; 00067 } 00068 // ===================================================================== 00069 public: 00070 // ----------------------------------------------------- 00071 // member functions not in user API 00072 // 00073 /// current sparsity setting 00074 option_enum sparsity(void) const 00075 { return sparsity_; } 00076 00077 /// Name corresponding to a base_atomic object 00078 const std::string& afun_name(void) const 00079 { return class_name()[index_]; } 00080 /* 00081 $begin atomic_ctor$$ 00082 $spell 00083 sq 00084 std 00085 afun 00086 arg 00087 CppAD 00088 bool 00089 ctor 00090 const 00091 matrix_mul.hpp 00092 $$ 00093 00094 $section Atomic Function Constructor$$ 00095 $index constructor, atomic function$$ 00096 $index atomic, function constructor$$ 00097 $index function, atomic constructor$$ 00098 00099 $head Syntax$$ 00100 $icode%atomic_user afun%(%ctor_arg_list%) 00101 %$$ 00102 $codei%atomic_base<%Base%>(%name%) 00103 %$$ 00104 00105 $head atomic_user$$ 00106 00107 $subhead ctor_arg_list$$ 00108 Is a list of arguments for the $icode atomic_user$$ constructor. 00109 00110 $subhead afun$$ 00111 The object $icode afun$$ must stay in scope for as long 00112 as the corresponding atomic function is used. 00113 This includes use by any $cref/ADFun<Base>/ADFun/$$ that 00114 has this $icode atomic_user$$ operation in its 00115 $cref/operation sequence/glossary/Operation/Sequence/$$. 00116 00117 $subhead Implementation$$ 00118 The user defined $icode atomic_user$$ class is a publicly derived class of 00119 $codei%atomic_base<%Base%>%$$. 00120 It should be declared as follows: 00121 $codei% 00122 class %atomic_user% : public CppAD::atomic_base<%Base%> { 00123 public: 00124 %atomic_user%(%ctor_arg_list%) : atomic_base<%Base%>(%name%) 00125 %...% 00126 }; 00127 %$$ 00128 where $icode ...$$ 00129 denotes the rest of the implementation of the derived class. 00130 This includes completing the constructor and 00131 all the virtual functions that have their 00132 $code atomic_base$$ implementations replaced by 00133 $icode atomic_user$$ implementations. 00134 00135 $head atomic_base$$ 00136 00137 $subhead Restrictions$$ 00138 The $code atomic_base$$ constructor cannot be called in 00139 $cref/parallel/ta_in_parallel/$$ mode. 00140 00141 $subhead Base$$ 00142 The template parameter determines the 00143 $icode Base$$ type for this $codei%AD<%Base%>%$$ atomic operation. 00144 00145 $subhead name$$ 00146 This $icode atomic_base$$ constructor argument has either of the 00147 following prototypes 00148 $codei% 00149 const char* %name% 00150 const std::string& %name% 00151 %$$ 00152 It is the name for this atomic function and is used for error reporting. 00153 The suggested value for $icode name$$ is $icode afun$$ or $icode atomic_user$$, 00154 i.e., the name of the corresponding atomic object or class. 00155 00156 $head Examples$$ 00157 00158 $subhead Define Constructor$$ 00159 The following are links to user atomic function constructor definitions: 00160 $cref%get_started.cpp%atomic_get_started.cpp%Constructor%$$, 00161 $cref%norm_sq.cpp%atomic_norm_sq.cpp%Constructor%$$, 00162 $cref%reciprocal.cpp%atomic_reciprocal.cpp%Constructor%$$, 00163 $cref%tangent.cpp%atomic_tangent.cpp%Constructor%$$, 00164 $cref%matrix_mul.hpp%atomic_matrix_mul.hpp%Constructor%$$. 00165 00166 $subhead Use Constructor$$ 00167 The following are links to user atomic function constructor uses: 00168 $cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Constructor%$$, 00169 $cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%Constructor%$$, 00170 $cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%Constructor%$$, 00171 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%Constructor%$$, 00172 $cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%Constructor%$$. 00173 00174 00175 $end 00176 */ 00177 /*! 00178 Base class for atomic_user functions. 00179 00180 \tparam Base 00181 This class is used for defining an AD<Base> atomic operation y = f(x). 00182 */ 00183 /// make sure user does not invoke the default constructor 00184 atomic_base(void) 00185 { CPPAD_ASSERT_KNOWN(false, 00186 "Attempt to use the atomic_base default constructor" 00187 ); 00188 } 00189 /*! 00190 Constructor 00191 00192 \param name 00193 name used for error reporting 00194 */ 00195 atomic_base( const std::string& name) : 00196 index_( class_object().size() ) , 00197 sparsity_( set_sparsity_enum ) 00198 { CPPAD_ASSERT_KNOWN( 00199 ! thread_alloc::in_parallel() , 00200 "atomic_base: constructor cannot be called in parallel mode." 00201 ); 00202 class_object().push_back(this); 00203 class_name().push_back(name); 00204 CPPAD_ASSERT_UNKNOWN( class_object().size() == class_name().size() ); 00205 } 00206 /// destructor informs CppAD that this atomic function with this index 00207 /// has dropped out of scope by setting its pointer to null 00208 virtual ~atomic_base(void) 00209 { CPPAD_ASSERT_UNKNOWN( class_object().size() > index_ ); 00210 // change object pointer to null, but leave name for error reporting 00211 class_object()[index_] = CPPAD_NULL; 00212 } 00213 /// atomic_base function object corresponding to a certain index 00214 static atomic_base* class_object(size_t index) 00215 { CPPAD_ASSERT_UNKNOWN( class_object().size() > index ); 00216 return class_object()[index]; 00217 } 00218 /// atomic_base function name corresponding to a certain index 00219 static const std::string& class_name(size_t index) 00220 { CPPAD_ASSERT_UNKNOWN( class_name().size() > index ); 00221 return class_name()[index]; 00222 } 00223 /* 00224 $begin atomic_option$$ 00225 $spell 00226 sq 00227 enum 00228 afun 00229 bool 00230 CppAD 00231 std 00232 typedef 00233 $$ 00234 00235 $section Set Atomic Function Options$$ 00236 $index atomic, options$$ 00237 $index options, atomic$$ 00238 00239 $head Syntax$$ 00240 $icode%afun%.option(%option_value%)%$$ 00241 00242 $head atomic_sparsity$$ 00243 $index atomic_sparsity$$ 00244 $index sparsity, atomic$$ 00245 You can used this option to set to type used for 00246 $icode afun$$ sparsity patterns. 00247 This does not apply individual calls to $icode afun$$, 00248 but rather all its uses between when the sparsity pattern is set and when 00249 it is changed. 00250 If neither the $code set_sparsity_enum$$ or 00251 $code bool_sparsity_enum$$ option is set, 00252 the type for $icode atomic_sparsity$$ is one of the two choices below 00253 (and otherwise unspecified). 00254 00255 $subhead bool_sparsity_enum$$ 00256 $index bool_sparsity_enum$$ 00257 If $icode option_value$$ is $codei%atomic_base<%Base%>::bool_sparsity_enum%$$, 00258 then the type used by $icode afun$$ for 00259 $cref/sparsity patterns/glossary/Sparsity Pattern/$$, 00260 (after the option is set) will be 00261 $codei% 00262 typedef CppAD::vector<bool> %atomic_sparsity% 00263 %$$ 00264 If $icode r$$ is a sparsity pattern 00265 for a matrix $latex R \in B^{p \times q}$$: 00266 $icode%r%.size() == %p% * %q%$$. 00267 00268 $subhead set_sparsity_enum$$ 00269 $index set_sparsity_enum$$ 00270 If $icode option_value$$ is $icode%atomic_base<%Base%>::set_sparsity_enum%$$, 00271 then the type used by $icode afun$$ for 00272 $cref/sparsity patterns/glossary/Sparsity Pattern/$$, 00273 (after the option is set) will be 00274 $codei% 00275 typedef CppAD::vector< std::set<size_t> > %atomic_sparsity% 00276 %$$ 00277 If $icode r$$ is a sparsity pattern 00278 for a matrix $latex R \in B^{p \times q}$$: 00279 $icode%r%.size() == %p%$$, and for $latex i = 0 , \ldots , p-1$$, 00280 the elements of $icode%r%[%i%]%$$ are between zero and $latex q-1$$ inclusive. 00281 00282 $end 00283 */ 00284 void option(enum option_enum option_value) 00285 { switch( option_value ) 00286 { case bool_sparsity_enum: 00287 case set_sparsity_enum: 00288 sparsity_ = option_value; 00289 break; 00290 00291 default: 00292 CPPAD_ASSERT_KNOWN( 00293 false, 00294 "atoic_base::option: option_value is not valid" 00295 ); 00296 } 00297 return; 00298 } 00299 /* 00300 ----------------------------------------------------------------------------- 00301 $begin atomic_afun$$ 00302 00303 $spell 00304 sq 00305 mul 00306 afun 00307 const 00308 CppAD 00309 $$ 00310 00311 $section Using AD Version of Atomic Function$$ 00312 $index atomic, use function$$ 00313 00314 $head Syntax$$ 00315 $icode%afun%(%ax%, %ay%)%$$ 00316 00317 $head Purpose$$ 00318 Given $icode ax$$, 00319 this call computes the corresponding value of $icode ay$$. 00320 If $codei%AD<%Base%>%$$ operations are being recorded, 00321 it enters the computation as an atomic operation in the recording; 00322 see $cref/start recording/Independent/Start Recording/$$. 00323 00324 $head ADVector$$ 00325 The type $icode ADVector$$ must be a 00326 $cref/simple vector class/SimpleVector/$$ with elements of type 00327 $codei%AD<%Base%>%$$; see $cref/Base/atomic_ctor/atomic_base/Base/$$. 00328 00329 $head afun$$ 00330 is a $cref/atomic_user/atomic_ctor/atomic_user/$$ object 00331 and this $icode afun$$ function call is implemented by the 00332 $cref/atomic_base/atomic_ctor/atomic_base/$$ class. 00333 00334 $head ax$$ 00335 This argument has prototype 00336 $codei% 00337 const %ADVector%& %ax% 00338 %$$ 00339 and size must be equal to $icode n$$. 00340 It specifies vector $latex x \in B^n$$ 00341 at which an $codei%AD<%Base%>%$$ version of 00342 $latex y = f(x)$$ is to be evaluated; see 00343 $cref/Base/atomic_ctor/atomic_base/Base/$$. 00344 00345 $head ay$$ 00346 This argument has prototype 00347 $codei% 00348 %ADVector%& %ay% 00349 %$$ 00350 and size must be equal to $icode m$$. 00351 The input values of its elements 00352 are not specified (must not matter). 00353 Upon return, it is an $codei%AD<%Base%>%$$ version of 00354 $latex y = f(x)$$. 00355 00356 $head Examples$$ 00357 The following files contain example uses of 00358 the AD version of atomic functions during recording: 00359 $cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Recording%$$, 00360 $cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%Recording%$$, 00361 $cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%Recording%$$, 00362 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%Recording%$$, 00363 $cref%matrix_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%Recording%$$. 00364 00365 $end 00366 ----------------------------------------------------------------------------- 00367 */ 00368 /*! 00369 Implement the user call to <tt>afun(ax, ay)</tt> and old_atomic call to 00370 <tt>afun(ax, ay, id)</tt>. 00371 00372 \tparam ADVector 00373 A simple vector class with elements of type <code>AD<Base></code>. 00374 00375 \param id 00376 optional extra information vector that is just passed through by CppAD, 00377 and used by old_atomic derived class (not other derived classes). 00378 This is an extra parameter to the virtual callbacks for old_atomic; 00379 see the set_id member function. 00380 00381 \param ax 00382 is the argument vector for this call, 00383 <tt>ax.size()</tt> determines the number of arguments. 00384 00385 \param ay 00386 is the result vector for this call, 00387 <tt>ay.size()</tt> determines the number of results. 00388 */ 00389 template <class ADVector> 00390 void operator()( 00391 const ADVector& ax , 00392 ADVector& ay , 00393 size_t id = 0 ) 00394 { size_t i, j; 00395 size_t n = ax.size(); 00396 size_t m = ay.size(); 00397 # ifndef NDEBUG 00398 bool ok; 00399 std::string msg = "atomic_base: " + afun_name() + ".eval: "; 00400 if( (n == 0) | (m == 0) ) 00401 { msg += "ax.size() or ay.size() is zero"; 00402 CPPAD_ASSERT_KNOWN(false, msg.c_str() ); 00403 } 00404 # endif 00405 size_t thread = thread_alloc::thread_num(); 00406 vector <Base>& tx = afun_tx_[thread]; 00407 vector <Base>& ty = afun_ty_[thread]; 00408 vector <bool>& vx = afun_vx_[thread]; 00409 vector <bool>& vy = afun_vy_[thread]; 00410 // 00411 if( vx.size() != n ) 00412 { vx.resize(n); 00413 tx.resize(n); 00414 } 00415 if( vy.size() != m ) 00416 { vy.resize(m); 00417 ty.resize(m); 00418 } 00419 // 00420 // Determine tape corresponding to variables in ax 00421 tape_id_t tape_id = 0; 00422 ADTape<Base>* tape = CPPAD_NULL; 00423 for(j = 0; j < n; j++) 00424 { tx[j] = ax[j].value_; 00425 vx[j] = Variable( ax[j] ); 00426 if( vx[j] ) 00427 { 00428 if( tape_id == 0 ) 00429 { tape = ax[j].tape_this(); 00430 tape_id = ax[j].tape_id_; 00431 CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL ); 00432 } 00433 # ifndef NDEBUG 00434 if( tape_id != ax[j].tape_id_ ) 00435 { msg += afun_name() + 00436 ": ax contains variables from different threads."; 00437 CPPAD_ASSERT_KNOWN(false, msg.c_str()); 00438 } 00439 # endif 00440 } 00441 } 00442 // Use zero order forward mode to compute values 00443 size_t p = 0, q = 0; 00444 set_id(id); 00445 # ifdef NDEBUG 00446 forward(p, q, vx, vy, tx, ty); 00447 # else 00448 ok = forward(p, q, vx, vy, tx, ty); 00449 if( ! ok ) 00450 { msg += afun_name() + ": ok is false for " 00451 "zero order forward mode calculation."; 00452 CPPAD_ASSERT_KNOWN(false, msg.c_str()); 00453 } 00454 # endif 00455 bool record_operation = false; 00456 for(i = 0; i < m; i++) 00457 { 00458 // pass back values 00459 ay[i].value_ = ty[i]; 00460 00461 // initialize entire vector parameters (not in tape) 00462 ay[i].tape_id_ = 0; 00463 ay[i].taddr_ = 0; 00464 00465 // we need to record this operation if 00466 // any of the eleemnts of ay are variables, 00467 record_operation |= vy[i]; 00468 } 00469 # ifndef NDEBUG 00470 if( record_operation & (tape == CPPAD_NULL) ) 00471 { msg += 00472 "all elements of vx are false but vy contains a true element"; 00473 CPPAD_ASSERT_KNOWN(false, msg.c_str() ); 00474 } 00475 # endif 00476 // if tape is not null, ay is on the tape 00477 if( record_operation ) 00478 { 00479 // Operator that marks beginning of this atomic operation 00480 CPPAD_ASSERT_UNKNOWN( NumRes(UserOp) == 0 ); 00481 CPPAD_ASSERT_UNKNOWN( NumArg(UserOp) == 4 ); 00482 tape->Rec_.PutArg(index_, id, n, m); 00483 tape->Rec_.PutOp(UserOp); 00484 00485 // Now put n operators, one for each element of arugment vector 00486 CPPAD_ASSERT_UNKNOWN( NumRes(UsravOp) == 0 ); 00487 CPPAD_ASSERT_UNKNOWN( NumRes(UsrapOp) == 0 ); 00488 CPPAD_ASSERT_UNKNOWN( NumArg(UsravOp) == 1 ); 00489 CPPAD_ASSERT_UNKNOWN( NumArg(UsrapOp) == 1 ); 00490 for(j = 0; j < n; j++) 00491 { if( vx[j] ) 00492 { // information for an argument that is a variable 00493 tape->Rec_.PutArg(ax[j].taddr_); 00494 tape->Rec_.PutOp(UsravOp); 00495 } 00496 else 00497 { // information for an arugment that is parameter 00498 addr_t par = tape->Rec_.PutPar(ax[j].value_); 00499 tape->Rec_.PutArg(par); 00500 tape->Rec_.PutOp(UsrapOp); 00501 } 00502 } 00503 00504 // Now put m operators, one for each element of result vector 00505 CPPAD_ASSERT_UNKNOWN( NumArg(UsrrpOp) == 1 ); 00506 CPPAD_ASSERT_UNKNOWN( NumRes(UsrrpOp) == 0 ); 00507 CPPAD_ASSERT_UNKNOWN( NumArg(UsrrvOp) == 0 ); 00508 CPPAD_ASSERT_UNKNOWN( NumRes(UsrrvOp) == 1 ); 00509 for(i = 0; i < m; i++) 00510 { if( vy[i] ) 00511 { ay[i].taddr_ = tape->Rec_.PutOp(UsrrvOp); 00512 ay[i].tape_id_ = tape_id; 00513 } 00514 else 00515 { addr_t par = tape->Rec_.PutPar(ay[i].value_); 00516 tape->Rec_.PutArg(par); 00517 tape->Rec_.PutOp(UsrrpOp); 00518 } 00519 } 00520 00521 // Put a duplicate UserOp at end of UserOp sequence 00522 tape->Rec_.PutArg(index_, id, n, m); 00523 tape->Rec_.PutOp(UserOp); 00524 } 00525 return; 00526 } 00527 /* 00528 ----------------------------------------------------------------------------- 00529 $begin atomic_forward$$ 00530 $spell 00531 sq 00532 mul.hpp 00533 hes 00534 afun 00535 vx 00536 vy 00537 ty 00538 Taylor 00539 const 00540 CppAD 00541 bool 00542 $$ 00543 00544 $section Atomic Forward Mode$$ 00545 $index atomic, forward callback$$ 00546 $index forward, atomic callback$$ 00547 $index forward, atomic virtual$$ 00548 00549 00550 $head Syntax$$ 00551 $icode%ok% = %afun%.forward(%p%, %q%, %vx%, %vy%, %tx%, %ty%)%$$ 00552 00553 $head Purpose$$ 00554 This virtual function is used by $cref atomic_afun$$ 00555 to evaluate function values. 00556 It is also used buy $cref/forward/Forward/$$ 00557 to compute function vales and derivatives. 00558 00559 $head Implementation$$ 00560 This virtual function must be defined by the 00561 $cref/atomic_user/atomic_ctor/atomic_user/$$ class. 00562 It can just return $icode%ok% == false%$$ 00563 (and not compute anything) for values 00564 of $icode%q% > 0%$$ that are greater than those used by your 00565 $cref/forward/Forward/$$ mode calculations. 00566 00567 $head p$$ 00568 The argument $icode p$$ has prototype 00569 $codei% 00570 size_t %p% 00571 %$$ 00572 It specifies the lowest order Taylor coefficient that we are evaluating. 00573 During calls to $cref atomic_afun$$, $icode%p% == 0%$$. 00574 00575 $head q$$ 00576 The argument $icode q$$ has prototype 00577 $codei% 00578 size_t %q% 00579 %$$ 00580 It specifies the highest order Taylor coefficient that we are evaluating. 00581 During calls to $cref atomic_afun$$, $icode%q% == 0%$$. 00582 00583 $head vx$$ 00584 The $code forward$$ argument $icode vx$$ has prototype 00585 $codei% 00586 const CppAD::vector<bool>& %vx% 00587 %$$ 00588 The case $icode%vx%.size() > 0%$$ only occurs while evaluating a call to 00589 $cref atomic_afun$$. 00590 In this case, 00591 $icode%p% == %q% == 0%$$, 00592 $icode%vx%.size() == %n%$$, and 00593 for $latex j = 0 , \ldots , n-1$$, 00594 $icode%vx%[%j%]%$$ is true if and only if 00595 $icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$ 00596 in the corresponding call to 00597 $codei% 00598 %afun%(%ax%, %ay%, %id%) 00599 %$$ 00600 If $icode%vx%.size() == 0%$$, 00601 then $icode%vy%.size() == 0%$$ and neither of these vectors 00602 should be used. 00603 00604 $head vy$$ 00605 The $code forward$$ argument $icode vy$$ has prototype 00606 $codei% 00607 CppAD::vector<bool>& %vy% 00608 %$$ 00609 If $icode%vy%.size() == 0%$$, it should not be used. 00610 Otherwise, 00611 $icode%q% == 0%$$ and $icode%vy%.size() == %m%$$. 00612 The input values of the elements of $icode vy$$ 00613 are not specified (must not matter). 00614 Upon return, for $latex j = 0 , \ldots , m-1$$, 00615 $icode%vy%[%i%]%$$ is true if and only if 00616 $icode%ay%[%i%]%$$ is a variable 00617 (CppAD uses $icode vy$$ to reduce the necessary computations). 00618 00619 $head tx$$ 00620 The argument $icode tx$$ has prototype 00621 $codei% 00622 const CppAD::vector<%Base%>& %tx% 00623 %$$ 00624 and $icode%tx%.size() == (%q%+1)*%n%$$. 00625 For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$, 00626 we use the Taylor coefficient notation 00627 $latex \[ 00628 \begin{array}{rcl} 00629 x_j^k & = & tx [ j * ( q + 1 ) + k ] 00630 \\ 00631 X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q 00632 \end{array} 00633 \] $$ 00634 Note that superscripts represent an index for $latex x_j^k$$ 00635 and an exponent for $latex t^k$$. 00636 Also note that the Taylor coefficients for $latex X(t)$$ correspond 00637 to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way: 00638 $latex \[ 00639 x_j^k = \frac{1}{ k ! } X_j^{(k)} (0) 00640 \] $$ 00641 00642 $head ty$$ 00643 The argument $icode ty$$ has prototype 00644 $codei% 00645 CppAD::vector<%Base%>& %ty% 00646 %$$ 00647 and $icode%tx%.size() == (%q%+1)*%m%$$. 00648 Upon return, 00649 For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$, 00650 $latex \[ 00651 \begin{array}{rcl} 00652 Y_i (t) & = & f_i [ X(t) ] 00653 \\ 00654 Y_i (t) & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q ) 00655 \\ 00656 ty [ i * ( q + 1 ) + k ] & = & y_i^k 00657 \end{array} 00658 \] $$ 00659 where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$. 00660 Note that superscripts represent an index for $latex y_j^k$$ 00661 and an exponent for $latex t^k$$. 00662 Also note that the Taylor coefficients for $latex Y(t)$$ correspond 00663 to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way: 00664 $latex \[ 00665 y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0) 00666 \] $$ 00667 If $latex p > 0$$, 00668 for $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p-1$$, 00669 the input of $icode ty$$ satisfies 00670 $latex \[ 00671 ty [ i * ( q + 1 ) + k ] = y_i^k 00672 \]$$ 00673 and hence the corresponding elements need not be recalculated. 00674 00675 $head ok$$ 00676 If the required results are calculated, $icode ok$$ should be true. 00677 Otherwise, it should be false. 00678 00679 $head Discussion$$ 00680 For example, suppose that $icode%q% == 2%$$, 00681 and you know how to compute the function $latex f(x)$$, 00682 its first derivative $latex f^{(1)} (x)$$, 00683 and it component wise Hessian $latex f_i^{(2)} (x)$$. 00684 Then you can compute $icode ty$$ using the following formulas: 00685 $latex \[ 00686 \begin{array}{rcl} 00687 y_i^0 & = & Y(0) 00688 = f_i ( x^0 ) 00689 \\ 00690 y_i^1 & = & Y^{(1)} ( 0 ) 00691 = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 ) 00692 = f_i^{(1)} ( x^0 ) x^1 00693 \\ 00694 y_i^2 00695 & = & \frac{1}{2 !} Y^{(2)} (0) 00696 \\ 00697 & = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 ) 00698 + \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 ) 00699 \\ 00700 & = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1 00701 + f_i^{(1)} ( x^0 ) x^2 00702 \end{array} 00703 \] $$ 00704 For $latex i = 0 , \ldots , m-1$$, and $latex k = 0 , 1 , 2$$, 00705 $latex \[ 00706 ty [ i * (q + 1) + k ] = y_i^k 00707 \] $$ 00708 00709 $head Examples$$ 00710 00711 $subhead Define forward$$ 00712 The following files contain example atomic $code forward$$ functions: 00713 $cref%get_started.cpp%atomic_get_started.cpp%forward%$$, 00714 $cref%norm_sq.cpp%atomic_norm_sq.cpp%forward%$$, 00715 $cref%reciprocal.cpp%atomic_reciprocal.cpp%forward%$$, 00716 $cref%tangent.cpp%atomic_tangent.cpp%forward%$$, 00717 $cref%matrix_mul.hpp%atomic_matrix_mul.hpp%forward%$$. 00718 00719 $subhead Use forward$$ 00720 The following are links to user atomic function forward uses: 00721 $cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%forward%$$, 00722 $cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%forward%$$, 00723 $cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%forward%$$, 00724 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%forward%$$, 00725 $cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%forward%$$. 00726 00727 00728 $end 00729 ----------------------------------------------------------------------------- 00730 */ 00731 /*! 00732 Link from atomic_base to forward mode 00733 00734 \param p [in] 00735 lowerest order for this forward mode calculation. 00736 00737 \param q [in] 00738 highest order for this forward mode calculation. 00739 00740 \param vx [in] 00741 if size not zero, which components of \c x are variables 00742 00743 \param vy [out] 00744 if size not zero, which components of \c y are variables 00745 00746 \param tx [in] 00747 Taylor coefficients corresponding to \c x for this calculation. 00748 00749 \param ty [out] 00750 Taylor coefficient corresponding to \c y for this calculation 00751 00752 See the forward mode in user's documentation for base_atomic 00753 */ 00754 virtual bool forward( 00755 size_t p , 00756 size_t q , 00757 const vector<bool>& vx , 00758 vector<bool>& vy , 00759 const vector<Base>& tx , 00760 vector<Base>& ty ) 00761 { return false; } 00762 /* 00763 ----------------------------------------------------------------------------- 00764 $begin atomic_reverse$$ 00765 $spell 00766 sq 00767 mul.hpp 00768 afun 00769 ty 00770 px 00771 py 00772 Taylor 00773 const 00774 CppAD 00775 $$ 00776 00777 $section Atomic Reverse Mode$$ 00778 $index atomic, reverse callback$$ 00779 $index reverse, atomic callback$$ 00780 $index reverse, atomic virtual$$ 00781 $spell 00782 bool 00783 $$ 00784 00785 $head Syntax$$ 00786 $icode%ok% = %afun%.reverse(%q%, %tx%, %ty%, %px%, %py%)%$$ 00787 00788 $head Purpose$$ 00789 This function is used by $cref/reverse/Reverse/$$ 00790 to compute derivatives. 00791 00792 $head Implementation$$ 00793 If you are using 00794 $cref/reverse/Reverse/$$ mode, 00795 this virtual function must be defined by the 00796 $cref/atomic_user/atomic_ctor/atomic_user/$$ class. 00797 It can just return $icode%ok% == false%$$ 00798 (and not compute anything) for values 00799 of $icode q$$ that are greater than those used by your 00800 $cref/reverse/Reverse/$$ mode calculations. 00801 00802 $head q$$ 00803 The argument $icode q$$ has prototype 00804 $codei% 00805 size_t %q% 00806 %$$ 00807 It specifies the highest order Taylor coefficient that 00808 computing the derivative of. 00809 00810 $head tx$$ 00811 The argument $icode tx$$ has prototype 00812 $codei% 00813 const CppAD::vector<%Base%>& %tx% 00814 %$$ 00815 and $icode%tx%.size() == (%q%+1)*%n%$$. 00816 For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$, 00817 we use the Taylor coefficient notation 00818 $latex \[ 00819 \begin{array}{rcl} 00820 x_j^k & = & tx [ j * ( q + 1 ) + k ] 00821 \\ 00822 X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q 00823 \end{array} 00824 \] $$ 00825 Note that superscripts represent an index for $latex x_j^k$$ 00826 and an exponent for $latex t^k$$. 00827 Also note that the Taylor coefficients for $latex X(t)$$ correspond 00828 to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way: 00829 $latex \[ 00830 x_j^k = \frac{1}{ k ! } X_j^{(k)} (0) 00831 \] $$ 00832 00833 $head ty$$ 00834 The argument $icode ty$$ has prototype 00835 $codei% 00836 const CppAD::vector<%Base%>& %ty% 00837 %$$ 00838 and $icode%tx%.size() == (%q%+1)*%m%$$. 00839 For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$, 00840 we use the Taylor coefficient notation 00841 $latex \[ 00842 \begin{array}{rcl} 00843 Y_i (t) & = & f_i [ X(t) ] 00844 \\ 00845 Y_i (t) & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q ) 00846 \\ 00847 y_i^k & = & ty [ i * ( q + 1 ) + k ] 00848 \end{array} 00849 \] $$ 00850 where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$. 00851 Note that superscripts represent an index for $latex y_j^k$$ 00852 and an exponent for $latex t^k$$. 00853 Also note that the Taylor coefficients for $latex Y(t)$$ correspond 00854 to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way: 00855 $latex \[ 00856 y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0) 00857 \] $$ 00858 00859 00860 $head F, G, H$$ 00861 We use the notation $latex \{ x_j^k \} \in B^{n \times (q+1)}$$ for 00862 $latex \[ 00863 \{ x_j^k \W{:} j = 0 , \ldots , n-1, k = 0 , \ldots , q \} 00864 \]$$ 00865 We use the notation $latex \{ y_i^k \} \in B^{m \times (q+1)}$$ for 00866 $latex \[ 00867 \{ y_i^k \W{:} i = 0 , \ldots , m-1, k = 0 , \ldots , q \} 00868 \]$$ 00869 We define the function 00870 $latex F : B^{n \times (q+1)} \rightarrow B^{m \times (q+1)}$$ by 00871 $latex \[ 00872 y_i^k = F_i^k [ \{ x_j^k \} ] 00873 \] $$ 00874 We use $latex G : B^{m \times (q+1)} \rightarrow B$$ 00875 to denote an arbitrary scalar valued function of $latex \{ y_i^k \}$$. 00876 We use $latex H : B^{n \times (q+1)} \rightarrow B$$ 00877 defined by 00878 $latex \[ 00879 H ( \{ x_j^k \} ) = G[ F( \{ x_j^k \} ) ] 00880 \] $$ 00881 00882 $head py$$ 00883 The argument $icode py$$ has prototype 00884 $codei% 00885 const CppAD::vector<%Base%>& %py% 00886 %$$ 00887 and $icode%py%.size() == m * (%q%+1)%$$. 00888 For $latex i = 0 , \ldots , m-1$$, $latex k = 0 , \ldots , q$$, 00889 $latex \[ 00890 py[ i * (q + 1 ) + k ] = \partial G / \partial y_i^k 00891 \] $$ 00892 00893 $subhead px$$ 00894 The $icode px$$ has prototype 00895 $codei% 00896 CppAD::vector<%Base%>& %px% 00897 %$$ 00898 and $icode%px%.size() == n * (%q%+1)%$$. 00899 The input values of the elements of $icode px$$ 00900 are not specified (must not matter). 00901 Upon return, 00902 for $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , q$$, 00903 $latex \[ 00904 \begin{array}{rcl} 00905 px [ j * (q + 1) + \ell ] & = & \partial H / \partial x_j^\ell 00906 \\ 00907 & = & 00908 ( \partial G / \partial \{ y_i^k \} ) 00909 ( \partial \{ y_i^k \} / \partial x_j^\ell ) 00910 \\ 00911 & = & 00912 \sum_{i=0}^{m-1} \sum_{k=0}^q 00913 ( \partial G / \partial y_i^k ) ( \partial y_i^k / \partial x_j^\ell ) 00914 \\ 00915 & = & 00916 \sum_{i=0}^{m-1} \sum_{k=\ell}^q 00917 py[ i * (q + 1 ) + k ] ( \partial F_i^k / \partial x_j^\ell ) 00918 \end{array} 00919 \] $$ 00920 Note that we have used the fact that for $latex k < \ell$$, 00921 $latex \partial F_i^k / \partial x_j^\ell = 0$$. 00922 00923 $head ok$$ 00924 The return value $icode ok$$ has prototype 00925 $codei% 00926 bool %ok% 00927 %$$ 00928 If it is $code true$$, the corresponding evaluation succeeded, 00929 otherwise it failed. 00930 00931 $head Examples$$ 00932 00933 $subhead Define reverse$$ 00934 The following files contain example atomic $code reverse$$ functions: 00935 $cref%norm_sq.cpp%atomic_norm_sq.cpp%reverse%$$, 00936 $cref%reciprocal.cpp%atomic_reciprocal.cpp%reverse%$$, 00937 $cref%tangent.cpp%atomic_tangent.cpp%reverse%$$, 00938 $cref%matrix_mul.hpp%atomic_matrix_mul.hpp%reverse%$$. 00939 00940 $subhead Use reverse$$ 00941 The following are links to user atomic function constructor uses: 00942 $cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%reverse%$$, 00943 $cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%reverse%$$, 00944 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%reverse%$$, 00945 $cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%reverse%$$. 00946 00947 $end 00948 ----------------------------------------------------------------------------- 00949 */ 00950 /*! 00951 Link from reverse mode sweep to users routine. 00952 00953 \param q [in] 00954 highest order for this reverse mode calculation. 00955 00956 \param tx [in] 00957 Taylor coefficients corresponding to \c x for this calculation. 00958 00959 \param ty [in] 00960 Taylor coefficient corresponding to \c y for this calculation 00961 00962 \param px [out] 00963 Partials w.r.t. the \c x Taylor coefficients. 00964 00965 \param py [in] 00966 Partials w.r.t. the \c y Taylor coefficients. 00967 00968 See atomic_reverse mode use documentation 00969 */ 00970 virtual bool reverse( 00971 size_t q , 00972 const vector<Base>& tx , 00973 const vector<Base>& ty , 00974 vector<Base>& px , 00975 const vector<Base>& py ) 00976 { return false; } 00977 /* 00978 -------------------------------------- --------------------------------------- 00979 $begin atomic_for_sparse_jac$$ 00980 $spell 00981 sq 00982 mul.hpp 00983 afun 00984 Jacobian 00985 jac 00986 const 00987 CppAD 00988 std 00989 bool 00990 std 00991 $$ 00992 00993 $section Atomic Forward Jacobian Sparsity Patterns$$ 00994 $index atomic, for_sparse_jac callback$$ 00995 $index for_sparse_jac, atomic callback$$ 00996 $index for_sparse_jac, atomic virtual$$ 00997 00998 $head Syntax$$ 00999 $icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%)%$$ 01000 01001 $head Purpose$$ 01002 This function is used by $cref ForSparseJac$$ to compute 01003 Jacobian sparsity patterns. 01004 For a fixed matrix $latex R \in B^{n \times q}$$, 01005 the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is 01006 $latex \[ 01007 S(x) = f^{(1)} (x) * R 01008 \] $$ 01009 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$, 01010 $code for_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$. 01011 01012 $head Implementation$$ 01013 If you are using $cref ForSparseJac$$, 01014 this virtual function must be defined by the 01015 $cref/atomic_user/atomic_ctor/atomic_user/$$ class. 01016 01017 $subhead q$$ 01018 The argument $icode q$$ has prototype 01019 $codei% 01020 size_t %q% 01021 %$$ 01022 It specifies the number of columns in 01023 $latex R \in B^{n \times q}$$ and the Jacobian 01024 $latex S(x) \in B^{m \times q}$$. 01025 01026 $subhead r$$ 01027 This argument has prototype 01028 $codei% 01029 const %atomic_sparsity%& %r% 01030 %$$ 01031 and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for 01032 $latex R \in B^{n \times q}$$. 01033 01034 $subhead s$$ 01035 This argument has prototype 01036 $codei% 01037 %atomic_sparsity%& %s% 01038 %$$ 01039 The input values of its elements 01040 are not specified (must not matter). 01041 Upon return, $icode s$$ is a 01042 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for 01043 $latex S(x) \in B^{m \times q}$$. 01044 01045 $head ok$$ 01046 The return value $icode ok$$ has prototype 01047 $codei% 01048 bool %ok% 01049 %$$ 01050 If it is $code true$$, the corresponding evaluation succeeded, 01051 otherwise it failed. 01052 01053 $head Examples$$ 01054 01055 $subhead Define for_sparse_jac$$ 01056 The following files contain example atomic $code for_sparse_jac$$ functions: 01057 $cref%norm_sq.cpp%atomic_norm_sq.cpp%for_sparse_jac%$$, 01058 $cref%reciprocal.cpp%atomic_reciprocal.cpp%for_sparse_jac%$$, 01059 $cref%tangent.cpp%atomic_tangent.cpp%for_sparse_jac%$$, 01060 $cref%matrix_mul.hpp%atomic_matrix_mul.hpp%for_sparse_jac%$$. 01061 01062 $subhead Use for_sparse_jac$$ 01063 The following are links to user atomic function constructor uses: 01064 $cref%norm_sq.cpp% 01065 atomic_norm_sq.cpp%Use Atomic Function%for_sparse_jac%$$, 01066 $cref%reciprocal.cpp% 01067 atomic_reciprocal.cpp%Use Atomic Function%for_sparse_jac%$$, 01068 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%for_sparse_jac%$$, 01069 $cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%for_sparse_jac%$$. 01070 01071 $end 01072 ----------------------------------------------------------------------------- 01073 */ 01074 /*! 01075 Link from forward Jacobian sparsity sweep to atomic_base 01076 01077 \param q 01078 is the column dimension for the Jacobian sparsity partterns. 01079 01080 \param r 01081 is the Jacobian sparsity pattern for the argument vector x 01082 01083 \param s 01084 is the Jacobian sparsity pattern for the result vector y 01085 */ 01086 virtual bool for_sparse_jac( 01087 size_t q , 01088 const vector< std::set<size_t> >& r , 01089 vector< std::set<size_t> >& s ) 01090 { return false; } 01091 virtual bool for_sparse_jac( 01092 size_t q , 01093 const vector<bool>& r , 01094 vector<bool>& s ) 01095 { return false; } 01096 /* 01097 -------------------------------------- --------------------------------------- 01098 $begin atomic_rev_sparse_jac$$ 01099 $spell 01100 sq 01101 mul.hpp 01102 rt 01103 afun 01104 Jacobian 01105 jac 01106 CppAD 01107 std 01108 bool 01109 const 01110 hes 01111 $$ 01112 01113 $section Atomic Reverse Jacobian Sparsity Patterns$$ 01114 $index atomic, rev_sparse_jac callback$$ 01115 $index rev_sparse_jac, atomic callback$$ 01116 $index rev_sparse_jac, atomic virtual$$ 01117 01118 $head Syntax$$ 01119 $icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%)%$$ 01120 01121 $head Purpose$$ 01122 This function is used by $cref RevSparseJac$$ to compute 01123 Jacobian sparsity patterns. 01124 For a fixed matrix $latex R \in B^{q \times m}$$, 01125 the Jacobian of $latex R * f( x )$$ with respect to $latex x \in B^q$$ is 01126 $latex \[ 01127 S(x) = R * f^{(1)} (x) 01128 \] $$ 01129 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$, 01130 $code rev_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$. 01131 01132 $head Implementation$$ 01133 If you are using $cref RevSparseJac$$, 01134 this virtual function must be defined by the 01135 $cref/atomic_user/atomic_ctor/atomic_user/$$ class. 01136 01137 $subhead q$$ 01138 The argument $icode q$$ has prototype 01139 $codei% 01140 size_t %q% 01141 %$$ 01142 It specifies the number of rows in 01143 $latex R \in B^{q \times m}$$ and the Jacobian 01144 $latex S(x) \in B^{q \times n}$$. 01145 01146 $subhead rt$$ 01147 This argument has prototype 01148 $codei% 01149 const %atomic_sparsity%& %rt% 01150 %$$ 01151 and is a 01152 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for 01153 $latex R^\R{T} \in B^{m \times q}$$. 01154 01155 $subhead st$$ 01156 This argument has prototype 01157 $codei% 01158 %atomic_sparsity%& %st% 01159 %$$ 01160 The input value of its elements 01161 are not specified (must not matter). 01162 Upon return, $icode s$$ is a 01163 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for 01164 $latex S(x)^\R{T} \in B^{n \times q}$$. 01165 01166 $head ok$$ 01167 The return value $icode ok$$ has prototype 01168 $codei% 01169 bool %ok% 01170 %$$ 01171 If it is $code true$$, the corresponding evaluation succeeded, 01172 otherwise it failed. 01173 01174 $head Examples$$ 01175 01176 $subhead Define rev_sparse_jac$$ 01177 The following files contain example atomic $code rev_sparse_jac$$ functions: 01178 $cref%norm_sq.cpp%atomic_norm_sq.cpp%rev_sparse_jac%$$, 01179 $cref%reciprocal.cpp%atomic_reciprocal.cpp%rev_sparse_jac%$$, 01180 $cref%tangent.cpp%atomic_tangent.cpp%rev_sparse_jac%$$, 01181 $cref%matrix_mul.hpp%atomic_matrix_mul.hpp%rev_sparse_jac%$$. 01182 01183 $subhead Use rev_sparse_jac$$ 01184 The following are links to user atomic function constructor uses: 01185 $cref%norm_sq.cpp% 01186 atomic_norm_sq.cpp%Use Atomic Function%rev_sparse_jac%$$, 01187 $cref%reciprocal.cpp% 01188 atomic_reciprocal.cpp%Use Atomic Function%rev_sparse_jac%$$, 01189 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%rev_sparse_jac%$$, 01190 $cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%rev_sparse_jac%$$. 01191 01192 $end 01193 ----------------------------------------------------------------------------- 01194 */ 01195 /*! 01196 Link from reverse Jacobian sparsity sweep to atomic_base 01197 01198 \param q [in] 01199 is the row dimension for the Jacobian sparsity partterns 01200 01201 \param rt [out] 01202 is the tansposed Jacobian sparsity pattern w.r.t to range variables y 01203 01204 \param st [in] 01205 is the tansposed Jacobian sparsity pattern for the argument variables x 01206 */ 01207 virtual bool rev_sparse_jac( 01208 size_t q , 01209 const vector< std::set<size_t> >& rt , 01210 vector< std::set<size_t> >& st ) 01211 { return false; } 01212 virtual bool rev_sparse_jac( 01213 size_t q , 01214 const vector<bool>& rt , 01215 vector<bool>& st ) 01216 { return false; } 01217 /* 01218 -------------------------------------- --------------------------------------- 01219 $begin atomic_rev_sparse_hes$$ 01220 $spell 01221 sq 01222 mul.hpp 01223 vx 01224 afun 01225 Jacobian 01226 jac 01227 CppAD 01228 std 01229 bool 01230 hes 01231 const 01232 $$ 01233 01234 $section Atomic Reverse Hessian Sparsity Patterns$$ 01235 $index atomic, rev_sparse_hes callback$$ 01236 $index rev_sparse_hes, atomic callback$$ 01237 $index rev_sparse_hes, atomic virtual$$ 01238 01239 $head Syntax$$ 01240 $icode%ok% = %afun%.rev_sparse_hes(%vx%, %s%, %t%, %q%, %r%, %u%, %v%)%$$ 01241 01242 $head Purpose$$ 01243 This function is used by $cref RevSparseHes$$ to compute 01244 Hessian sparsity patterns. 01245 There is an unspecified scalar valued function 01246 $latex g : B^m \rightarrow B$$. 01247 Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for 01248 $latex R \in B^{n \times q}$$, 01249 and information about the function $latex z = g(y)$$, 01250 this routine computes the sparsity pattern for 01251 $latex \[ 01252 V(x) = (g \circ f)^{(2)}( x ) R 01253 \] $$ 01254 01255 $head Implementation$$ 01256 If you are using and $cref RevSparseHes$$, 01257 this virtual function must be defined by the 01258 $cref/atomic_user/atomic_ctor/atomic_user/$$ class. 01259 01260 $subhead vx$$ 01261 The argument $icode vx$$ has prototype 01262 $codei% 01263 const CppAD:vector<bool>& %vx% 01264 %$$ 01265 $icode%vx%.size() == %n%$$, and 01266 for $latex j = 0 , \ldots , n-1$$, 01267 $icode%vx%[%j%]%$$ is true if and only if 01268 $icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$ 01269 in the corresponding call to 01270 $codei% 01271 %afun%(%ax%, %ay%, %id%) 01272 %$$ 01273 01274 $subhead s$$ 01275 The argument $icode s$$ has prototype 01276 $codei% 01277 const CppAD:vector<bool>& %s% 01278 %$$ 01279 and its size is $icode m$$. 01280 It is a sparsity pattern for 01281 $latex S(x) = g^{(1)} (y) \in B^{1 \times m}$$. 01282 01283 $subhead t$$ 01284 This argument has prototype 01285 $codei% 01286 CppAD:vector<bool>& %t% 01287 %$$ 01288 and its size is $icode m$$. 01289 The input values of its elements 01290 are not specified (must not matter). 01291 Upon return, $icode t$$ is a 01292 sparsity pattern for 01293 $latex T(x) \in B^{1 \times n}$$ where 01294 $latex \[ 01295 T(x) = (g \circ f)^{(1)} (x) = S(x) * f^{(1)} (x) 01296 \]$$ 01297 01298 $subhead q$$ 01299 The argument $icode q$$ has prototype 01300 $codei% 01301 size_t %q% 01302 %$$ 01303 It specifies the number of columns in 01304 $latex R \in B^{n \times q}$$, 01305 $latex U(x) \in B^{m \times q}$$, and 01306 $latex V(x) \in B^{n \times q}$$. 01307 01308 $subhead r$$ 01309 This argument has prototype 01310 $codei% 01311 const %atomic_sparsity%& %r% 01312 %$$ 01313 and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for 01314 $latex R \in B^{n \times q}$$. 01315 01316 $head u$$ 01317 This argument has prototype 01318 $codei% 01319 const %atomic_sparsity%& %u% 01320 %$$ 01321 and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for 01322 $latex U(x) \in B^{m \times q}$$ which is defined by 01323 $latex \[ 01324 \begin{array}{rcl} 01325 U(x) 01326 & = & 01327 \partial_u \{ \partial_y g[ y + f^{(1)} (x) R u ] \}_{u=0} 01328 \\ 01329 & = & 01330 \partial_u \{ g^{(1)} [ y + f^{(1)} (x) R u ] \}_{u=0} 01331 \\ 01332 & = & 01333 g^{(2)} (y) f^{(1)} (x) R 01334 \end{array} 01335 \] $$ 01336 01337 $subhead v$$ 01338 This argument has prototype 01339 $codei% 01340 %atomic_sparsity%& %v% 01341 %$$ 01342 The input value of its elements 01343 are not specified (must not matter). 01344 Upon return, $icode v$$ is a 01345 $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for 01346 $latex V(x) \in B^{n \times q}$$ which is defined by 01347 $latex \[ 01348 \begin{array}{rcl} 01349 V(x) 01350 & = & 01351 \partial_u [ \partial_x (g \circ f) ( x + R u ) ]_{u=0} 01352 \\ 01353 & = & 01354 \partial_u [ (g \circ f)^{(1)}( x + R u ) ]_{u=0} 01355 \\ 01356 & = & 01357 (g \circ f)^{(2)}( x ) R 01358 \\ 01359 & = & 01360 f^{(1)} (x)^\R{T} g^{(2)} ( y ) f^{(1)} (x) R 01361 + 01362 \sum_{i=1}^m g_i^{(1)} (y) \; f_i^{(2)} (x) R 01363 \\ 01364 & = & 01365 f^{(1)} (x)^\R{T} U(x) 01366 + 01367 \sum_{i=1}^m S_i (x) \; f_i^{(2)} (x) R 01368 \end{array} 01369 \] $$ 01370 01371 $head Examples$$ 01372 01373 $subhead Define rev_sparse_hes$$ 01374 The following files contain example atomic $code rev_sparse_hes$$ functions: 01375 $cref%norm_sq.cpp%atomic_norm_sq.cpp%rev_sparse_hes%$$, 01376 $cref%reciprocal.cpp%atomic_reciprocal.cpp%rev_sparse_hes%$$, 01377 $cref%tangent.cpp%atomic_tangent.cpp%rev_sparse_hes%$$, 01378 $cref%matrix_mul.hpp%atomic_matrix_mul.hpp%rev_sparse_hes%$$. 01379 01380 $subhead Use rev_sparse_hes$$ 01381 The following are links to user atomic function constructor uses: 01382 $cref%norm_sq.cpp% 01383 atomic_norm_sq.cpp%Use Atomic Function%rev_sparse_hes%$$, 01384 $cref%reciprocal.cpp% 01385 atomic_reciprocal.cpp%Use Atomic Function%rev_sparse_hes%$$, 01386 $cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%rev_sparse_hes%$$, 01387 $cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%rev_sparse_hes%$$. 01388 01389 $end 01390 ----------------------------------------------------------------------------- 01391 */ 01392 /*! 01393 Link from reverse Hessian sparsity sweep to base_atomic 01394 01395 \param vx [in] 01396 which componens of x are variables. 01397 01398 \param s [in] 01399 is the reverse Jacobian sparsity pattern w.r.t the result vector y. 01400 01401 \param t [out] 01402 is the reverse Jacobian sparsity pattern w.r.t the argument vector x. 01403 01404 \param q [in] 01405 is the column dimension for the sparsity partterns. 01406 01407 \param r [in] 01408 is the forward Jacobian sparsity pattern w.r.t the argument vector x 01409 01410 \param u [in] 01411 is the Hessian sparsity pattern w.r.t the result vector y. 01412 01413 \param v [out] 01414 is the Hessian sparsity pattern w.r.t the argument vector x. 01415 */ 01416 virtual bool rev_sparse_hes( 01417 const vector<bool>& vx , 01418 const vector<bool>& s , 01419 vector<bool>& t , 01420 size_t q , 01421 const vector< std::set<size_t> >& r , 01422 const vector< std::set<size_t> >& u , 01423 vector< std::set<size_t> >& v ) 01424 { return false; } 01425 virtual bool rev_sparse_hes( 01426 const vector<bool>& vx , 01427 const vector<bool>& s , 01428 vector<bool>& t , 01429 size_t q , 01430 const vector<bool>& r , 01431 const vector<bool>& u , 01432 vector<bool>& v ) 01433 { return false; } 01434 /* 01435 ------------------------------------------------------------------------------ 01436 $begin atomic_base_clear$$ 01437 $spell 01438 sq 01439 alloc 01440 $$ 01441 01442 $section Free Static Variables$$ 01443 $index free, atomic static$$ 01444 $index atomic, free static$$ 01445 $index static, free atomic$$ 01446 $index clear, atomic static$$ 01447 01448 $head Syntax$$ 01449 $codei%atomic_base<%Base%>::clear()%$$ 01450 01451 $head Purpose$$ 01452 Each $code atomic_base$$ objects holds onto work space in order to 01453 avoid repeated memory allocation calls and thereby increase speed 01454 (until it is deleted). 01455 If an the $code atomic_base$$ object is global or static because, 01456 the it does not get deleted. 01457 This is a problem when using 01458 $code thread_alloc$$ $cref/free_all/ta_free_all/$$ 01459 to check that all allocated memory has been freed. 01460 Calling this $code clear$$ function will free all the 01461 memory currently being held onto by the 01462 $codei%atomic_base<%Base%>%$$ class. 01463 01464 $head Future Use$$ 01465 If there is future use of an $code atomic_base$$ object, 01466 after a call to $code clear$$, 01467 the work space will be reallocated and held onto. 01468 01469 $head Restriction$$ 01470 This routine cannot be called 01471 while in $cref/parallel/ta_in_parallel/$$ execution mode. 01472 01473 $end 01474 ------------------------------------------------------------------------------ 01475 */ 01476 /*! 01477 Free all thread_alloc static memory held by atomic_base (avoids reallocations). 01478 (This does not include class_object() which is an std::vector.) 01479 */ 01480 /// Free vector memory used by this class (work space) 01481 static void clear(void) 01482 { CPPAD_ASSERT_KNOWN( 01483 ! thread_alloc::in_parallel() , 01484 "cannot use atomic_base clear during parallel execution" 01485 ); 01486 size_t i = class_object().size(); 01487 while(i--) 01488 { size_t thread = CPPAD_MAX_NUM_THREADS; 01489 while(thread--) 01490 { 01491 atomic_base* op = class_object()[i]; 01492 if( op != CPPAD_NULL ) 01493 { op->afun_vx_[thread].clear(); 01494 op->afun_vy_[thread].clear(); 01495 op->afun_tx_[thread].clear(); 01496 op->afun_ty_[thread].clear(); 01497 } 01498 } 01499 } 01500 return; 01501 } 01502 // ------------------------------------------------------------------------- 01503 /*! 01504 Set value of id (used by deprecated old_atomic class) 01505 01506 This function is called just before calling any of the virtual funcitons 01507 and has the corresponding id of the corresponding virtual call. 01508 */ 01509 virtual void set_id(size_t id) 01510 { } 01511 // --------------------------------------------------------------------------- 01512 }; 01513 } // END_CPPAD_NAMESPACE 01514 # endif