CppAD: A C++ Algorithmic Differentiation Package 20110419
for_two.hpp
Go to the documentation of this file.
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