CppAD: A C++ Algorithmic Differentiation Package
20130102
|
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