CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_REV_HES_SWEEP_INCLUDED 00003 # define CPPAD_REV_HES_SWEEP_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00017 /*! 00018 \file rev_hes_sweep.hpp 00019 Compute Reverse mode Hessian sparsity patterns. 00020 */ 00021 00022 /*! 00023 \def CPPAD_REV_HES_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 rev_hes_sweep computation is printed. 00027 */ 00028 # define CPPAD_REV_HES_SWEEP_TRACE 0 00029 00030 /*! 00031 Given the forward Jacobian sparsity pattern for all the variables, 00032 and the reverse Jacobian sparsity pattern for the dependent variables, 00033 RevHesSweep computes the Hessian sparsity pattern for all the independent 00034 variables. 00035 00036 \tparam Base 00037 base type for the operator; i.e., this operation sequence was recorded 00038 using AD< \a Base > and computations by this routine are done using type 00039 \a Base. 00040 00041 \tparam Vector_set 00042 is the type used for vectors of sets. It can be either 00043 \c sparse_pack, \c sparse_set, or \c sparse_list. 00044 00045 \param n 00046 is the number of independent variables on the tape. 00047 00048 \param numvar 00049 is the total number of variables on the tape; i.e., 00050 \a play->num_var_rec(). 00051 This is also the number of rows in the entire sparsity pattern 00052 \a rev_hes_sparse. 00053 00054 \param play 00055 The information stored in \a play 00056 is a recording of the operations corresponding to a 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 00061 and \f$ m \f$ is the number of dependent variables. 00062 The object \a play is effectly constant. 00063 It is not declared const because while playing back the tape 00064 the object \a play holds information about the currentl location 00065 with in the tape and this changes during playback. 00066 00067 \param for_jac_sparse 00068 For i = 0 , ... , \a numvar - 1, 00069 (for all the variables on the tape), 00070 the forward Jacobian sparsity pattern for the variable with index i 00071 corresponds to the set with index i in \a for_jac_sparse. 00072 00073 \param RevJac 00074 \b Input: 00075 For i = 0, ... , \a numvar - 1 00076 the if the variable with index i on the tape is an dependent variable and 00077 included in the Hessian, \a RevJac[ i ] is equal to true, 00078 otherwise it is equal to false. 00079 \n 00080 \n 00081 \b Output: The values in \a RevJac upon return are not specified; i.e., 00082 it is used for temporary work space. 00083 00084 \param rev_hes_sparse 00085 The reverse Hessian sparsity pattern for the variable with index i 00086 corresponds to the set with index i in \a rev_hes_sparse. 00087 \n 00088 \n 00089 \b Input: For i = 0 , ... , \a numvar - 1 00090 the reverse Hessian sparsity pattern for the variable with index i is empty. 00091 \n 00092 \n 00093 \b Output: For j = 1 , ... , \a n, 00094 the reverse Hessian sparsity pattern for the independent dependent variable 00095 with index (j-1) is given by the set with index j 00096 in \a rev_hes_sparse. 00097 The values in the rest of \a rev_hes_sparse are not specified; i.e., 00098 they are used for temporary work space. 00099 */ 00100 00101 template <class Base, class Vector_set> 00102 void RevHesSweep( 00103 size_t n, 00104 size_t numvar, 00105 player<Base> *play, 00106 Vector_set& for_jac_sparse, // should be const 00107 bool* RevJac, 00108 Vector_set& rev_hes_sparse 00109 ) 00110 { 00111 OpCode op; 00112 size_t i_op; 00113 size_t i_var; 00114 00115 const addr_t* arg = CPPAD_NULL; 00116 00117 // length of the parameter vector (used by CppAD assert macros) 00118 const size_t num_par = play->num_par_rec(); 00119 00120 size_t i, j, k; 00121 00122 // check numvar argument 00123 CPPAD_ASSERT_UNKNOWN( play->num_var_rec() == numvar ); 00124 CPPAD_ASSERT_UNKNOWN( for_jac_sparse.n_set() == numvar ); 00125 CPPAD_ASSERT_UNKNOWN( rev_hes_sparse.n_set() == numvar ); 00126 CPPAD_ASSERT_UNKNOWN( numvar > 0 ); 00127 00128 // upper limit exclusive for set elements 00129 size_t limit = rev_hes_sparse.end(); 00130 CPPAD_ASSERT_UNKNOWN( for_jac_sparse.end() == limit ); 00131 00132 // check number of sets match 00133 CPPAD_ASSERT_UNKNOWN( 00134 for_jac_sparse.n_set() == rev_hes_sparse.n_set() 00135 ); 00136 00137 // vecad_sparsity contains a sparsity pattern for each VecAD object. 00138 // vecad_ind maps a VecAD index (beginning of the VecAD object) 00139 // to the index for the corresponding set in vecad_sparsity. 00140 size_t num_vecad_ind = play->num_vec_ind_rec(); 00141 size_t num_vecad_vec = play->num_vecad_vec_rec(); 00142 Vector_set vecad_sparse; 00143 vecad_sparse.resize(num_vecad_vec, limit); 00144 pod_vector<size_t> vecad_ind; 00145 pod_vector<bool> vecad_jac; 00146 if( num_vecad_vec > 0 ) 00147 { size_t length; 00148 vecad_ind.extend(num_vecad_ind); 00149 vecad_jac.extend(num_vecad_vec); 00150 j = 0; 00151 for(i = 0; i < num_vecad_vec; i++) 00152 { // length of this VecAD 00153 length = play->GetVecInd(j); 00154 // set vecad_ind to proper index for this VecAD 00155 vecad_ind[j] = i; 00156 // make all other values for this vector invalid 00157 for(k = 1; k <= length; k++) 00158 vecad_ind[j+k] = num_vecad_vec; 00159 // start of next VecAD 00160 j += length + 1; 00161 // initialize this vector's reverse jacobian value 00162 vecad_jac[i] = false; 00163 } 00164 CPPAD_ASSERT_UNKNOWN( j == play->num_vec_ind_rec() ); 00165 } 00166 00167 // work space used by UserOp. 00168 vector<size_t> user_ix; // variable indices for argument vector x 00169 typedef std::set<size_t> size_set; 00170 size_set::iterator set_itr; // iterator for a standard set 00171 size_set::iterator set_end; // end of iterator sequence 00172 vector< size_set > set_r; // forward Jacobian sparsity for x 00173 vector< size_set > set_u; // reverse Hessian sparsity for y 00174 vector< size_set > set_v; // reverse Hessian sparsity for x 00175 // 00176 vector<bool> bool_r; // bool forward Jacobian sparsity for x 00177 vector<bool> bool_u; // bool reverse Hessian sparsity for y 00178 vector<bool> bool_v; // bool reverse Hessian sparsity for x 00179 // 00180 vector<bool> user_vx; // which components of x are variables 00181 vector<bool> user_s; // reverse Jacobian sparsity for y 00182 vector<bool> user_t; // reverse Jacobian sparsity for x 00183 const size_t user_q = limit; // maximum element plus one 00184 size_t user_index = 0; // indentifier for this atomic operation 00185 size_t user_id = 0; // user identifier for this call to operator 00186 size_t user_i = 0; // index in result vector 00187 size_t user_j = 0; // index in argument vector 00188 size_t user_m = 0; // size of result vector 00189 size_t user_n = 0; // size of arugment vector 00190 // 00191 atomic_base<Base>* user_atom = CPPAD_NULL; // user's atomic op calculator 00192 bool user_bool = false; // use bool or set sparsity ? 00193 # ifndef NDEBUG 00194 bool user_ok = false; // atomic op return value 00195 # endif 00196 // next expected operator in a UserOp sequence 00197 enum { user_start, user_arg, user_ret, user_end } user_state = user_end; 00198 00199 00200 // Initialize 00201 play->reverse_start(op, arg, i_op, i_var); 00202 CPPAD_ASSERT_UNKNOWN( op == EndOp ); 00203 # if CPPAD_REV_HES_SWEEP_TRACE 00204 std::cout << std::endl; 00205 CppAD::vectorBool zf_value(limit); 00206 CppAD::vectorBool zh_value(limit); 00207 # endif 00208 bool more_operators = true; 00209 while(more_operators) 00210 { 00211 // next op 00212 play->reverse_next(op, arg, i_op, i_var); 00213 # ifndef NDEBUG 00214 if( i_op <= n ) 00215 { CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp)); 00216 } 00217 else CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp)); 00218 # endif 00219 00220 // rest of information depends on the case 00221 switch( op ) 00222 { 00223 case AbsOp: 00224 CPPAD_ASSERT_NARG_NRES(op, 1, 1) 00225 reverse_sparse_hessian_linear_unary_op( 00226 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00227 ); 00228 break; 00229 // ------------------------------------------------- 00230 00231 case AddvvOp: 00232 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00233 reverse_sparse_hessian_addsub_op( 00234 i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse 00235 ); 00236 break; 00237 // ------------------------------------------------- 00238 00239 case AddpvOp: 00240 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00241 reverse_sparse_hessian_linear_unary_op( 00242 i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse 00243 ); 00244 break; 00245 // ------------------------------------------------- 00246 00247 case AcosOp: 00248 // sqrt(1 - x * x), acos(x) 00249 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00250 reverse_sparse_hessian_nonlinear_unary_op( 00251 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00252 ); 00253 break; 00254 // ------------------------------------------------- 00255 00256 case AsinOp: 00257 // sqrt(1 - x * x), asin(x) 00258 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00259 reverse_sparse_hessian_nonlinear_unary_op( 00260 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00261 ); 00262 break; 00263 // ------------------------------------------------- 00264 00265 case AtanOp: 00266 // 1 + x * x, atan(x) 00267 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00268 reverse_sparse_hessian_nonlinear_unary_op( 00269 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00270 ); 00271 break; 00272 // ------------------------------------------------- 00273 00274 case BeginOp: 00275 CPPAD_ASSERT_NARG_NRES(op, 1, 1) 00276 more_operators = false; 00277 break; 00278 // ------------------------------------------------- 00279 00280 case CSkipOp: 00281 // CSkipOp has a variable number of arguments and 00282 // reverse_next thinks it one has one argument. 00283 // We must inform reverse_next of this special case. 00284 play->reverse_cskip(op, arg, i_op, i_var); 00285 break; 00286 // ------------------------------------------------- 00287 00288 case CSumOp: 00289 // CSumOp has a variable number of arguments and 00290 // reverse_next thinks it one has one argument. 00291 // We must inform reverse_next of this special case. 00292 play->reverse_csum(op, arg, i_op, i_var); 00293 reverse_sparse_hessian_csum_op( 00294 i_var, arg, RevJac, rev_hes_sparse 00295 ); 00296 break; 00297 // ------------------------------------------------- 00298 00299 case CExpOp: 00300 reverse_sparse_hessian_cond_op( 00301 i_var, arg, num_par, RevJac, rev_hes_sparse 00302 ); 00303 break; 00304 // --------------------------------------------------- 00305 00306 case ComOp: 00307 CPPAD_ASSERT_NARG_NRES(op, 4, 0) 00308 CPPAD_ASSERT_UNKNOWN( arg[1] > 1 ); 00309 break; 00310 // -------------------------------------------------- 00311 00312 case CosOp: 00313 // sin(x), cos(x) 00314 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00315 reverse_sparse_hessian_nonlinear_unary_op( 00316 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00317 ); 00318 break; 00319 // --------------------------------------------------- 00320 00321 case CoshOp: 00322 // sinh(x), cosh(x) 00323 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00324 reverse_sparse_hessian_nonlinear_unary_op( 00325 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00326 ); 00327 break; 00328 // ------------------------------------------------- 00329 00330 case DisOp: 00331 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00332 // derivativve is identically zero 00333 break; 00334 // ------------------------------------------------- 00335 00336 case DivvvOp: 00337 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00338 reverse_sparse_hessian_div_op( 00339 i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse 00340 ); 00341 break; 00342 // ------------------------------------------------- 00343 00344 case DivpvOp: 00345 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00346 reverse_sparse_hessian_nonlinear_unary_op( 00347 i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse 00348 ); 00349 break; 00350 // ------------------------------------------------- 00351 00352 case DivvpOp: 00353 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00354 reverse_sparse_hessian_linear_unary_op( 00355 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00356 ); 00357 break; 00358 // ------------------------------------------------- 00359 00360 case ExpOp: 00361 CPPAD_ASSERT_NARG_NRES(op, 1, 1) 00362 reverse_sparse_hessian_nonlinear_unary_op( 00363 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00364 ); 00365 break; 00366 // ------------------------------------------------- 00367 00368 case InvOp: 00369 CPPAD_ASSERT_NARG_NRES(op, 0, 1) 00370 // Z is already defined 00371 break; 00372 // ------------------------------------------------- 00373 00374 case LdpOp: 00375 reverse_sparse_hessian_load_op( 00376 op, 00377 i_var, 00378 arg, 00379 num_vecad_ind, 00380 vecad_ind.data(), 00381 rev_hes_sparse, 00382 vecad_sparse, 00383 RevJac, 00384 vecad_jac.data() 00385 ); 00386 break; 00387 // ------------------------------------------------- 00388 00389 case LdvOp: 00390 reverse_sparse_hessian_load_op( 00391 op, 00392 i_var, 00393 arg, 00394 num_vecad_ind, 00395 vecad_ind.data(), 00396 rev_hes_sparse, 00397 vecad_sparse, 00398 RevJac, 00399 vecad_jac.data() 00400 ); 00401 break; 00402 // ------------------------------------------------- 00403 00404 case LogOp: 00405 CPPAD_ASSERT_NARG_NRES(op, 1, 1) 00406 reverse_sparse_hessian_nonlinear_unary_op( 00407 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00408 ); 00409 break; 00410 // ------------------------------------------------- 00411 00412 case MulvvOp: 00413 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00414 reverse_sparse_hessian_mul_op( 00415 i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse 00416 ); 00417 break; 00418 // ------------------------------------------------- 00419 00420 case MulpvOp: 00421 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00422 reverse_sparse_hessian_linear_unary_op( 00423 i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse 00424 ); 00425 break; 00426 // ------------------------------------------------- 00427 00428 case ParOp: 00429 CPPAD_ASSERT_NARG_NRES(op, 1, 1) 00430 00431 break; 00432 // ------------------------------------------------- 00433 00434 case PowpvOp: 00435 CPPAD_ASSERT_NARG_NRES(op, 2, 3) 00436 reverse_sparse_hessian_nonlinear_unary_op( 00437 i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse 00438 ); 00439 break; 00440 // ------------------------------------------------- 00441 00442 case PowvpOp: 00443 CPPAD_ASSERT_NARG_NRES(op, 2, 3) 00444 reverse_sparse_hessian_nonlinear_unary_op( 00445 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00446 ); 00447 break; 00448 // ------------------------------------------------- 00449 00450 case PowvvOp: 00451 CPPAD_ASSERT_NARG_NRES(op, 2, 3) 00452 reverse_sparse_hessian_pow_op( 00453 i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse 00454 ); 00455 break; 00456 // ------------------------------------------------- 00457 00458 case PriOp: 00459 CPPAD_ASSERT_NARG_NRES(op, 5, 0); 00460 break; 00461 // ------------------------------------------------- 00462 00463 case SignOp: 00464 CPPAD_ASSERT_NARG_NRES(op, 1, 1); 00465 // Derivative is identiaclly zero 00466 break; 00467 // ------------------------------------------------- 00468 00469 case SinOp: 00470 // cos(x), sin(x) 00471 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00472 reverse_sparse_hessian_nonlinear_unary_op( 00473 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00474 ); 00475 break; 00476 // ------------------------------------------------- 00477 00478 case SinhOp: 00479 // cosh(x), sinh(x) 00480 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00481 reverse_sparse_hessian_nonlinear_unary_op( 00482 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00483 ); 00484 break; 00485 // ------------------------------------------------- 00486 00487 case SqrtOp: 00488 CPPAD_ASSERT_NARG_NRES(op, 1, 1) 00489 reverse_sparse_hessian_nonlinear_unary_op( 00490 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00491 ); 00492 break; 00493 // ------------------------------------------------- 00494 00495 case StppOp: 00496 // sparsity cannot propagate through a parameter 00497 CPPAD_ASSERT_NARG_NRES(op, 3, 0) 00498 break; 00499 // ------------------------------------------------- 00500 00501 case StpvOp: 00502 reverse_sparse_hessian_store_op( 00503 op, 00504 arg, 00505 num_vecad_ind, 00506 vecad_ind.data(), 00507 rev_hes_sparse, 00508 vecad_sparse, 00509 RevJac, 00510 vecad_jac.data() 00511 ); 00512 break; 00513 // ------------------------------------------------- 00514 00515 case StvpOp: 00516 // sparsity cannot propagate through a parameter 00517 CPPAD_ASSERT_NARG_NRES(op, 3, 0) 00518 break; 00519 // ------------------------------------------------- 00520 00521 case StvvOp: 00522 reverse_sparse_hessian_store_op( 00523 op, 00524 arg, 00525 num_vecad_ind, 00526 vecad_ind.data(), 00527 rev_hes_sparse, 00528 vecad_sparse, 00529 RevJac, 00530 vecad_jac.data() 00531 ); 00532 break; 00533 // ------------------------------------------------- 00534 00535 case SubvvOp: 00536 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00537 reverse_sparse_hessian_addsub_op( 00538 i_var, arg, RevJac, for_jac_sparse, rev_hes_sparse 00539 ); 00540 break; 00541 // ------------------------------------------------- 00542 00543 case SubpvOp: 00544 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00545 reverse_sparse_hessian_linear_unary_op( 00546 i_var, arg[1], RevJac, for_jac_sparse, rev_hes_sparse 00547 ); 00548 break; 00549 // ------------------------------------------------- 00550 00551 case SubvpOp: 00552 CPPAD_ASSERT_NARG_NRES(op, 2, 1) 00553 reverse_sparse_hessian_linear_unary_op( 00554 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00555 ); 00556 break; 00557 // ------------------------------------------------- 00558 00559 case TanOp: 00560 // tan(x)^2, tan(x) 00561 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00562 reverse_sparse_hessian_nonlinear_unary_op( 00563 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00564 ); 00565 break; 00566 // ------------------------------------------------- 00567 00568 case TanhOp: 00569 // tanh(x)^2, tanh(x) 00570 CPPAD_ASSERT_NARG_NRES(op, 1, 2) 00571 reverse_sparse_hessian_nonlinear_unary_op( 00572 i_var, arg[0], RevJac, for_jac_sparse, rev_hes_sparse 00573 ); 00574 break; 00575 // ------------------------------------------------- 00576 00577 case UserOp: 00578 // start or end an atomic operation sequence 00579 CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 ); 00580 CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 ); 00581 if( user_state == user_end ) 00582 { user_index = arg[0]; 00583 user_id = arg[1]; 00584 user_n = arg[2]; 00585 user_m = arg[3]; 00586 user_atom = atomic_base<Base>::class_object(user_index); 00587 # ifndef NDEBUG 00588 if( user_atom == CPPAD_NULL ) 00589 { std::string msg = 00590 atomic_base<Base>::class_name(user_index) 00591 + ": atomic_base function has been deleted"; 00592 CPPAD_ASSERT_KNOWN(false, msg.c_str() ); 00593 } 00594 # endif 00595 user_bool = user_atom->sparsity() == 00596 atomic_base<Base>::bool_sparsity_enum; 00597 user_ix.resize(user_n); 00598 user_vx.resize(user_n); 00599 user_s.resize(user_m); 00600 user_t.resize(user_n); 00601 00602 // simpler to initialize all sparsity patterns as empty 00603 for(i = 0; i < user_m; i++) 00604 user_s[i] = false; 00605 for(i = 0; i < user_n; i++) 00606 user_t[i] = false; 00607 if( user_bool ) 00608 { bool_r.resize(user_n * user_q); 00609 bool_u.resize(user_m * user_q); 00610 bool_v.resize(user_n * user_q); 00611 // simpler to initialize all patterns as empty 00612 for(i = 0; i < user_m; i++) 00613 { 00614 for(j = 0; j < user_q; j++) 00615 bool_u[ i * user_q + j] = false; 00616 } 00617 for(i = 0; i < user_n; i++) 00618 { 00619 for(j = 0; j < user_q; j++) 00620 { bool_r[ i * user_q + j] = false; 00621 bool_v[ i * user_q + j] = false; 00622 } 00623 } 00624 } 00625 else 00626 { set_r.resize(user_n); 00627 set_u.resize(user_m); 00628 set_v.resize(user_n); 00629 for(i = 0; i < user_m; i++) 00630 set_u[i].clear(); 00631 for(i = 0; i < user_n; i++) 00632 { set_r[i].clear(); 00633 set_v[i].clear(); 00634 } 00635 } 00636 user_j = user_n; 00637 user_i = user_m; 00638 user_state = user_ret; 00639 } 00640 else 00641 { CPPAD_ASSERT_UNKNOWN( user_state == user_start ); 00642 CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) ); 00643 CPPAD_ASSERT_UNKNOWN( user_id == size_t(arg[1]) ); 00644 CPPAD_ASSERT_UNKNOWN( user_n == size_t(arg[2]) ); 00645 CPPAD_ASSERT_UNKNOWN( user_m == size_t(arg[3]) ); 00646 user_state = user_end; 00647 00648 // call users function for this operation 00649 user_atom->set_id(user_id); 00650 # ifdef NDEBUG 00651 if( user_bool ) 00652 user_atom->rev_sparse_hes(user_vx, 00653 user_s, user_t, user_q, bool_r, bool_u, bool_v 00654 ); 00655 else 00656 user_atom->rev_sparse_hes(user_vx, 00657 user_s, user_t, user_q, set_r, set_u, set_v 00658 ); 00659 # else 00660 if( user_bool ) 00661 user_ok = user_atom->rev_sparse_hes(user_vx, 00662 user_s, user_t, user_q, bool_r, bool_u, bool_v 00663 ); 00664 else 00665 user_ok = user_atom->rev_sparse_hes(user_vx, 00666 user_s, user_t, user_q, set_r, set_u, set_v 00667 ); 00668 if( ! user_ok ) 00669 { std::string msg = 00670 atomic_base<Base>::class_name(user_index) 00671 + ": atomic_base.rev_sparse_hes: returned false"; 00672 CPPAD_ASSERT_KNOWN(false, msg.c_str() ); 00673 } 00674 # endif 00675 for(i = 0; i < user_n; i++) if( user_ix[i] > 0 ) 00676 { 00677 size_t i_x = user_ix[i]; 00678 if( user_t[i] ) 00679 RevJac[i_x] = true; 00680 if( user_bool ) 00681 { 00682 for(j = 0; j < user_q; j++) 00683 if( bool_v[ i * user_q + j ] ) 00684 rev_hes_sparse.add_element(i_x, j); 00685 } 00686 else 00687 { 00688 set_itr = set_v[i].begin(); 00689 set_end = set_v[i].end(); 00690 while( set_itr != set_end ) 00691 rev_hes_sparse.add_element(i_x, *set_itr++); 00692 } 00693 } 00694 } 00695 break; 00696 00697 case UsrapOp: 00698 // parameter argument in an atomic operation sequence 00699 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00700 CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n ); 00701 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00702 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00703 --user_j; 00704 user_ix[user_j] = 0; 00705 user_vx[user_j] = false; 00706 if( user_j == 0 ) 00707 user_state = user_start; 00708 break; 00709 00710 case UsravOp: 00711 // variable argument in an atomic operation sequence 00712 CPPAD_ASSERT_UNKNOWN( user_state == user_arg ); 00713 CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n ); 00714 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00715 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var ); 00716 CPPAD_ASSERT_UNKNOWN( 0 < arg[0] ); 00717 --user_j; 00718 user_ix[user_j] = arg[0]; 00719 user_vx[user_j] = true; 00720 for_jac_sparse.begin(arg[0]); 00721 i = for_jac_sparse.next_element(); 00722 while( i < user_q ) 00723 { if( user_bool ) 00724 bool_r[ user_j * user_q + i ] = true; 00725 else 00726 set_r[user_j].insert(i); 00727 i = for_jac_sparse.next_element(); 00728 } 00729 if( user_j == 0 ) 00730 user_state = user_start; 00731 break; 00732 00733 case UsrrpOp: 00734 // parameter result in an atomic operation sequence 00735 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00736 CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m ); 00737 CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 ); 00738 CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par ); 00739 --user_i; 00740 if( user_i == 0 ) 00741 user_state = user_arg; 00742 break; 00743 00744 case UsrrvOp: 00745 // variable result in an atomic operation sequence 00746 CPPAD_ASSERT_UNKNOWN( user_state == user_ret ); 00747 CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m ); 00748 --user_i; 00749 if( RevJac[i_var] ) 00750 { 00751 user_s[user_i] = true; 00752 } 00753 rev_hes_sparse.begin(i_var); 00754 j = rev_hes_sparse.next_element(); 00755 while( j < user_q ) 00756 { if( user_bool ) 00757 bool_u[user_i * user_q + j] = true; 00758 else 00759 set_u[user_i].insert(j); 00760 j = rev_hes_sparse.next_element(); 00761 } 00762 if( user_i == 0 ) 00763 user_state = user_arg; 00764 break; 00765 00766 // ------------------------------------------------- 00767 00768 default: 00769 CPPAD_ASSERT_UNKNOWN(0); 00770 } 00771 # if CPPAD_REV_HES_SWEEP_TRACE 00772 for(j = 0; j < limit; j++) 00773 { zf_value[j] = false; 00774 zh_value[j] = false; 00775 } 00776 for_jac_sparse.begin(i_var);; 00777 j = for_jac_sparse.next_element();; 00778 while( j < limit ) 00779 { zf_value[j] = true; 00780 j = for_jac_sparse.next_element(); 00781 } 00782 rev_hes_sparse.begin(i_var);; 00783 j = rev_hes_sparse.next_element();; 00784 while( j < limit ) 00785 { zh_value[j] = true; 00786 j = rev_hes_sparse.next_element(); 00787 } 00788 // should also print RevJac[i_var], but printOp does not 00789 // yet allow for this. 00790 printOp( 00791 std::cout, 00792 play, 00793 i_op, 00794 i_var, 00795 op, 00796 arg, 00797 1, 00798 &zf_value, 00799 1, 00800 &zh_value 00801 ); 00802 # endif 00803 } 00804 // values corresponding to BeginOp 00805 CPPAD_ASSERT_UNKNOWN( i_op == 0 ); 00806 CPPAD_ASSERT_UNKNOWN( i_var == 0 ); 00807 00808 return; 00809 } 00810 } // END_CPPAD_NAMESPACE 00811 00812 // preprocessor symbols that are local to this file 00813 # undef CPPAD_REV_HES_SWEEP_TRACE 00814 00815 # endif