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