CppAD: A C++ Algorithmic Differentiation Package 20110419
forward0sweep.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_FORWARD0SWEEP_INCLUDED
00003 # define CPPAD_FORWARD0SWEEP_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 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 CPPAD_BEGIN_NAMESPACE
00017 /*!
00018 \file forward0sweep.hpp
00019 Compute zero order forward mode Taylor coefficients.
00020 */
00021 
00022 /*!
00023 \def CPPAD_FORWARD0SWEEP_TRACE
00024 This value is either zero or one. 
00025 Zero is the normal operational value.
00026 If it is one, a trace of every forward0sweep computation is printed.
00027 (Note that forward0sweep is not used if CPPAD_USE_FORWARD0SWEEP is zero).
00028 */
00029 # define CPPAD_FORWARD0SWEEP_TRACE 0
00030 
00031 /*!
00032 Compute zero order forward mode Taylor coefficients.
00033 
00034 \tparam Base
00035 base type for the operator; i.e., this operation was recorded
00036 using AD< \a Base > and computations by this routine are done using type 
00037 \a Base.
00038 
00039 \param print
00040 If \a print is false,
00041 suppress the output that is otherwise generated 
00042 by the PripOp and PrivOp instructions.
00043 
00044 \param n
00045 is the number of independent variables on the tape.
00046 
00047 \param numvar
00048 is the total number of variables on the tape.
00049 This is also equal to the number of rows in the matrix \a Taylor; i.e.,
00050 \a Rec->num_rec_var().
00051 
00052 \param Rec
00053 The information stored in \a Rec
00054 is a recording of the operations corresponding to the function
00055 \f[
00056         F : {\bf R}^n \rightarrow {\bf R}^m
00057 \f]
00058 where \f$ n \f$ is the number of independent variables and
00059 \f$ m \f$ is the number of dependent variables.
00060 \n
00061 \n
00062 The object \a Rec is effectly constant.
00063 There are two exceptions to this.
00064 The first exception is that while palying back the tape
00065 the object \a Rec holds information about the current location
00066 with in the tape and this changes during palyback. 
00067 The second exception is the fact that the 
00068 zero order ( \a d = 0 ) versions of the VecAD operators LdpOp and LdvOp 
00069 modify the corresponding \a op_arg values returned by 
00070 \ref player::next_forward and \ref player::next_reverse; see the
00071 \link load_op.hpp LdpOp and LdvOp \endlink operations.
00072 
00073 \param J
00074 Is the number of columns in the coefficient matrix \a Taylor.
00075 This must be greater than or equal one.
00076 
00077 \param Taylor
00078 \b Input: For j = 1 , ... , \a n,
00079 \a Taylor [ j * J + 0 ]
00080 variable with index i on the tape 
00081 (independent variable with index (j-1) in the independent variable vector).
00082 \n
00083 \n
00084 \b Output: For i = \a n + 1, ... , \a numvar - 1,
00085 \a Taylor [ i * J + 0 ]
00086 is the zero order Taylor coefficient for the variable with 
00087 index i on the tape.
00088 
00089 \a return
00090 The return value is equal to the number of ComOp operations
00091 that have a different result from when the information in 
00092 \a Rec was recorded.
00093 (Note that if NDEBUG is true, there are no ComOp operations
00094 in Rec and hence this return value is always zero.)
00095 */
00096 
00097 template <class Base>
00098 size_t forward0sweep(
00099         bool                  print,
00100         size_t                n,
00101         size_t                numvar,
00102         player<Base>         *Rec,
00103         size_t                J,
00104         Base                 *Taylor
00105 )
00106 {       CPPAD_ASSERT_UNKNOWN( J >= 1 );
00107 
00108         // op code for current instruction
00109         OpCode           op;
00110 
00111         // index for current instruction
00112         size_t         i_op;
00113 
00114         // next variables 
00115         size_t        i_var;
00116 
00117         // constant and non-constant version of the operation argument indices
00118         size_t*         non_const_arg;
00119         const size_t   *arg = 0;
00120 
00121         // initialize the comparision operator (ComOp) counter
00122         size_t compareCount = 0;
00123 
00124         // This is an order zero calculation, initialize vector indices
00125         size_t *VectorInd = CPPAD_NULL;  // address for each element
00126         bool   *VectorVar = CPPAD_NULL;  // is element a variable
00127         size_t  i = Rec->num_rec_vecad_ind();
00128         if( i > 0 )
00129         {       VectorInd = CPPAD_TRACK_NEW_VEC(i, VectorInd);
00130                 VectorVar = CPPAD_TRACK_NEW_VEC(i, VectorVar);
00131                 while(i--)
00132                 {       VectorInd[i] = Rec->GetVecInd(i);
00133                         VectorVar[i] = false;
00134                 }
00135         }
00136 
00137         // work space used by UserOp.
00138         const size_t user_k = 0;     // order of this forward mode calculation
00139         vector<Base> user_tx;        // argument vector Taylor coefficients 
00140         vector<Base> user_ty;        // result vector Taylor coefficients 
00141         size_t user_index = 0;       // indentifier for this user_atomic operation
00142         size_t user_id    = 0;       // user identifier for this call to operator
00143         size_t user_i     = 0;       // index in result vector
00144         size_t user_j     = 0;       // index in argument vector
00145         size_t user_m     = 0;       // size of result vector
00146         size_t user_n     = 0;       // size of arugment vector
00147         // next expected operator in a UserOp sequence
00148         enum { user_start, user_arg, user_ret, user_end } user_state = user_start;
00149 
00150         // check numvar argument
00151         CPPAD_ASSERT_UNKNOWN( Rec->num_rec_var() == numvar );
00152 
00153         // length of the parameter vector (used by CppAD assert macros)
00154         const size_t num_par = Rec->num_rec_par();
00155 
00156         // length of the text vector (used by CppAD assert macros)
00157         const size_t num_text = Rec->num_rec_text();
00158 
00159         // pointer to the beginning of the parameter vector
00160         const Base* parameter = CPPAD_NULL;
00161         if( num_par > 0 )
00162                 parameter = Rec->GetPar();
00163 
00164         // pointer to the beginning of the text vector
00165         const char* text = CPPAD_NULL;
00166         if( num_text > 0 )
00167                 text = Rec->GetTxt(0);
00168 
00169         // skip the BeginOp at the beginning of the recording
00170         Rec->start_forward(op, arg, i_op, i_var);
00171         CPPAD_ASSERT_UNKNOWN( op == BeginOp );
00172 # if CPPAD_FORWARD0SWEEP_TRACE
00173         std::cout << std::endl;
00174 # endif
00175         while(op != EndOp)
00176         {
00177                 // this op
00178                 Rec->next_forward(op, arg, i_op, i_var);
00179 # ifndef NDEBUG
00180                 if( i_op <= n )
00181                 {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
00182                 }
00183                 else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
00184 # endif
00185 
00186                 // action to take depends on the case
00187                 switch( op )
00188                 {
00189                         case AbsOp:
00190                         forward_abs_op_0(i_var, arg[0], J, Taylor);
00191                         break;
00192                         // -------------------------------------------------
00193 
00194                         case AddvvOp:
00195                         forward_addvv_op_0(i_var, arg, parameter, J, Taylor);
00196                         break;
00197                         // -------------------------------------------------
00198 
00199                         case AddpvOp:
00200                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00201                         forward_addpv_op_0(i_var, arg, parameter, J, Taylor);
00202                         break;
00203                         // -------------------------------------------------
00204 
00205                         case AcosOp:
00206                         // results: sqrt(1 - x * x), acos(x) 
00207                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
00208                         forward_acos_op_0(i_var, arg[0], J, Taylor);
00209                         break;
00210                         // -------------------------------------------------
00211 
00212                         case AsinOp:
00213                         // results: sqrt(1 - x * x), asin(x) 
00214                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
00215                         forward_asin_op_0(i_var, arg[0], J, Taylor);
00216                         break;
00217                         // -------------------------------------------------
00218 
00219                         case AtanOp:
00220                         // results: 1 + x * x, atan(x) 
00221                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
00222                         forward_atan_op_0(i_var, arg[0], J, Taylor);
00223                         break;
00224                         // -------------------------------------------------
00225 
00226                         case CSumOp:
00227                         // CSumOp has a variable number of arguments and
00228                         // next_forward thinks it one has one argument.
00229                         // we must inform next_forward of this special case.
00230                         Rec->forward_csum(op, arg, i_op, i_var);
00231                         forward_csum_op(
00232                                 0, i_var, arg, num_par, parameter, J, Taylor
00233                         );
00234                         break;
00235 
00236                         // -------------------------------------------------
00237                         case CExpOp:
00238                         // Use the general case with d == 0 
00239                         // (could create an optimzied verison for this case)
00240                         forward_cond_op_0(
00241                                 i_var, arg, num_par, parameter, J, Taylor
00242                         );
00243                         break;
00244                         // ---------------------------------------------------
00245                         case ComOp:
00246                         forward_comp_op_0(
00247                         compareCount, arg, num_par, parameter, J, Taylor
00248                         );
00249                         break;
00250                         // ---------------------------------------------------
00251 
00252                         case CosOp:
00253                         // results: sin(x), cos(x)
00254                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
00255                         forward_cos_op_0(i_var, arg[0], J, Taylor);
00256                         break;
00257                         // ---------------------------------------------------
00258 
00259                         case CoshOp:
00260                         // variables: sinh(x), cosh(x)
00261                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
00262                         forward_cosh_op_0(i_var, arg[0], J, Taylor);
00263                         break;
00264                         // -------------------------------------------------
00265 
00266                         case DisOp:
00267                         forward_dis_op_0(i_var, arg, J, Taylor);
00268                         break;
00269                         // -------------------------------------------------
00270 
00271                         case DivvvOp:
00272                         forward_divvv_op_0(i_var, arg, parameter, J, Taylor);
00273                         break;
00274                         // -------------------------------------------------
00275 
00276                         case DivpvOp:
00277                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00278                         forward_divpv_op_0(i_var, arg, parameter, J, Taylor);
00279                         break;
00280                         // -------------------------------------------------
00281 
00282                         case DivvpOp:
00283                         CPPAD_ASSERT_UNKNOWN( arg[1] < num_par );
00284                         forward_divvp_op_0(i_var, arg, parameter, J, Taylor);
00285                         break;
00286                         // -------------------------------------------------
00287 
00288                         case EndOp:
00289                         CPPAD_ASSERT_NARG_NRES(op, 0, 0);
00290                         break;
00291                         // -------------------------------------------------
00292 
00293                         case ExpOp:
00294                         forward_exp_op_0(i_var, arg[0], J, Taylor);
00295                         break;
00296                         // -------------------------------------------------
00297 
00298                         case InvOp:
00299                         break;
00300                         // -------------------------------------------------
00301 
00302                         case LdpOp:
00303                         CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00304                         CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00305                         non_const_arg = Rec->forward_non_const_arg();
00306                         forward_load_p_op_0(
00307                                 i_var, 
00308                                 non_const_arg, 
00309                                 num_par, 
00310                                 parameter, 
00311                                 J, 
00312                                 Taylor,
00313                                 Rec->num_rec_vecad_ind(),
00314                                 VectorVar,
00315                                 VectorInd
00316                         );
00317                         break;
00318                         // -------------------------------------------------
00319 
00320                         case LdvOp:
00321                         CPPAD_ASSERT_UNKNOWN( VectorInd != CPPAD_NULL );
00322                         CPPAD_ASSERT_UNKNOWN( VectorVar != CPPAD_NULL );
00323                         non_const_arg = Rec->forward_non_const_arg();
00324                         forward_load_v_op_0(
00325                                 i_var, 
00326                                 non_const_arg, 
00327                                 num_par, 
00328                                 parameter, 
00329                                 J, 
00330                                 Taylor,
00331                                 Rec->num_rec_vecad_ind(),
00332                                 VectorVar,
00333                                 VectorInd
00334                         );
00335                         break;
00336                         // -------------------------------------------------
00337 
00338                         case LogOp:
00339                         forward_log_op_0(i_var, arg[0], J, Taylor);
00340                         break;
00341                         // -------------------------------------------------
00342 
00343                         case MulvvOp:
00344                         forward_mulvv_op_0(i_var, arg, parameter, J, Taylor);
00345                         break;
00346                         // -------------------------------------------------
00347 
00348                         case MulpvOp:
00349                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00350                         forward_mulpv_op_0(i_var, arg, parameter, J, Taylor);
00351                         break;
00352                         // -------------------------------------------------
00353 
00354                         case ParOp:
00355                         forward_par_op_0(
00356                                 i_var, arg, num_par, parameter, J, Taylor
00357                         );
00358                         break;
00359                         // -------------------------------------------------
00360 
00361                         case PowvpOp:
00362                         CPPAD_ASSERT_UNKNOWN( arg[1] < num_par );
00363                         forward_powvp_op_0(i_var, arg, parameter, J, Taylor);
00364                         break;
00365                         // -------------------------------------------------
00366 
00367                         case PowpvOp:
00368                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00369                         forward_powpv_op_0(i_var, arg, parameter, J, Taylor);
00370                         break;
00371                         // -------------------------------------------------
00372 
00373                         case PowvvOp:
00374                         forward_powvv_op_0(i_var, arg, parameter, J, Taylor);
00375                         break;
00376                         // -------------------------------------------------
00377 
00378                         case PripOp:
00379                         if( print ) forward_prip_0(
00380                                 arg, num_text, text, num_par, parameter
00381                         );
00382                         break;
00383                         // -------------------------------------------------
00384 
00385                         case PrivOp:
00386                         if( print ) forward_priv_0(
00387                                 i_var, arg, num_text, text, J, Taylor
00388                         );
00389                         break;
00390                         // -------------------------------------------------
00391 
00392                         case SinOp:
00393                         // variables: cos(x), sin(x)
00394                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
00395                         forward_sin_op_0(i_var, arg[0], J, Taylor);
00396                         break;
00397                         // -------------------------------------------------
00398 
00399                         case SinhOp:
00400                         // variables: cosh(x), sinh(x)
00401                         CPPAD_ASSERT_UNKNOWN( i_var < numvar  );
00402                         forward_sinh_op_0(i_var, arg[0], J, Taylor);
00403                         break;
00404                         // -------------------------------------------------
00405 
00406                         case SqrtOp:
00407                         forward_sqrt_op_0(i_var, arg[0], J, Taylor);
00408                         break;
00409                         // -------------------------------------------------
00410 
00411                         case StppOp:
00412                         forward_store_pp_op_0(
00413                                 i_var, 
00414                                 arg, 
00415                                 num_par, 
00416                                 J, 
00417                                 Taylor,
00418                                 Rec->num_rec_vecad_ind(),
00419                                 VectorVar,
00420                                 VectorInd
00421                         );
00422                         break;
00423                         // -------------------------------------------------
00424 
00425                         case StpvOp:
00426                         forward_store_pv_op_0(
00427                                 i_var, 
00428                                 arg, 
00429                                 num_par, 
00430                                 J, 
00431                                 Taylor,
00432                                 Rec->num_rec_vecad_ind(),
00433                                 VectorVar,
00434                                 VectorInd
00435                         );
00436                         break;
00437                         // -------------------------------------------------
00438 
00439                         case StvpOp:
00440                         forward_store_vp_op_0(
00441                                 i_var, 
00442                                 arg, 
00443                                 num_par, 
00444                                 J, 
00445                                 Taylor,
00446                                 Rec->num_rec_vecad_ind(),
00447                                 VectorVar,
00448                                 VectorInd
00449                         );
00450                         break;
00451                         // -------------------------------------------------
00452 
00453                         case StvvOp:
00454                         forward_store_vv_op_0(
00455                                 i_var, 
00456                                 arg, 
00457                                 num_par, 
00458                                 J, 
00459                                 Taylor,
00460                                 Rec->num_rec_vecad_ind(),
00461                                 VectorVar,
00462                                 VectorInd
00463                         );
00464                         break;
00465                         // -------------------------------------------------
00466 
00467                         case SubvvOp:
00468                         forward_subvv_op_0(i_var, arg, parameter, J, Taylor);
00469                         break;
00470                         // -------------------------------------------------
00471 
00472                         case SubpvOp:
00473                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00474                         forward_subpv_op_0(i_var, arg, parameter, J, Taylor);
00475                         break;
00476                         // -------------------------------------------------
00477 
00478                         case SubvpOp:
00479                         CPPAD_ASSERT_UNKNOWN( arg[1] < num_par );
00480                         forward_subvp_op_0(i_var, arg, parameter, J, Taylor);
00481                         break;
00482                         // -------------------------------------------------
00483 
00484                         case UserOp:
00485                         // start or end an atomic operation sequence
00486                         CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
00487                         CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
00488                         if( user_state == user_start )
00489                         {       user_index = arg[0];
00490                                 user_id    = arg[1];
00491                                 user_n     = arg[2];
00492                                 user_m     = arg[3];
00493                                 if(user_tx.size() < user_n)
00494                                         user_tx.resize(user_n);
00495                                 if(user_ty.size() < user_m)
00496                                         user_ty.resize(user_m);
00497                                 user_j     = 0;
00498                                 user_i     = 0;
00499                                 user_state = user_arg;
00500                         }
00501                         else
00502                         {       CPPAD_ASSERT_UNKNOWN( user_state == user_end );
00503                                 CPPAD_ASSERT_UNKNOWN( user_index == arg[0] );
00504                                 CPPAD_ASSERT_UNKNOWN( user_id    == arg[1] );
00505                                 CPPAD_ASSERT_UNKNOWN( user_n     == arg[2] );
00506                                 CPPAD_ASSERT_UNKNOWN( user_m     == arg[3] );
00507                                 user_state = user_start;
00508                         }
00509                         break;
00510 
00511                         case UsrapOp:
00512                         // parameter argument in an atomic operation sequence
00513                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00514                         CPPAD_ASSERT_UNKNOWN( user_j < user_n );
00515                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00516                         user_tx[user_j++] = parameter[ arg[0] ];
00517                         if( user_j == user_n )
00518                         {       // call users function for this operation
00519                                 user_atomic<Base>::forward(user_index, user_id, 
00520                                         user_k, user_n, user_m, user_tx, user_ty
00521                                 );
00522                                 user_state = user_ret;
00523                         }
00524                         break;
00525 
00526                         case UsravOp:
00527                         // variable argument in an atomic operation sequence
00528                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00529                         CPPAD_ASSERT_UNKNOWN( user_j < user_n );
00530                         CPPAD_ASSERT_UNKNOWN( arg[0] <= i_var );
00531                         user_tx[user_j++] = Taylor[ arg[0] * J + 0 ];
00532                         if( user_j == user_n )
00533                         {       // call users function for this operation
00534                                 user_atomic<Base>::forward(user_index, user_id,
00535                                         user_k, user_n, user_m, user_tx, user_ty
00536                                 );
00537                                 user_state = user_ret;
00538                         }
00539                         break;
00540 
00541                         case UsrrpOp:
00542                         // parameter result in an atomic operation sequence
00543                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00544                         CPPAD_ASSERT_UNKNOWN( user_i < user_m );
00545                         user_i++;
00546                         if( user_i == user_m )
00547                                 user_state = user_end;
00548                         break;
00549 
00550                         case UsrrvOp:
00551                         // variable result in an atomic operation sequence
00552                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00553                         CPPAD_ASSERT_UNKNOWN( user_i < user_m );
00554                         Taylor[ i_var * J + 0 ] = user_ty[user_i++];
00555                         if( user_i == user_m )
00556                                 user_state = user_end;
00557                         break;
00558                         // -------------------------------------------------
00559 
00560                         default:
00561                         CPPAD_ASSERT_UNKNOWN(0);
00562                 }
00563 # if CPPAD_FORWARD0SWEEP_TRACE
00564                 size_t       d      = 0;
00565                 size_t       i_tmp  = i_var;
00566                 Base*        Z_tmp  = Taylor + i_var * J;
00567 
00568                 printOp(
00569                         std::cout, 
00570                         Rec,
00571                         i_tmp,
00572                         op, 
00573                         arg,
00574                         d + 1, 
00575                         Z_tmp, 
00576                         0, 
00577                         (Base *) CPPAD_NULL
00578                 );
00579         }
00580         std::cout << std::endl;
00581 # else
00582         }
00583 # endif
00584         CPPAD_ASSERT_UNKNOWN( user_state == user_start );
00585         CPPAD_ASSERT_UNKNOWN( i_var + 1 == Rec->num_rec_var() );
00586 
00587         if( VectorInd != CPPAD_NULL )
00588                 CPPAD_TRACK_DEL_VEC(VectorInd);
00589         if( VectorVar != CPPAD_NULL )
00590                 CPPAD_TRACK_DEL_VEC(VectorVar);
00591 
00592         return compareCount;
00593 }
00594 
00595 # undef CPPAD_FORWARD0SWEEP_TRACE
00596 
00597 CPPAD_END_NAMESPACE
00598 # endif