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