CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_FOR_TWO_INCLUDED 00003 # define CPPAD_FOR_TWO_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 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 $begin ForTwo$$ 00018 $spell 00019 ddy 00020 typename 00021 Taylor 00022 const 00023 $$ 00024 00025 00026 $index partial, second order driver$$ 00027 $index second, order partial driver$$ 00028 $index driver, second order partial$$ 00029 00030 $index easy, partial$$ 00031 $index driver, easy partial$$ 00032 $index partial, easy$$ 00033 00034 00035 $section Forward Mode Second Partial Derivative Driver$$ 00036 00037 $head Syntax$$ 00038 $syntax%%ddy% = %f%.ForTwo(%x%, %j%, %k%)%$$ 00039 00040 00041 $head Purpose$$ 00042 We use $latex F : B^n \rightarrow B^m$$ to denote the 00043 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$. 00044 The syntax above sets 00045 $latex \[ 00046 ddy [ i * p + \ell ] 00047 = 00048 \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x) 00049 \] $$ 00050 for $latex i = 0 , \ldots , m-1$$ 00051 and $latex \ell = 0 , \ldots , p$$, 00052 where $latex p$$ is the size of the vectors $italic j$$ and $italic k$$. 00053 00054 $head f$$ 00055 The object $italic f$$ has prototype 00056 $syntax% 00057 ADFun<%Base%> %f% 00058 %$$ 00059 Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$ 00060 (see $xref/ForTwo/ForTwo Uses Forward/ForTwo Uses Forward/$$ below). 00061 00062 $head x$$ 00063 The argument $italic x$$ has prototype 00064 $syntax% 00065 const %VectorBase% &%x% 00066 %$$ 00067 (see $xref/ForTwo/VectorBase/VectorBase/$$ below) 00068 and its size 00069 must be equal to $italic n$$, the dimension of the 00070 $xref/seq_property/Domain/domain/$$ space for $italic f$$. 00071 It specifies 00072 that point at which to evaluate the partial derivatives listed above. 00073 00074 $head j$$ 00075 The argument $italic j$$ has prototype 00076 $syntax% 00077 const %VectorSize_t% &%j% 00078 %$$ 00079 (see $xref/ForTwo/VectorSize_t/VectorSize_t/$$ below) 00080 We use $italic p$$ to denote the size of the vector $italic j$$. 00081 All of the indices in $italic j$$ 00082 must be less than $italic n$$; i.e., 00083 for $latex \ell = 0 , \ldots , p-1$$, $latex j[ \ell ] < n$$. 00084 00085 $head k$$ 00086 The argument $italic k$$ has prototype 00087 $syntax% 00088 const %VectorSize_t% &%k% 00089 %$$ 00090 (see $xref/ForTwo/VectorSize_t/VectorSize_t/$$ below) 00091 and its size must be equal to $italic p$$, 00092 the size of the vector $italic j$$. 00093 All of the indices in $italic k$$ 00094 must be less than $italic n$$; i.e., 00095 for $latex \ell = 0 , \ldots , p-1$$, $latex k[ \ell ] < n$$. 00096 00097 $head ddy$$ 00098 The result $italic ddy$$ has prototype 00099 $syntax% 00100 %VectorBase% %ddy% 00101 %$$ 00102 (see $xref/ForTwo/VectorBase/VectorBase/$$ below) 00103 and its size is $latex m * p$$. 00104 It contains the requested partial derivatives; to be specific, 00105 for $latex i = 0 , \ldots , m - 1 $$ 00106 and $latex \ell = 0 , \ldots , p - 1$$ 00107 $latex \[ 00108 ddy [ i * p + \ell ] 00109 = 00110 \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x) 00111 \] $$ 00112 00113 $head VectorBase$$ 00114 The type $italic VectorBase$$ must be a $xref/SimpleVector/$$ class with 00115 $xref/SimpleVector/Elements of Specified Type/elements of type Base/$$. 00116 The routine $xref/CheckSimpleVector/$$ will generate an error message 00117 if this is not the case. 00118 00119 $head VectorSize_t$$ 00120 The type $italic VectorSize_t$$ must be a $xref/SimpleVector/$$ class with 00121 $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$. 00122 The routine $xref/CheckSimpleVector/$$ will generate an error message 00123 if this is not the case. 00124 00125 $head ForTwo Uses Forward$$ 00126 After each call to $xref/Forward/$$, 00127 the object $italic f$$ contains the corresponding 00128 $xref/glossary/Taylor Coefficient/Taylor coefficients/$$. 00129 After $code ForTwo$$, 00130 the previous calls to $xref/Forward/$$ are undefined. 00131 00132 00133 $head Examples$$ 00134 $children% 00135 example/for_two.cpp 00136 %$$ 00137 The routine 00138 $xref/ForTwo.cpp//ForTwo/$$ is both an example and test. 00139 It returns $code true$$, if it succeeds and $code false$$ otherwise. 00140 00141 $end 00142 ----------------------------------------------------------------------------- 00143 */ 00144 00145 // BEGIN CppAD namespace 00146 namespace CppAD { 00147 00148 template <typename Base> 00149 template <typename VectorBase, typename VectorSize_t> 00150 VectorBase ADFun<Base>::ForTwo( 00151 const VectorBase &x, 00152 const VectorSize_t &j, 00153 const VectorSize_t &k) 00154 { size_t i; 00155 size_t j1; 00156 size_t k1; 00157 size_t l; 00158 00159 size_t n = Domain(); 00160 size_t m = Range(); 00161 size_t p = j.size(); 00162 00163 // check VectorBase is Simple Vector class with Base type elements 00164 CheckSimpleVector<Base, VectorBase>(); 00165 00166 // check VectorSize_t is Simple Vector class with size_t elements 00167 CheckSimpleVector<size_t, VectorSize_t>(); 00168 00169 CPPAD_ASSERT_KNOWN( 00170 x.size() == n, 00171 "ForTwo: Length of x not equal domain dimension for f." 00172 ); 00173 CPPAD_ASSERT_KNOWN( 00174 j.size() == k.size(), 00175 "ForTwo: Lenght of the j and k vectors are not equal." 00176 ); 00177 // point at which we are evaluating the second partials 00178 Forward(0, x); 00179 00180 00181 // dimension the return value 00182 VectorBase ddy(m * p); 00183 00184 // allocate memory to hold all possible diagonal Taylor coefficients 00185 // (for large sparse cases, this is not efficient) 00186 VectorBase D(m * n); 00187 00188 // boolean flag for which diagonal coefficients are computed 00189 CppAD::vector<bool> c(n); 00190 for(j1 = 0; j1 < n; j1++) 00191 c[j1] = false; 00192 00193 // direction vector in argument space 00194 VectorBase dx(n); 00195 for(j1 = 0; j1 < n; j1++) 00196 dx[j1] = Base(0); 00197 00198 // result vector in range space 00199 VectorBase dy(m); 00200 00201 // compute the diagonal coefficients that are needed 00202 for(l = 0; l < p; l++) 00203 { j1 = j[l]; 00204 k1 = k[l]; 00205 CPPAD_ASSERT_KNOWN( 00206 j1 < n, 00207 "ForTwo: an element of j not less than domain dimension for f." 00208 ); 00209 CPPAD_ASSERT_KNOWN( 00210 k1 < n, 00211 "ForTwo: an element of k not less than domain dimension for f." 00212 ); 00213 size_t count = 2; 00214 while(count) 00215 { count--; 00216 if( ! c[j1] ) 00217 { // diagonal term in j1 direction 00218 c[j1] = true; 00219 dx[j1] = Base(1); 00220 Forward(1, dx); 00221 00222 dx[j1] = Base(0); 00223 dy = Forward(2, dx); 00224 for(i = 0; i < m; i++) 00225 D[i * n + j1 ] = dy[i]; 00226 } 00227 j1 = k1; 00228 } 00229 } 00230 // compute all the requested cross partials 00231 for(l = 0; l < p; l++) 00232 { j1 = j[l]; 00233 k1 = k[l]; 00234 if( j1 == k1 ) 00235 { for(i = 0; i < m; i++) 00236 ddy[i * p + l] = Base(2) * D[i * n + j1]; 00237 } 00238 else 00239 { 00240 // cross term in j1 and k1 directions 00241 dx[j1] = Base(1); 00242 dx[k1] = Base(1); 00243 Forward(1, dx); 00244 00245 dx[j1] = Base(0); 00246 dx[k1] = Base(0); 00247 dy = Forward(2, dx); 00248 00249 // place result in return value 00250 for(i = 0; i < m; i++) 00251 ddy[i * p + l] = dy[i] - D[i*n+j1] - D[i*n+k1]; 00252 00253 } 00254 } 00255 return ddy; 00256 } 00257 00258 } // END CppAD namespace 00259 00260 # endif