CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_DIV_OP_INCLUDED 00003 # define CPPAD_DIV_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 div_op.hpp 00019 Forward and reverse mode calculations for z = x / y. 00020 */ 00021 00022 // --------------------------- Divvv ----------------------------------------- 00023 /*! 00024 Compute forward mode Taylor coefficients for result of op = DivvvOp. 00025 00026 The C++ source code corresponding to this operation is 00027 \verbatim 00028 z = x / y 00029 \endverbatim 00030 In the documentation below, 00031 this operations is for the case where both x and y are variables 00032 and the argument \a parameter is not used. 00033 00034 \copydetails forward_binary_op 00035 */ 00036 00037 template <class Base> 00038 inline void forward_divvv_op( 00039 size_t p , 00040 size_t q , 00041 size_t i_z , 00042 const addr_t* arg , 00043 const Base* parameter , 00044 size_t nc_taylor , 00045 Base* taylor ) 00046 { 00047 // check assumptions 00048 CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 ); 00049 CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 ); 00050 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z ); 00051 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z ); 00052 CPPAD_ASSERT_UNKNOWN( q < nc_taylor ); 00053 CPPAD_ASSERT_UNKNOWN( p <= q ); 00054 00055 // Taylor coefficients corresponding to arguments and result 00056 Base* x = taylor + arg[0] * nc_taylor; 00057 Base* y = taylor + arg[1] * nc_taylor; 00058 Base* z = taylor + i_z * nc_taylor; 00059 00060 00061 // Using CondExp, it can make sense to divide by zero, 00062 // so do not make it an error. 00063 size_t k; 00064 for(size_t d = p; d <= q; d++) 00065 { z[d] = x[d]; 00066 for(k = 1; k <= d; k++) 00067 z[d] -= z[d-k] * y[k]; 00068 z[d] /= y[0]; 00069 } 00070 } 00071 00072 00073 /*! 00074 Compute zero order forward mode Taylor coefficients for result of op = DivvvOp. 00075 00076 The C++ source code corresponding to this operation is 00077 \verbatim 00078 z = x / y 00079 \endverbatim 00080 In the documentation below, 00081 this operations is for the case where both x and y are variables 00082 and the argument \a parameter is not used. 00083 00084 \copydetails forward_binary_op_0 00085 */ 00086 00087 template <class Base> 00088 inline void forward_divvv_op_0( 00089 size_t i_z , 00090 const addr_t* arg , 00091 const Base* parameter , 00092 size_t nc_taylor , 00093 Base* taylor ) 00094 { 00095 // check assumptions 00096 CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 ); 00097 CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 ); 00098 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z ); 00099 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z ); 00100 00101 // Taylor coefficients corresponding to arguments and result 00102 Base* x = taylor + arg[0] * nc_taylor; 00103 Base* y = taylor + arg[1] * nc_taylor; 00104 Base* z = taylor + i_z * nc_taylor; 00105 00106 z[0] = x[0] / y[0]; 00107 } 00108 00109 /*! 00110 Compute reverse mode partial derivatives for result of op = DivvvOp. 00111 00112 The C++ source code corresponding to this operation is 00113 \verbatim 00114 z = x / y 00115 \endverbatim 00116 In the documentation below, 00117 this operations is for the case where both x and y are variables 00118 and the argument \a parameter is not used. 00119 00120 \copydetails reverse_binary_op 00121 */ 00122 00123 template <class Base> 00124 inline void reverse_divvv_op( 00125 size_t d , 00126 size_t i_z , 00127 const addr_t* arg , 00128 const Base* parameter , 00129 size_t nc_taylor , 00130 const Base* taylor , 00131 size_t nc_partial , 00132 Base* partial ) 00133 { 00134 // check assumptions 00135 CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 ); 00136 CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 ); 00137 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z ); 00138 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z ); 00139 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00140 CPPAD_ASSERT_UNKNOWN( d < nc_partial ); 00141 00142 // Arguments 00143 const Base* y = taylor + arg[1] * nc_taylor; 00144 const Base* z = taylor + i_z * nc_taylor; 00145 00146 // Partial derivatives corresponding to arguments and result 00147 Base* px = partial + arg[0] * nc_partial; 00148 Base* py = partial + arg[1] * nc_partial; 00149 Base* pz = partial + i_z * nc_partial; 00150 00151 // Using CondExp, it can make sense to divide by zero 00152 // so do not make it an error. 00153 00154 size_t k; 00155 // number of indices to access 00156 size_t j = d + 1; 00157 while(j) 00158 { --j; 00159 // scale partial w.r.t. z[j] 00160 pz[j] /= y[0]; 00161 00162 px[j] += pz[j]; 00163 for(k = 1; k <= j; k++) 00164 { pz[j-k] -= pz[j] * y[k]; 00165 py[k] -= pz[j] * z[j-k]; 00166 } 00167 py[0] -= pz[j] * z[j]; 00168 } 00169 } 00170 00171 // --------------------------- Divpv ----------------------------------------- 00172 /*! 00173 Compute forward mode Taylor coefficients for result of op = DivpvOp. 00174 00175 The C++ source code corresponding to this operation is 00176 \verbatim 00177 z = x / y 00178 \endverbatim 00179 In the documentation below, 00180 this operations is for the case where x is a parameter and y is a variable. 00181 00182 \copydetails forward_binary_op 00183 */ 00184 00185 template <class Base> 00186 inline void forward_divpv_op( 00187 size_t p , 00188 size_t q , 00189 size_t i_z , 00190 const addr_t* arg , 00191 const Base* parameter , 00192 size_t nc_taylor , 00193 Base* taylor ) 00194 { 00195 // check assumptions 00196 CPPAD_ASSERT_UNKNOWN( NumArg(DivpvOp) == 2 ); 00197 CPPAD_ASSERT_UNKNOWN( NumRes(DivpvOp) == 1 ); 00198 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z ); 00199 CPPAD_ASSERT_UNKNOWN( q < nc_taylor ); 00200 CPPAD_ASSERT_UNKNOWN( p <= q ); 00201 00202 // Taylor coefficients corresponding to arguments and result 00203 Base* y = taylor + arg[1] * nc_taylor; 00204 Base* z = taylor + i_z * nc_taylor; 00205 00206 // Paraemter value 00207 Base x = parameter[ arg[0] ]; 00208 00209 // Using CondExp, it can make sense to divide by zero, 00210 // so do not make it an error. 00211 size_t k; 00212 if( p == 0 ) 00213 { z[0] = x / y[0]; 00214 p++; 00215 } 00216 for(size_t d = p; d <= q; d++) 00217 { z[d] = Base(0); 00218 for(k = 1; k <= d; k++) 00219 z[d] -= z[d-k] * y[k]; 00220 z[d] /= y[0]; 00221 } 00222 } 00223 00224 /*! 00225 Compute zero order forward mode Taylor coefficient for result of op = DivpvOp. 00226 00227 The C++ source code corresponding to this operation is 00228 \verbatim 00229 z = x / y 00230 \endverbatim 00231 In the documentation below, 00232 this operations is for the case where x is a parameter and y is a variable. 00233 00234 \copydetails forward_binary_op_0 00235 */ 00236 00237 template <class Base> 00238 inline void forward_divpv_op_0( 00239 size_t i_z , 00240 const addr_t* arg , 00241 const Base* parameter , 00242 size_t nc_taylor , 00243 Base* taylor ) 00244 { 00245 // check assumptions 00246 CPPAD_ASSERT_UNKNOWN( NumArg(DivpvOp) == 2 ); 00247 CPPAD_ASSERT_UNKNOWN( NumRes(DivpvOp) == 1 ); 00248 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z ); 00249 00250 // Paraemter value 00251 Base x = parameter[ arg[0] ]; 00252 00253 // Taylor coefficients corresponding to arguments and result 00254 Base* y = taylor + arg[1] * nc_taylor; 00255 Base* z = taylor + i_z * nc_taylor; 00256 00257 z[0] = x / y[0]; 00258 } 00259 00260 /*! 00261 Compute reverse mode partial derivative for result of op = DivpvOp. 00262 00263 The C++ source code corresponding to this operation is 00264 \verbatim 00265 z = x / y 00266 \endverbatim 00267 In the documentation below, 00268 this operations is for the case where x is a parameter and y is a variable. 00269 00270 \copydetails reverse_binary_op 00271 */ 00272 00273 template <class Base> 00274 inline void reverse_divpv_op( 00275 size_t d , 00276 size_t i_z , 00277 const addr_t* arg , 00278 const Base* parameter , 00279 size_t nc_taylor , 00280 const Base* taylor , 00281 size_t nc_partial , 00282 Base* partial ) 00283 { 00284 // check assumptions 00285 CPPAD_ASSERT_UNKNOWN( NumArg(DivvvOp) == 2 ); 00286 CPPAD_ASSERT_UNKNOWN( NumRes(DivvvOp) == 1 ); 00287 CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < i_z ); 00288 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00289 CPPAD_ASSERT_UNKNOWN( d < nc_partial ); 00290 00291 // Arguments 00292 const Base* y = taylor + arg[1] * nc_taylor; 00293 const Base* z = taylor + i_z * nc_taylor; 00294 00295 // Partial derivatives corresponding to arguments and result 00296 Base* py = partial + arg[1] * nc_partial; 00297 Base* pz = partial + i_z * nc_partial; 00298 00299 // Using CondExp, it can make sense to divide by zero so do not 00300 // make it an error. 00301 00302 size_t k; 00303 // number of indices to access 00304 size_t j = d + 1; 00305 while(j) 00306 { --j; 00307 // scale partial w.r.t z[j] 00308 pz[j] /= y[0]; 00309 00310 for(k = 1; k <= j; k++) 00311 { pz[j-k] -= pz[j] * y[k]; 00312 py[k] -= pz[j] * z[j-k]; 00313 } 00314 py[0] -= pz[j] * z[j]; 00315 } 00316 } 00317 00318 00319 // --------------------------- Divvp ----------------------------------------- 00320 /*! 00321 Compute forward mode Taylor coefficients for result of op = DivvvOp. 00322 00323 The C++ source code corresponding to this operation is 00324 \verbatim 00325 z = x / y 00326 \endverbatim 00327 In the documentation below, 00328 this operations is for the case where x is a variable and y is a parameter. 00329 00330 \copydetails forward_binary_op 00331 */ 00332 00333 template <class Base> 00334 inline void forward_divvp_op( 00335 size_t p , 00336 size_t q , 00337 size_t i_z , 00338 const addr_t* arg , 00339 const Base* parameter , 00340 size_t nc_taylor , 00341 Base* taylor ) 00342 { 00343 // check assumptions 00344 CPPAD_ASSERT_UNKNOWN( NumArg(DivvpOp) == 2 ); 00345 CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 ); 00346 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z ); 00347 CPPAD_ASSERT_UNKNOWN( q < nc_taylor ); 00348 CPPAD_ASSERT_UNKNOWN( p <= q ); 00349 00350 // Taylor coefficients corresponding to arguments and result 00351 Base* x = taylor + arg[0] * nc_taylor; 00352 Base* z = taylor + i_z * nc_taylor; 00353 00354 // Parameter value 00355 Base y = parameter[ arg[1] ]; 00356 00357 // Using CondExp and multiple levels of AD, it can make sense 00358 // to divide by zero so do not make it an error. 00359 for(size_t d = p; d <= q; d++) 00360 z[d] = x[d] / y; 00361 } 00362 00363 00364 /*! 00365 Compute zero order forward mode Taylor coefficients for result of op = DivvvOp. 00366 00367 The C++ source code corresponding to this operation is 00368 \verbatim 00369 z = x / y 00370 \endverbatim 00371 In the documentation below, 00372 this operations is for the case where x is a variable and y is a parameter. 00373 00374 \copydetails forward_binary_op_0 00375 */ 00376 00377 template <class Base> 00378 inline void forward_divvp_op_0( 00379 size_t i_z , 00380 const addr_t* arg , 00381 const Base* parameter , 00382 size_t nc_taylor , 00383 Base* taylor ) 00384 { 00385 // check assumptions 00386 CPPAD_ASSERT_UNKNOWN( NumArg(DivvpOp) == 2 ); 00387 CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 ); 00388 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z ); 00389 00390 // Parameter value 00391 Base y = parameter[ arg[1] ]; 00392 00393 // Taylor coefficients corresponding to arguments and result 00394 Base* x = taylor + arg[0] * nc_taylor; 00395 Base* z = taylor + i_z * nc_taylor; 00396 00397 z[0] = x[0] / y; 00398 } 00399 00400 /*! 00401 Compute reverse mode partial derivative for result of op = DivvpOp. 00402 00403 The C++ source code corresponding to this operation is 00404 \verbatim 00405 z = x / y 00406 \endverbatim 00407 In the documentation below, 00408 this operations is for the case where x is a variable and y is a parameter. 00409 00410 \copydetails reverse_binary_op 00411 */ 00412 00413 template <class Base> 00414 inline void reverse_divvp_op( 00415 size_t d , 00416 size_t i_z , 00417 const addr_t* arg , 00418 const Base* parameter , 00419 size_t nc_taylor , 00420 const Base* taylor , 00421 size_t nc_partial , 00422 Base* partial ) 00423 { 00424 // check assumptions 00425 CPPAD_ASSERT_UNKNOWN( NumArg(DivvpOp) == 2 ); 00426 CPPAD_ASSERT_UNKNOWN( NumRes(DivvpOp) == 1 ); 00427 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < i_z ); 00428 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00429 CPPAD_ASSERT_UNKNOWN( d < nc_partial ); 00430 00431 // Argument values 00432 Base y = parameter[ arg[1] ]; 00433 00434 // Partial derivatives corresponding to arguments and result 00435 Base* px = partial + arg[0] * nc_partial; 00436 Base* pz = partial + i_z * nc_partial; 00437 00438 // Using CondExp, it can make sense to divide by zero 00439 // so do not make it an error. 00440 00441 // number of indices to access 00442 size_t j = d + 1; 00443 while(j) 00444 { --j; 00445 px[j] += pz[j] / y; 00446 } 00447 } 00448 00449 } // END_CPPAD_NAMESPACE 00450 # endif