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