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