CppAD: A C++ Algorithmic Differentiation Package
20130102
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_LOAD_OP_INCLUDED 00003 # define CPPAD_LOAD_OP_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 CPPAD_BEGIN_NAMESPACE 00018 /*! 00019 \defgroup load_op_hpp load_op.hpp 00020 \{ 00021 \file load_op.hpp 00022 Setting a variable so that it corresponds to current value of a VecAD element. 00023 */ 00024 00025 /*! 00026 Zero order forward mode implementation of op = LdpOp. 00027 00028 \copydetails forward_load_op_0 00029 */ 00030 template <class Base> 00031 inline void forward_load_p_op_0( 00032 size_t i_z , 00033 addr_t* arg , 00034 size_t num_par , 00035 const Base* parameter , 00036 size_t nc_taylor , 00037 Base* taylor , 00038 size_t nc_combined , 00039 const bool* variable , 00040 const size_t* combined ) 00041 { size_t i_vec = arg[1]; 00042 00043 // Because the index is a parameter, this indexing error should be 00044 // caught and reported to the user at an when the tape is recording. 00045 CPPAD_ASSERT_UNKNOWN( i_vec < combined[ arg[0] - 1 ] ); 00046 00047 00048 CPPAD_ASSERT_UNKNOWN( variable != CPPAD_NULL ); 00049 CPPAD_ASSERT_UNKNOWN( combined != CPPAD_NULL ); 00050 CPPAD_ASSERT_UNKNOWN( NumArg(LdpOp) == 3 ); 00051 CPPAD_ASSERT_UNKNOWN( NumRes(LdpOp) == 1 ); 00052 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00053 CPPAD_ASSERT_UNKNOWN( arg[0] + i_vec < nc_combined ); 00054 00055 size_t combined_index = arg[0] + i_vec; 00056 size_t i_y_x = combined[ combined_index ]; 00057 Base* z = taylor + i_z * nc_taylor; 00058 if( variable[ combined_index ] ) 00059 { CPPAD_ASSERT_UNKNOWN( i_y_x < i_z ); 00060 Base* y_x = taylor + i_y_x * nc_taylor; 00061 arg[2] = i_y_x; 00062 z[0] = y_x[0]; 00063 } 00064 else 00065 { CPPAD_ASSERT_UNKNOWN( i_y_x < num_par ); 00066 Base y_x = parameter[i_y_x]; 00067 arg[2] = 0; 00068 z[0] = y_x; 00069 } 00070 } 00071 00072 /*! 00073 Zero order forward mode implementation of op = LdvOp. 00074 00075 \copydetails forward_load_op_0 00076 */ 00077 template <class Base> 00078 inline void forward_load_v_op_0( 00079 size_t i_z , 00080 addr_t* arg , 00081 size_t num_par , 00082 const Base* parameter , 00083 size_t nc_taylor , 00084 Base* taylor , 00085 size_t nc_combined , 00086 const bool* variable , 00087 const size_t* combined ) 00088 { 00089 CPPAD_ASSERT_UNKNOWN( variable != CPPAD_NULL ); 00090 CPPAD_ASSERT_UNKNOWN( combined != CPPAD_NULL ); 00091 CPPAD_ASSERT_UNKNOWN( NumArg(LdvOp) == 3 ); 00092 CPPAD_ASSERT_UNKNOWN( NumRes(LdvOp) == 1 ); 00093 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00094 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z ); 00095 00096 size_t i_vec = Integer( taylor[ arg[1] * nc_taylor + 0 ] ); 00097 CPPAD_ASSERT_KNOWN( 00098 i_vec < combined[ arg[0] - 1 ] , 00099 "VecAD: index during zero order forward sweep is out of range" 00100 ); 00101 CPPAD_ASSERT_UNKNOWN( arg[0] + i_vec < nc_combined ); 00102 00103 size_t combined_index = arg[0] + i_vec; 00104 size_t i_y_x = combined[ combined_index ]; 00105 00106 00107 Base* z = taylor + i_z * nc_taylor; 00108 if( variable[ combined_index ] ) 00109 { CPPAD_ASSERT_UNKNOWN( i_y_x < i_z ); 00110 Base* y_x = taylor + i_y_x * nc_taylor; 00111 arg[2] = i_y_x; 00112 z[0] = y_x[0]; 00113 } 00114 else 00115 { CPPAD_ASSERT_UNKNOWN( i_y_x < num_par ); 00116 Base y_x = parameter[i_y_x]; 00117 arg[2] = 0; 00118 z[0] = y_x; 00119 } 00120 } 00121 00122 /*! 00123 Forward mode, except for zero order, for op = LdpOp or op = LdvOp 00124 00125 The C++ source code corresponding to this operation is 00126 \verbatim 00127 z = y[x] 00128 \endverbatim 00129 where y is a VecAD<Base> vector and x is an AD<Base> or Base index. 00130 00131 \tparam Base 00132 base type for the operator; i.e., this operation was recorded 00133 using AD< \a Base > and computations by this routine are done using type 00134 \a Base. 00135 00136 \param op 00137 is the code corresponding to this operator; i.e., LdpOp or LdvOp 00138 (only used for error checking). 00139 00140 \param d 00141 is the order of the Taylor coefficient that we are computing. 00142 00143 \param i_z 00144 is the AD variable index corresponding to the variable z. 00145 00146 \param arg 00147 \a arg[2] 00148 If y[x] is a parameter, \a arg[2] is zero 00149 (which is not a valid variable index). 00150 If y[x] is a variable, 00151 \a arg[2] is the variable index corresponding to y[x]. 00152 00153 \param nc_taylor 00154 number of columns in the matrix containing the Taylor coefficients. 00155 00156 \param taylor 00157 \b Input: if y[x] is a variable, \a taylor[ \a arg[2] * nc_taylor + d ] 00158 is the d-order Taylor coefficient corresponding to y[x]. 00159 \n 00160 \b Output: \a taylor[ i_z * nc_taylor + d ] 00161 is the d-order Taylor coefficient for the variable z. 00162 00163 \par Checked Assertions 00164 \li NumArg(op) == 3 00165 \li NumRes(op) == 1 00166 \li 0 < d < nc_taylor 00167 \li size_t(arg[2]) < i_z 00168 */ 00169 template <class Base> 00170 inline void forward_load_op( 00171 OpCode op , 00172 size_t d , 00173 size_t i_z , 00174 const addr_t* arg , 00175 size_t nc_taylor , 00176 Base* taylor ) 00177 { 00178 00179 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00180 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 ); 00181 CPPAD_ASSERT_UNKNOWN( d > 0 ) 00182 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00183 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z ); 00184 00185 Base* z = taylor + i_z * nc_taylor; 00186 if( arg[2] > 0 ) 00187 { Base* y_x = taylor + arg[2] * nc_taylor; 00188 z[d] = y_x[d]; 00189 } 00190 else z[d] = Base(0); 00191 } 00192 00193 /*! 00194 Reverse mode for op = LdpOp or LdvOp. 00195 00196 The C++ source code corresponding to this operation is 00197 \verbatim 00198 z = y[x] 00199 \endverbatim 00200 where y is a VecAD<Base> vector and x is an AD<Base> or Base index. 00201 00202 This routine is given the partial derivatives of a function 00203 G(z , y[x] , w , u ... ) 00204 and it uses them to compute the partial derivatives of 00205 \verbatim 00206 H( y[x] , w , u , ... ) = G[ z( y[x] ) , y[x] , w , u , ... ] 00207 \endverbatim 00208 00209 \tparam Base 00210 base type for the operator; i.e., this operation was recorded 00211 using AD< \a Base > and computations by this routine are done using type 00212 \a Base. 00213 00214 \param op 00215 is the code corresponding to this operator; i.e., LdpOp or LdvOp 00216 (only used for error checking). 00217 00218 \param d 00219 highest order the Taylor coefficient that we are computing the partial 00220 derivative with respect to. 00221 00222 \param i_z 00223 is the AD variable index corresponding to the variable z. 00224 00225 \param arg 00226 \a arg[2] 00227 If y[x] is a parameter, \a arg[2] is zero 00228 (which is not a valid variable index). 00229 If y[x] is a variable, 00230 \a arg[2] is the variable index corresponding to y[x]. 00231 00232 \param nc_taylor 00233 number of columns in the matrix containing the Taylor coefficients 00234 (not used). 00235 00236 \param taylor 00237 matrix of Taylor coefficients (not used). 00238 00239 \param nc_partial 00240 number of colums in the matrix containing all the partial derivatives 00241 (not used if \a arg[2] is zero). 00242 00243 \param partial 00244 If \a arg[2] is zero, y[x] is a parameter 00245 and no values need to be modified; i.e., \a partial is not used. 00246 Otherwise, y[x] is a variable and: 00247 \n 00248 \n 00249 \a partial [ \a i_z * \a nc_partial + k ] 00250 for k = 0 , ... , \a d 00251 is the partial derivative of G 00252 with respect to the k-th order Taylor coefficient for z. 00253 \n 00254 \n 00255 If \a arg[2] is not zero, 00256 \a partial [ \a arg[2] * \a nc_partial + k ] 00257 for k = 0 , ... , \a d 00258 is the partial derivative with respect to 00259 the k-th order Taylor coefficient for x. 00260 On input, it corresponds to the function G, 00261 and on output it corresponds to the the function H. 00262 00263 \par Checked Assertions 00264 \li NumArg(op) == 3 00265 \li NumRes(op) == 1 00266 \li d < nc_taylor 00267 \li size_t(arg[2]) < i_z 00268 */ 00269 template <class Base> 00270 inline void reverse_load_op( 00271 OpCode op , 00272 size_t d , 00273 size_t i_z , 00274 const addr_t* arg , 00275 size_t nc_taylor , 00276 const Base* taylor , 00277 size_t nc_partial , 00278 Base* partial ) 00279 { 00280 00281 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00282 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 ); 00283 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00284 CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z ); 00285 00286 if( arg[2] > 0 ) 00287 { 00288 Base* pz = partial + i_z * nc_partial; 00289 Base* py_x = partial + arg[2] * nc_partial; 00290 size_t j = d + 1; 00291 while(j--) 00292 py_x[j] += pz[j]; 00293 } 00294 } 00295 00296 00297 /*! 00298 Forward mode sparsity operations for LdpOp and LdvOp 00299 00300 \copydetails sparse_load_op 00301 */ 00302 template <class Vector_set> 00303 inline void forward_sparse_load_op( 00304 OpCode op , 00305 size_t i_z , 00306 const addr_t* arg , 00307 size_t num_combined , 00308 const size_t* combined , 00309 Vector_set& var_sparsity , 00310 Vector_set& vecad_sparsity ) 00311 { 00312 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00313 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 ); 00314 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00315 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined ); 00316 size_t i_v = combined[ arg[0] - 1 ]; 00317 CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() ); 00318 00319 var_sparsity.assignment(i_z, i_v, vecad_sparsity); 00320 00321 return; 00322 } 00323 00324 00325 /*! 00326 Reverse mode Jacobian sparsity operations for LdpOp and LdvOp 00327 00328 \copydetails sparse_load_op 00329 */ 00330 template <class Vector_set> 00331 inline void reverse_sparse_jacobian_load_op( 00332 OpCode op , 00333 size_t i_z , 00334 const addr_t* arg , 00335 size_t num_combined , 00336 const size_t* combined , 00337 Vector_set& var_sparsity , 00338 Vector_set& vecad_sparsity ) 00339 { 00340 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00341 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 ); 00342 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00343 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined ); 00344 size_t i_v = combined[ arg[0] - 1 ]; 00345 CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() ); 00346 00347 vecad_sparsity.binary_union(i_v, i_v, i_z, var_sparsity); 00348 00349 return; 00350 } 00351 00352 00353 /*! 00354 Reverse mode Hessian sparsity operations for LdpOp and LdvOp 00355 00356 This routine is given the sparsity patterns for 00357 G(z , v[x] , w , u ... ) 00358 and it uses them to compute the sparsity patterns for 00359 \verbatim 00360 H( v[x] , w , u , ... ) = G[ z( v[x] ) , v[x] , w , u , ... ] 00361 \endverbatim 00362 00363 \copydetails sparse_load_op 00364 00365 \param var_jacobian 00366 \a var_jacobian[i_z] 00367 is false (true) if the Jacobian of G with respect to z is always zero 00368 (many be non-zero). 00369 00370 \param vecad_jacobian 00371 \a vecad_jacobian[i_v] 00372 is false (true) if the Jacobian with respect to x is always zero 00373 (may be non-zero). 00374 On input, it corresponds to the function G, 00375 and on output it corresponds to the function H. 00376 00377 */ 00378 template <class Vector_set> 00379 inline void reverse_sparse_hessian_load_op( 00380 OpCode op , 00381 size_t i_z , 00382 const addr_t* arg , 00383 size_t num_combined , 00384 const size_t* combined , 00385 Vector_set& var_sparsity , 00386 Vector_set& vecad_sparsity , 00387 bool* var_jacobian , 00388 bool* vecad_jacobian ) 00389 { 00390 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 ); 00391 CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 ); 00392 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00393 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_combined ); 00394 size_t i_v = combined[ arg[0] - 1 ]; 00395 CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() ); 00396 00397 vecad_sparsity.binary_union(i_v, i_v, i_z, var_sparsity); 00398 00399 vecad_jacobian[i_v] |= var_jacobian[i_z]; 00400 00401 return; 00402 } 00403 00404 00405 00406 /*! \} */ 00407 CPPAD_END_NAMESPACE 00408 # endif