CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_SQRT_OP_INCLUDED 00003 # define CPPAD_SQRT_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 00017 CPPAD_BEGIN_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 j , 00037 size_t i_z , 00038 size_t i_x , 00039 size_t nc_taylor , 00040 Base* taylor ) 00041 { 00042 // check assumptions 00043 CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 ); 00044 CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 ); 00045 CPPAD_ASSERT_UNKNOWN( i_x < i_z ); 00046 CPPAD_ASSERT_UNKNOWN( j < nc_taylor ); 00047 00048 // Taylor coefficients corresponding to argument and result 00049 Base* x = taylor + i_x * nc_taylor; 00050 Base* z = taylor + i_z * nc_taylor; 00051 00052 size_t k; 00053 if( j == 0 ) 00054 z[j] = sqrt( x[0] ); 00055 else 00056 { 00057 z[j] = Base(0); 00058 for(k = 1; k < j; k++) 00059 z[j] -= Base(k) * z[k] * z[j-k]; 00060 z[j] /= Base(j); 00061 z[j] += x[j] / Base(2); 00062 z[j] /= z[0]; 00063 } 00064 } 00065 00066 /*! 00067 Compute zero order forward mode Taylor coefficient for result of op = SqrtOp. 00068 00069 The C++ source code corresponding to this operation is 00070 \verbatim 00071 z = sqrt(x) 00072 \endverbatim 00073 00074 \copydetails forward_unary1_op_0 00075 */ 00076 template <class Base> 00077 inline void forward_sqrt_op_0( 00078 size_t i_z , 00079 size_t i_x , 00080 size_t nc_taylor , 00081 Base* taylor ) 00082 { 00083 // check assumptions 00084 CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 ); 00085 CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 ); 00086 CPPAD_ASSERT_UNKNOWN( i_x < i_z ); 00087 CPPAD_ASSERT_UNKNOWN( 0 < nc_taylor ); 00088 00089 // Taylor coefficients corresponding to argument and result 00090 Base* x = taylor + i_x * nc_taylor; 00091 Base* z = taylor + i_z * nc_taylor; 00092 00093 z[0] = sqrt( x[0] ); 00094 } 00095 /*! 00096 Compute reverse mode partial derivatives for result of op = SqrtOp. 00097 00098 The C++ source code corresponding to this operation is 00099 \verbatim 00100 z = sqrt(x) 00101 \endverbatim 00102 00103 \copydetails reverse_unary1_op 00104 */ 00105 00106 template <class Base> 00107 inline void reverse_sqrt_op( 00108 size_t d , 00109 size_t i_z , 00110 size_t i_x , 00111 size_t nc_taylor , 00112 const Base* taylor , 00113 size_t nc_partial , 00114 Base* partial ) 00115 { 00116 // check assumptions 00117 CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 ); 00118 CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 ); 00119 CPPAD_ASSERT_UNKNOWN( i_x < i_z ); 00120 CPPAD_ASSERT_UNKNOWN( d < nc_taylor ); 00121 CPPAD_ASSERT_UNKNOWN( d < nc_partial ); 00122 00123 // Taylor coefficients and partials corresponding to argument 00124 Base* px = partial + i_x * nc_partial; 00125 00126 // Taylor coefficients and partials corresponding to result 00127 const Base* z = taylor + i_z * nc_taylor; 00128 Base* pz = partial + i_z * nc_partial; 00129 00130 // number of indices to access 00131 size_t j = d; 00132 size_t k; 00133 while(j) 00134 { // scale partial w.r.t. z[j] 00135 pz[j] /= z[0]; 00136 00137 pz[0] -= pz[j] * z[j]; 00138 px[j] += pz[j] / Base(2); 00139 for(k = 1; k < j; k++) 00140 pz[k] -= pz[j] * z[j-k]; 00141 --j; 00142 } 00143 px[0] += pz[0] / (Base(2) * z[0]); 00144 } 00145 00146 CPPAD_END_NAMESPACE 00147 # endif