CppAD: A C++ Algorithmic Differentiation Package 20110419
load_op.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_LOAD_OP_INCLUDED
00003 # define CPPAD_LOAD_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 load_op.hpp
00020 Setting a variable so that it corresponds to current value of a VecAD element.
00021 */
00022 
00023 /*!
00024 Zero order forward mode implementation of op = LdpOp.
00025 
00026 \copydetails forward_load_op_0
00027 */
00028 template <class Base>
00029 inline void forward_load_p_op_0(
00030         size_t         i_z         ,
00031         size_t*        arg         , 
00032         size_t         num_par     ,
00033         const Base*    parameter   ,
00034         size_t         nc_taylor   ,
00035         Base*          taylor      ,
00036         size_t         nc_combined ,
00037         const bool*    variable    ,
00038         const size_t*  combined    )
00039 {       size_t i_vec = arg[1];
00040 
00041         // Because the index is a parameter, this indexing error should be
00042         // caught and reported to the user at an when the tape is recording.
00043         CPPAD_ASSERT_UNKNOWN( i_vec < combined[ arg[0] - 1 ] );
00044 
00045 
00046         CPPAD_ASSERT_UNKNOWN( variable != CPPAD_NULL );
00047         CPPAD_ASSERT_UNKNOWN( combined != CPPAD_NULL );
00048         CPPAD_ASSERT_UNKNOWN( NumArg(LdpOp) == 3 );
00049         CPPAD_ASSERT_UNKNOWN( NumRes(LdpOp) == 1 );
00050         CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
00051         CPPAD_ASSERT_UNKNOWN( arg[0] + i_vec < nc_combined );
00052 
00053         size_t combined_index = arg[0] + i_vec;
00054         size_t i_y_x          = combined[ combined_index ];     
00055         Base* z  = taylor + i_z * nc_taylor;
00056         if( variable[ combined_index ] )
00057         {       CPPAD_ASSERT_UNKNOWN( i_y_x < i_z );
00058                 Base* y_x = taylor + i_y_x * nc_taylor;
00059                 arg[2]    = i_y_x;
00060                 z[0]      = y_x[0];
00061         }
00062         else
00063         {       CPPAD_ASSERT_UNKNOWN( i_y_x < num_par );
00064                 Base y_x  = parameter[i_y_x];
00065                 arg[2]    = 0;
00066                 z[0]      = y_x;
00067         }
00068 }
00069 
00070 /*!
00071 Zero order forward mode implementation of op = LdvOp.
00072 
00073 \copydetails forward_load_op_0
00074 */
00075 template <class Base>
00076 inline void forward_load_v_op_0(
00077         size_t         i_z         ,
00078         size_t*        arg         , 
00079         size_t         num_par     ,
00080         const Base*    parameter   ,
00081         size_t         nc_taylor   ,
00082         Base*          taylor      ,
00083         size_t         nc_combined ,
00084         const bool*    variable    ,
00085         const size_t*  combined    )
00086 {
00087         CPPAD_ASSERT_UNKNOWN( variable != CPPAD_NULL );
00088         CPPAD_ASSERT_UNKNOWN( combined != CPPAD_NULL );
00089         CPPAD_ASSERT_UNKNOWN( NumArg(LdvOp) == 3 );
00090         CPPAD_ASSERT_UNKNOWN( NumRes(LdvOp) == 1 );
00091         CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
00092         CPPAD_ASSERT_UNKNOWN( arg[1] < i_z );
00093 
00094         size_t i_vec = Integer( taylor[ arg[1] * nc_taylor + 0 ] );
00095         CPPAD_ASSERT_KNOWN( 
00096                 i_vec < combined[ arg[0] - 1 ] ,
00097                 "VecAD: index during zero order forward sweep is out of range"
00098         );
00099         CPPAD_ASSERT_UNKNOWN( arg[0] + i_vec < nc_combined );
00100 
00101         size_t combined_index = arg[0] + i_vec;
00102         size_t i_y_x          = combined[ combined_index ];     
00103 
00104 
00105         Base* z  = taylor + i_z * nc_taylor;
00106         if( variable[ combined_index ] )
00107         {       CPPAD_ASSERT_UNKNOWN( i_y_x < i_z );
00108                 Base* y_x = taylor + i_y_x * nc_taylor;
00109                 arg[2]    = i_y_x;
00110                 z[0]      = y_x[0];
00111         }
00112         else
00113         {       CPPAD_ASSERT_UNKNOWN( i_y_x < num_par );
00114                 Base y_x  = parameter[i_y_x];
00115                 arg[2]    = 0;
00116                 z[0]      = y_x;
00117         }
00118 }
00119 
00120 /*!
00121 Forward mode, except for zero order, for op = LdpOp or op = LdvOp
00122 
00123 The C++ source code corresponding to this operation is
00124 \verbatim
00125         z = y[x]
00126 \endverbatim
00127 where y is a VecAD<Base> vector and x is an AD<Base> or Base index. 
00128 
00129 \tparam Base
00130 base type for the operator; i.e., this operation was recorded
00131 using AD< \a Base > and computations by this routine are done using type 
00132 \a Base.
00133 
00134 \param op
00135 is the code corresponding to this operator; i.e., LdpOp or LdvOp
00136 (only used for error checking).
00137 
00138 \param d
00139 is the order of the Taylor coefficient that we are computing.
00140 
00141 \param i_z
00142 is the AD variable index corresponding to the variable z.
00143 
00144 \param arg
00145 \a arg[2]
00146 If y[x] is a parameter, \a arg[2] is zero 
00147 (which is not a valid variable index).
00148 If y[x] is a variable, 
00149 \a arg[2] is the variable index corresponding to y[x].
00150 
00151 \param nc_taylor
00152 number of columns in the matrix containing the Taylor coefficients.
00153 
00154 \param taylor
00155 \b Input: if y[x] is a variable, \a taylor[ \a arg[2] * nc_taylor + d ]
00156 is the d-order Taylor coefficient corresponding to y[x].
00157 \n
00158 \b Output: \a taylor[ i_z * nc_taylor + d ]
00159 is the d-order Taylor coefficient for the variable z.
00160 
00161 \par Checked Assertions 
00162 \li NumArg(op) == 3
00163 \li NumRes(op) == 1
00164 \li 0 < d < nc_taylor
00165 \li arg[2] < i_z
00166 */
00167 template <class Base>
00168 inline void forward_load_op(
00169         OpCode         op          ,
00170         size_t         d           ,
00171         size_t         i_z         ,
00172         const size_t*  arg         , 
00173         size_t         nc_taylor   ,
00174         Base*          taylor      )
00175 {
00176 
00177         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
00178         CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 );
00179         CPPAD_ASSERT_UNKNOWN( d > 0 )
00180         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
00181         CPPAD_ASSERT_UNKNOWN( arg[2] < i_z );
00182 
00183         Base* z      = taylor + i_z * nc_taylor;
00184         if( arg[2] > 0 )
00185         {       Base* y_x = taylor + arg[2] * nc_taylor;
00186                 z[d]      = y_x[d];
00187         }
00188         else    z[d]      = Base(0);
00189 }
00190 
00191 /*!
00192 Reverse mode for op = LdpOp or LdvOp.
00193 
00194 The C++ source code corresponding to this operation is
00195 \verbatim
00196         z = y[x]
00197 \endverbatim
00198 where y is a VecAD<Base> vector and x is an AD<Base> or Base index. 
00199 
00200 This routine is given the partial derivatives of a function 
00201 G(z , y[x] , w , u ... )
00202 and it uses them to compute the partial derivatives of 
00203 \verbatim
00204         H( y[x] , w , u , ... ) = G[ z( y[x] ) , y[x] , w , u , ... ]
00205 \endverbatim
00206 
00207 \tparam Base
00208 base type for the operator; i.e., this operation was recorded
00209 using AD< \a Base > and computations by this routine are done using type 
00210 \a Base.
00211 
00212 \param op
00213 is the code corresponding to this operator; i.e., LdpOp or LdvOp
00214 (only used for error checking).
00215 
00216 \param d
00217 highest order the Taylor coefficient that we are computing the partial
00218 derivative with respect to.
00219 
00220 \param i_z
00221 is the AD variable index corresponding to the variable z.
00222 
00223 \param arg
00224 \a arg[2]
00225 If y[x] is a parameter, \a arg[2] is zero 
00226 (which is not a valid variable index).
00227 If y[x] is a variable, 
00228 \a arg[2] is the variable index corresponding to y[x].
00229 
00230 \param nc_taylor
00231 number of columns in the matrix containing the Taylor coefficients
00232 (not used).
00233 
00234 \param taylor
00235 matrix of Taylor coefficients (not used).
00236 
00237 \param nc_partial
00238 number of colums in the matrix containing all the partial derivatives
00239 (not used if \a arg[2] is zero).
00240 
00241 \param partial
00242 If \a arg[2] is zero, y[x] is a parameter
00243 and no values need to be modified; i.e., \a partial is not used.
00244 Otherwise, y[x] is a variable and:
00245 \n
00246 \n
00247 \a partial [ \a i_z * \a nc_partial + k ] 
00248 for k = 0 , ... , \a d
00249 is the partial derivative of G
00250 with respect to the k-th order Taylor coefficient for z.
00251 \n
00252 \n
00253 If \a arg[2] is not zero,
00254 \a partial [ \a arg[2] * \a nc_partial + k ]
00255 for k = 0 , ... , \a d
00256 is the partial derivative with respect to 
00257 the k-th order Taylor coefficient for x.
00258 On input, it corresponds to the function G,
00259 and on output it corresponds to the the function H. 
00260 
00261 \par Checked Assertions 
00262 \li NumArg(op) == 3
00263 \li NumRes(op) == 1
00264 \li d < nc_taylor
00265 \li arg[2] < i_z
00266 */
00267 template <class Base>
00268 inline void reverse_load_op(
00269         OpCode         op          ,
00270         size_t         d           ,
00271         size_t         i_z         ,
00272         const size_t*  arg         , 
00273         size_t         nc_taylor   ,
00274         const Base*    taylor      ,
00275         size_t         nc_partial  ,
00276         Base*          partial     )
00277 {
00278 
00279         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
00280         CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 );
00281         CPPAD_ASSERT_UNKNOWN( d < nc_taylor );
00282         CPPAD_ASSERT_UNKNOWN( arg[2] < i_z );
00283 
00284         if( arg[2] > 0 )
00285         {
00286                 Base* pz   = partial + i_z    * nc_partial;
00287                 Base* py_x = partial + arg[2] * nc_partial;
00288                 size_t j = d + 1;
00289                 while(j--)
00290                         py_x[j]   += pz[j];
00291         }
00292 }
00293 
00294 
00295 /*!
00296 Forward mode sparsity operations for LdpOp and LdvOp
00297 
00298 \copydetails sparse_load_op
00299 */
00300 template <class Vector_set>
00301 inline void forward_sparse_load_op(
00302         OpCode             op             ,
00303         size_t             i_z            ,
00304         const size_t*      arg            , 
00305         size_t             num_combined   ,
00306         const size_t*      combined       ,
00307         Vector_set&        var_sparsity   ,
00308         Vector_set&        vecad_sparsity )
00309 {
00310         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
00311         CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 );
00312         CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
00313         CPPAD_ASSERT_UNKNOWN( arg[0] < num_combined );
00314         size_t i_v = combined[ arg[0] - 1 ];
00315         CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() );
00316 
00317         var_sparsity.assignment(i_z, i_v, vecad_sparsity);
00318 
00319         return;
00320 }
00321 
00322 
00323 /*!
00324 Reverse mode Jacobian sparsity operations for LdpOp and LdvOp
00325 
00326 \copydetails sparse_load_op
00327 */
00328 template <class Vector_set>
00329 inline void reverse_sparse_jacobian_load_op(
00330         OpCode             op             ,
00331         size_t             i_z            ,
00332         const size_t*      arg            , 
00333         size_t             num_combined   ,
00334         const size_t*      combined       ,
00335         Vector_set&        var_sparsity   ,
00336         Vector_set&        vecad_sparsity )
00337 {
00338         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
00339         CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 );
00340         CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
00341         CPPAD_ASSERT_UNKNOWN( arg[0] < num_combined );
00342         size_t i_v = combined[ arg[0] - 1 ];
00343         CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() );
00344 
00345         vecad_sparsity.binary_union(i_v, i_v, i_z, var_sparsity);
00346 
00347         return;
00348 }
00349 
00350 
00351 /*!
00352 Reverse mode Hessian sparsity operations for LdpOp and LdvOp
00353 
00354 This routine is given the sparsity patterns for 
00355 G(z , v[x] , w , u ... )
00356 and it uses them to compute the sparsity patterns for
00357 \verbatim
00358         H( v[x] , w , u , ... ) = G[ z( v[x] ) , v[x] , w , u , ... ]
00359 \endverbatim
00360 
00361 \copydetails sparse_load_op
00362 
00363 \param var_jacobian
00364 \a var_jacobian[i_z] 
00365 is false (true) if the Jacobian of G with respect to z is always zero 
00366 (many be non-zero).
00367 
00368 \param vecad_jacobian
00369 \a vecad_jacobian[i_v] 
00370 is false (true) if the Jacobian with respect to x is always zero 
00371 (may be non-zero).
00372 On input, it corresponds to the function G,
00373 and on output it corresponds to the function H.
00374 
00375 */
00376 template <class Vector_set>
00377 inline void reverse_sparse_hessian_load_op(
00378         OpCode             op             ,
00379         size_t             i_z            ,
00380         const size_t*      arg            , 
00381         size_t             num_combined   ,
00382         const size_t*      combined       ,
00383         Vector_set&        var_sparsity   ,
00384         Vector_set&        vecad_sparsity ,
00385         bool*              var_jacobian   ,
00386         bool*              vecad_jacobian )
00387 {
00388         CPPAD_ASSERT_UNKNOWN( NumArg(op) == 3 );
00389         CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 );
00390         CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
00391         CPPAD_ASSERT_UNKNOWN( arg[0] < num_combined );
00392         size_t i_v = combined[ arg[0] - 1 ];
00393         CPPAD_ASSERT_UNKNOWN( i_v < vecad_sparsity.n_set() );
00394 
00395         vecad_sparsity.binary_union(i_v, i_v, i_z, var_sparsity);
00396 
00397         vecad_jacobian[i_v] |= var_jacobian[i_z];
00398 
00399         return;
00400 }
00401 
00402 
00403 
00404 CPPAD_END_NAMESPACE
00405 # endif