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