CppAD: A C++ Algorithmic Differentiation Package  20130918
sqrt_op.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_SQRT_OP_INCLUDED
00003 # define CPPAD_SQRT_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 
00017 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
00018 /*!
00019 \file sqrt_op.hpp
00020 Forward and reverse mode calculations for z = sqrt(x).
00021 */
00022 
00023 
00024 /*!
00025 Compute forward mode Taylor coefficient for result of op = SqrtOp.
00026 
00027 The C++ source code corresponding to this operation is
00028 \verbatim
00029      z = sqrt(x)
00030 \endverbatim
00031 
00032 \copydetails forward_unary1_op
00033 */
00034 template <class Base>
00035 inline void forward_sqrt_op(
00036      size_t p           ,
00037      size_t q           ,
00038      size_t i_z         ,
00039      size_t i_x         ,
00040      size_t nc_taylor   , 
00041      Base*  taylor      )
00042 {    
00043      // check assumptions
00044      CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
00045      CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
00046      CPPAD_ASSERT_UNKNOWN( i_x < i_z );
00047      CPPAD_ASSERT_UNKNOWN( q < nc_taylor );
00048      CPPAD_ASSERT_UNKNOWN( p <= q );
00049 
00050      // Taylor coefficients corresponding to argument and result
00051      Base* x = taylor + i_x * nc_taylor;
00052      Base* z = taylor + i_z * nc_taylor;
00053 
00054      size_t k;
00055      if( p == 0 )
00056      {    z[0] = sqrt( x[0] );
00057           p++;
00058      }
00059      for(size_t j = p; j <= q; j++)
00060      {
00061           CPPAD_ASSERT_KNOWN(
00062                x[0] != Base(0),
00063                "Forward: attempt to take derivatve of square root of zero"
00064           )
00065           z[j] = Base(0);
00066           for(k = 1; k < j; k++)
00067                z[j] -= Base(k) * z[k] * z[j-k];
00068           z[j] /= Base(j);
00069           z[j] += x[j] / Base(2);
00070           z[j] /= z[0];
00071      }
00072 }
00073 
00074 /*!
00075 Compute zero order forward mode Taylor coefficient for result of op = SqrtOp.
00076 
00077 The C++ source code corresponding to this operation is
00078 \verbatim
00079      z = sqrt(x)
00080 \endverbatim
00081 
00082 \copydetails forward_unary1_op_0
00083 */
00084 template <class Base>
00085 inline void forward_sqrt_op_0(
00086      size_t i_z         ,
00087      size_t i_x         ,
00088      size_t nc_taylor   , 
00089      Base*  taylor      )
00090 {
00091      // check assumptions
00092      CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
00093      CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
00094      CPPAD_ASSERT_UNKNOWN( i_x < i_z );
00095      CPPAD_ASSERT_UNKNOWN( 0 < nc_taylor );
00096 
00097      // Taylor coefficients corresponding to argument and result
00098      Base* x = taylor + i_x * nc_taylor;
00099      Base* z = taylor + i_z * nc_taylor;
00100 
00101      z[0] = sqrt( x[0] );
00102 }
00103 /*!
00104 Compute reverse mode partial derivatives for result of op = SqrtOp.
00105 
00106 The C++ source code corresponding to this operation is
00107 \verbatim
00108      z = sqrt(x)
00109 \endverbatim
00110 
00111 \copydetails reverse_unary1_op
00112 */
00113 
00114 template <class Base>
00115 inline void reverse_sqrt_op(
00116      size_t      d            ,
00117      size_t      i_z          ,
00118      size_t      i_x          ,
00119      size_t      nc_taylor    , 
00120      const Base* taylor       ,
00121      size_t      nc_partial   ,
00122      Base*       partial      )
00123 {
00124      // check assumptions
00125      CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
00126      CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
00127      CPPAD_ASSERT_UNKNOWN( i_x < i_z );
00128      CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
00129      CPPAD_ASSERT_UNKNOWN( d < nc_partial );
00130 
00131      // Taylor coefficients and partials corresponding to argument
00132      Base* px       = partial + i_x * nc_partial;
00133 
00134      // Taylor coefficients and partials corresponding to result
00135      const Base* z  = taylor  + i_z * nc_taylor;
00136      Base* pz       = partial + i_z * nc_partial;
00137 
00138      CPPAD_ASSERT_KNOWN(
00139           z[0] != Base(0),
00140           "Reverse: attempt to take derivatve of square root of zero"
00141      )
00142 
00143      // number of indices to access
00144      size_t j = d;
00145      size_t k;
00146      while(j)
00147      {
00148 
00149           // scale partial w.r.t. z[j]
00150           pz[j]   /= z[0];
00151 
00152           pz[0]   -= pz[j] * z[j];
00153           px[j]   += pz[j] / Base(2);
00154           for(k = 1; k < j; k++)
00155                pz[k]   -= pz[j] * z[j-k];
00156           --j;
00157      }
00158      px[0] += pz[0] / (Base(2) * z[0]);
00159 }
00160 
00161 } // END_CPPAD_NAMESPACE
00162 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines