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