CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_STORE_OP_INCLUDED 00003 # define CPPAD_STORE_OP_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 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00017 /*! 00018 \file store_op.hpp 00019 Changing the current value of a VecAD element. 00020 */ 00021 /*! 00022 Shared documentation for zero order forward implementation of 00023 op = StppOp, StpvOp, StvpOp, or StvvOp (not called). 00024 00025 The C++ source code corresponding to this operation is 00026 \verbatim 00027 v[x] = y 00028 \endverbatim 00029 where v is a VecAD<Base> vector, x is an AD<Base> object, 00030 and y is AD<Base> or Base objects. 00031 We define the index corresponding to v[x] by 00032 \verbatim 00033 i_v_x = index_by_ind[ arg[0] + i_vec ] 00034 \endverbatim 00035 where i_vec is defined under the heading arg[1] below: 00036 00037 \tparam Base 00038 base type for the operator; i.e., this operation was recorded 00039 using AD<Base> and computations by this routine are done using type Base. 00040 00041 \param i_z 00042 is the index corresponding to the previous variable on the tape 00043 (only used for error checking). 00044 00045 \param arg 00046 \n 00047 arg[0] 00048 \n 00049 is the offset of this VecAD vector relative to the beginning 00050 of the isvar_by_ind and index_by_ind arrays. 00051 \n 00052 \n 00053 arg[1] 00054 \n 00055 If this is a StppOp or StpvOp operation (if x is a parameter), 00056 i_vec is defined by 00057 \verbatim 00058 i_vec = arg[1] 00059 \endverbatim 00060 If this is a StvpOp or StvvOp operation (if x is a variable), 00061 i_vec is defined by 00062 \verbatim 00063 i_vec = floor( taylor[ arg[1] * nc_taylor + 0 ] ) 00064 \endverbatim 00065 where floor(c) is the greatest integer less that or equal c. 00066 \n 00067 \n 00068 arg[2] 00069 \n 00070 index corresponding to the third operand for this operator; 00071 i.e. the index corresponding to y. 00072 00073 \param num_par 00074 is the total number of parameters on the tape 00075 (only used for error checking). 00076 00077 \param nc_taylor 00078 number of columns in the matrix containing the Taylor coefficients. 00079 00080 \param taylor 00081 In StvpOp and StvvOp cases, <code><taylor[ arg[1] * nc_taylor + 0 ]</code> 00082 is used to compute the index in the definition of i_vec above. 00083 00084 \param isvar_by_ind 00085 If y is a varable (StpvOp and StvvOp cases), 00086 <code>isvar_by_ind[ arg[0] + i_vec ] </code> is set to true. 00087 Otherwise y is a paraemter (StppOp and StvpOp cases) and 00088 <code>isvar_by_ind[ arg[0] + i_vec ] </code> is set to false. 00089 00090 \param index_by_ind 00091 <code>index_by_ind[ arg[0] - 1 ]</code> 00092 is the number of elements in the user vector containing this element. 00093 The value <code>index_by_ind[ arg[0] + i_vec]</code> 00094 is set equal to arg[2]. 00095 00096 \par Check User Errors 00097 \li Check that the index is with in range; i.e. 00098 <code>i_vec < index_by_ind[ arg[0] - 1 ]</code> 00099 Note that, if x is a parameter, 00100 the corresponding vector index and it does not change. 00101 In this case, the error above should be detected during tape recording. 00102 00103 \par Checked Assertions 00104 \li NumArg(op) == 3 00105 \li NumRes(op) == 0 00106 \li 0 < arg[0] 00107 \li if y is a parameter, arg[2] < num_par 00108 \li if x is a variable, arg[1] <= i_z 00109 \li if y is a variable, arg[2] <= i_z 00110 */ 00111 template <class Base> 00112 inline void forward_store_op_0( 00113 size_t i_z , 00114 const addr_t* arg , 00115 size_t num_par , 00116 size_t nc_taylor , 00117 Base* taylor , 00118 bool* isvar_by_ind , 00119 size_t* index_by_ind ) 00120 { 00121 // This routine is only for documentaiton, it should not be used 00122 CPPAD_ASSERT_UNKNOWN( false ); 00123 } 00124 /*! 00125 Shared documnetation for sparsity operations corresponding to 00126 op = StpvOp or StvvOp (not called). 00127 00128 <!-- define sparse_store_op --> 00129 The C++ source code corresponding to this operation is 00130 \verbatim 00131 v[x] = y 00132 \endverbatim 00133 where v is a VecAD<Base> vector, x is an AD<Base> object, 00134 and y is AD<Base> or Base objects. 00135 We define the index corresponding to v[x] by 00136 \verbatim 00137 i_v_x = combined[ arg[0] + i_vec ] 00138 \endverbatim 00139 where i_vec is defined under the heading \a arg[1] below: 00140 00141 \tparam Vector_set 00142 is the type used for vectors of sets. It can be either 00143 \c sparse_pack, \c sparse_set, or \c sparse_list. 00144 00145 \param op 00146 is the code corresponding to this operator; i.e., StpvOp or StvvOp 00147 (only used for error checking). 00148 00149 \param arg 00150 \n 00151 \a arg[0] 00152 is the offset corresponding to this VecAD vector in the combined array. 00153 \n 00154 \n 00155 \a arg[2] 00156 \n 00157 The set with index \a arg[2] in \a var_sparsity 00158 is the sparsity pattern corresponding to y. 00159 (Note that \a arg[2] > 0 because y is a variable.) 00160 00161 \param num_combined 00162 is the total number of elements in the VecAD address array. 00163 00164 \param combined 00165 \a combined [ arg[0] - 1 ] 00166 is the index of the set in \a vecad_sparsity corresponding 00167 to the sparsity pattern for the vector v. 00168 We use the notation i_v below which is defined by 00169 \verbatim 00170 i_v = combined[ \a arg[0] - 1 ] 00171 \endverbatim 00172 00173 \param var_sparsity 00174 The set with index \a arg[2] in \a var_sparsity 00175 is the sparsity pattern for y. 00176 This is an input for forward mode operations. 00177 For reverse mode operations: 00178 The sparsity pattern for v is added to the spartisy pattern for y. 00179 00180 \param vecad_sparsity 00181 The set with index \a i_v in \a vecad_sparsity 00182 is the sparsity pattern for v. 00183 This is an input for reverse mode operations. 00184 For forward mode operations, the sparsity pattern for y is added 00185 to the sparsity pattern for the vector v. 00186 00187 \par Checked Assertions 00188 \li NumArg(op) == 3 00189 \li NumRes(op) == 0 00190 \li 0 < \a arg[0] 00191 \li \a arg[0] < \a num_combined 00192 \li \a arg[2] < \a var_sparsity.n_set() 00193 \li i_v < \a vecad_sparsity.n_set() 00194 <!-- end sparse_store_op --> 00195 */ 00196 template <class Vector_set> 00197 inline void sparse_store_op( 00198 OpCode op , 00199 const addr_t* arg , 00200 size_t num_combined , 00201 const size_t* combined , 00202 Vector_set& var_sparsity , 00203 Vector_set& vecad_sparsity ) 00204 { 00205 // This routine is only for documentaiton, it should not be used 00206 CPPAD_ASSERT_UNKNOWN( false ); 00207 } 00208 00209 00210 /*! 00211 Zero order forward mode implementation of op = StppOp. 00212 00213 \copydetails forward_store_op_0 00214 */ 00215 template <class Base> 00216 inline void forward_store_pp_op_0( 00217 size_t i_z , 00218 const addr_t* arg , 00219 size_t num_par , 00220 size_t nc_taylor , 00221 Base* taylor , 00222 bool* isvar_by_ind , 00223 size_t* index_by_ind ) 00224 { size_t i_vec = arg[1]; 00225 00226 // Because the index is a parameter, this indexing error should be 00227 // caught and reported to the user when the tape is recording. 00228 CPPAD_ASSERT_UNKNOWN( i_vec < index_by_ind[ arg[0] - 1 ] ); 00229 00230 CPPAD_ASSERT_UNKNOWN( NumArg(StppOp) == 3 ); 00231 CPPAD_ASSERT_UNKNOWN( NumRes(StppOp) == 0 ); 00232 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00233 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par ); 00234 00235 isvar_by_ind[ arg[0] + i_vec ] = false; 00236 index_by_ind[ arg[0] + i_vec ] = arg[2]; 00237 } 00238 00239 /*! 00240 Zero order forward mode implementation of op = StpvOp. 00241 00242 \copydetails forward_store_op_0 00243 */ 00244 template <class Base> 00245 inline void forward_store_pv_op_0( 00246 size_t i_z , 00247 const addr_t* arg , 00248 size_t num_par , 00249 size_t nc_taylor , 00250 Base* taylor , 00251 bool* isvar_by_ind , 00252 size_t* index_by_ind ) 00253 { size_t i_vec = arg[1]; 00254 00255 // Because the index is a parameter, this indexing error should be 00256 // caught and reported to the user when the tape is recording. 00257 CPPAD_ASSERT_UNKNOWN( i_vec < index_by_ind[ arg[0] - 1 ] ); 00258 00259 CPPAD_ASSERT_UNKNOWN( NumArg(StpvOp) == 3 ); 00260 CPPAD_ASSERT_UNKNOWN( NumRes(StpvOp) == 0 ); 00261 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00262 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) <= i_z ); 00263 00264 isvar_by_ind[ arg[0] + i_vec ] = true; 00265 index_by_ind[ arg[0] + i_vec ] = arg[2]; 00266 } 00267 00268 /*! 00269 Zero order forward mode implementation of op = StvpOp. 00270 00271 \copydetails forward_store_op_0 00272 */ 00273 template <class Base> 00274 inline void forward_store_vp_op_0( 00275 size_t i_z , 00276 const addr_t* arg , 00277 size_t num_par , 00278 size_t nc_taylor , 00279 Base* taylor , 00280 bool* isvar_by_ind , 00281 size_t* index_by_ind ) 00282 { 00283 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) <= i_z ); 00284 size_t i_vec = Integer( taylor[ arg[1] * nc_taylor + 0 ] ); 00285 CPPAD_ASSERT_KNOWN( 00286 i_vec < index_by_ind[ arg[0] - 1 ] , 00287 "VecAD: index during zero order forward sweep is out of range" 00288 ); 00289 00290 CPPAD_ASSERT_UNKNOWN( NumArg(StvpOp) == 3 ); 00291 CPPAD_ASSERT_UNKNOWN( NumRes(StvpOp) == 0 ); 00292 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00293 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par ); 00294 00295 isvar_by_ind[ arg[0] + i_vec ] = false; 00296 index_by_ind[ arg[0] + i_vec ] = arg[2]; 00297 } 00298 00299 /*! 00300 Zero order forward mode implementation of op = StvvOp. 00301 00302 \copydetails forward_store_op_0 00303 */ 00304 template <class Base> 00305 inline void forward_store_vv_op_0( 00306 size_t i_z , 00307 const addr_t* arg , 00308 size_t num_par , 00309 size_t nc_taylor , 00310 Base* taylor , 00311 bool* isvar_by_ind , 00312 size_t* index_by_ind ) 00313 { 00314 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) <= i_z ); 00315 size_t i_vec = Integer( taylor[ arg[1] * nc_taylor + 0 ] ); 00316 CPPAD_ASSERT_KNOWN( 00317 i_vec < index_by_ind[ arg[0] - 1 ] , 00318 "VecAD: index during zero order forward sweep is out of range" 00319 ); 00320 00321 CPPAD_ASSERT_UNKNOWN( NumArg(StvpOp) == 3 ); 00322 CPPAD_ASSERT_UNKNOWN( NumRes(StvpOp) == 0 ); 00323 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00324 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) <= i_z ); 00325 00326 isvar_by_ind[ arg[0] + i_vec ] = true; 00327 index_by_ind[ arg[0] + i_vec ] = arg[2]; 00328 } 00329 00330 /*! 00331 Forward mode sparsity operations for StpvOp and StvvOp 00332 00333 \copydetails sparse_store_op 00334 */ 00335 template <class Vector_set> 00336 inline void forward_sparse_store_op( 00337 OpCode op , 00338 const addr_t* arg , 00339 size_t num_combined , 00340 const size_t* combined , 00341 Vector_set& var_sparsity , 00342 Vector_set& vecad_sparsity ) 00343 { 00344 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00345 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 ); 00346 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00347 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined ); 00348 size_t i_v = combined[ arg[0] - 1 ]; 00349 CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() ); 00350 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < var_sparsity.n_set() ); 00351 00352 vecad_sparsity.binary_union(i_v, i_v, arg[2], var_sparsity); 00353 00354 return; 00355 } 00356 00357 /*! 00358 Reverse mode sparsity operations for StpvOp and StvvOp 00359 00360 This routine is given the sparsity patterns for 00361 G(v[x], y , w , u ... ) and it uses them to compute the 00362 sparsity patterns for 00363 \verbatim 00364 H(y , w , u , ... ) = G[ v[x], y , w , u , ... ] 00365 \endverbatim 00366 00367 <!-- replace sparse_store_op --> 00368 The C++ source code corresponding to this operation is 00369 \verbatim 00370 v[x] = y 00371 \endverbatim 00372 where v is a VecAD<Base> vector, x is an AD<Base> object, 00373 and y is AD<Base> or Base objects. 00374 We define the index corresponding to v[x] by 00375 \verbatim 00376 i_v_x = combined[ arg[0] + i_vec ] 00377 \endverbatim 00378 where i_vec is defined under the heading \a arg[1] below: 00379 00380 \tparam Vector_set 00381 is the type used for vectors of sets. It can be either 00382 \c sparse_pack, \c sparse_set, or \c sparse_list. 00383 00384 \param op 00385 is the code corresponding to this operator; i.e., StpvOp or StvvOp 00386 (only used for error checking). 00387 00388 \param arg 00389 \n 00390 \a arg[0] 00391 is the offset corresponding to this VecAD vector in the combined array. 00392 \n 00393 \n 00394 \a arg[2] 00395 \n 00396 The set with index \a arg[2] in \a var_sparsity 00397 is the sparsity pattern corresponding to y. 00398 (Note that \a arg[2] > 0 because y is a variable.) 00399 00400 \param num_combined 00401 is the total number of elements in the VecAD address array. 00402 00403 \param combined 00404 \a combined [ arg[0] - 1 ] 00405 is the index of the set in \a vecad_sparsity corresponding 00406 to the sparsity pattern for the vector v. 00407 We use the notation i_v below which is defined by 00408 \verbatim 00409 i_v = combined[ \a arg[0] - 1 ] 00410 \endverbatim 00411 00412 \param var_sparsity 00413 The set with index \a arg[2] in \a var_sparsity 00414 is the sparsity pattern for y. 00415 This is an input for forward mode operations. 00416 For reverse mode operations: 00417 The sparsity pattern for v is added to the spartisy pattern for y. 00418 00419 \param vecad_sparsity 00420 The set with index \a i_v in \a vecad_sparsity 00421 is the sparsity pattern for v. 00422 This is an input for reverse mode operations. 00423 For forward mode operations, the sparsity pattern for y is added 00424 to the sparsity pattern for the vector v. 00425 00426 \par Checked Assertions 00427 \li NumArg(op) == 3 00428 \li NumRes(op) == 0 00429 \li 0 < \a arg[0] 00430 \li \a arg[0] < \a num_combined 00431 \li \a arg[2] < \a var_sparsity.n_set() 00432 \li i_v < \a vecad_sparsity.n_set() 00433 <!-- end sparse_store_op --> 00434 */ 00435 template <class Vector_set> 00436 inline void reverse_sparse_jacobian_store_op( 00437 OpCode op , 00438 const addr_t* arg , 00439 size_t num_combined , 00440 const size_t* combined , 00441 Vector_set& var_sparsity , 00442 Vector_set& vecad_sparsity ) 00443 { 00444 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00445 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 ); 00446 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00447 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined ); 00448 size_t i_v = combined[ arg[0] - 1 ]; 00449 CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() ); 00450 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < var_sparsity.n_set() ); 00451 00452 var_sparsity.binary_union(arg[2], arg[2], i_v, vecad_sparsity); 00453 00454 return; 00455 } 00456 00457 /*! 00458 Reverse mode sparsity operations for StpvOp and StvvOp 00459 00460 This routine is given the sparsity patterns for 00461 G(v[x], y , w , u ... ) 00462 and it uses them to compute the sparsity patterns for 00463 \verbatim 00464 H(y , w , u , ... ) = G[ v[x], y , w , u , ... ] 00465 \endverbatim 00466 00467 <!-- replace sparse_store_op --> 00468 The C++ source code corresponding to this operation is 00469 \verbatim 00470 v[x] = y 00471 \endverbatim 00472 where v is a VecAD<Base> vector, x is an AD<Base> object, 00473 and y is AD<Base> or Base objects. 00474 We define the index corresponding to v[x] by 00475 \verbatim 00476 i_v_x = combined[ arg[0] + i_vec ] 00477 \endverbatim 00478 where i_vec is defined under the heading \a arg[1] below: 00479 00480 \tparam Vector_set 00481 is the type used for vectors of sets. It can be either 00482 \c sparse_pack, \c sparse_set, or \c sparse_list. 00483 00484 \param op 00485 is the code corresponding to this operator; i.e., StpvOp or StvvOp 00486 (only used for error checking). 00487 00488 \param arg 00489 \n 00490 \a arg[0] 00491 is the offset corresponding to this VecAD vector in the combined array. 00492 \n 00493 \n 00494 \a arg[2] 00495 \n 00496 The set with index \a arg[2] in \a var_sparsity 00497 is the sparsity pattern corresponding to y. 00498 (Note that \a arg[2] > 0 because y is a variable.) 00499 00500 \param num_combined 00501 is the total number of elements in the VecAD address array. 00502 00503 \param combined 00504 \a combined [ arg[0] - 1 ] 00505 is the index of the set in \a vecad_sparsity corresponding 00506 to the sparsity pattern for the vector v. 00507 We use the notation i_v below which is defined by 00508 \verbatim 00509 i_v = combined[ \a arg[0] - 1 ] 00510 \endverbatim 00511 00512 \param var_sparsity 00513 The set with index \a arg[2] in \a var_sparsity 00514 is the sparsity pattern for y. 00515 This is an input for forward mode operations. 00516 For reverse mode operations: 00517 The sparsity pattern for v is added to the spartisy pattern for y. 00518 00519 \param vecad_sparsity 00520 The set with index \a i_v in \a vecad_sparsity 00521 is the sparsity pattern for v. 00522 This is an input for reverse mode operations. 00523 For forward mode operations, the sparsity pattern for y is added 00524 to the sparsity pattern for the vector v. 00525 00526 \par Checked Assertions 00527 \li NumArg(op) == 3 00528 \li NumRes(op) == 0 00529 \li 0 < \a arg[0] 00530 \li \a arg[0] < \a num_combined 00531 \li \a arg[2] < \a var_sparsity.n_set() 00532 \li i_v < \a vecad_sparsity.n_set() 00533 <!-- end sparse_store_op --> 00534 00535 \param var_jacobian 00536 \a var_jacobian[ \a arg[2] ] 00537 is false (true) if the Jacobian of G with respect to y is always zero 00538 (may be non-zero). 00539 00540 \param vecad_jacobian 00541 \a vecad_jacobian[i_v] 00542 is false (true) if the Jacobian with respect to x is always zero 00543 (may be non-zero). 00544 On input, it corresponds to the function G, 00545 and on output it corresponds to the function H. 00546 */ 00547 template <class Vector_set> 00548 inline void reverse_sparse_hessian_store_op( 00549 OpCode op , 00550 const addr_t* arg , 00551 size_t num_combined , 00552 const size_t* combined , 00553 Vector_set& var_sparsity , 00554 Vector_set& vecad_sparsity , 00555 bool* var_jacobian , 00556 bool* vecad_jacobian ) 00557 { 00558 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00559 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 ); 00560 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00561 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined ); 00562 size_t i_v = combined[ arg[0] - 1 ]; 00563 CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() ); 00564 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < var_sparsity.n_set() ); 00565 00566 var_sparsity.binary_union(arg[2], arg[2], i_v, vecad_sparsity); 00567 00568 var_jacobian[ arg[2] ] |= vecad_jacobian[i_v]; 00569 00570 return; 00571 } 00572 00573 } // END_CPPAD_NAMESPACE 00574 # endif