CppAD: A C++ Algorithmic Differentiation Package 20110419
|
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