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