diff --git a/appendix/whats_new/2024.xrst b/appendix/whats_new/2024.xrst index 22022687f..b69ae1af0 100644 --- a/appendix/whats_new/2024.xrst +++ b/appendix/whats_new/2024.xrst @@ -34,6 +34,13 @@ Release Notes for 2024 mm-dd ***** +11-16 +===== +The :ref:`for_hes_sparsity-name` and :ref:`ForSparseHes-name` routines +did not handle the cumulative summation operators properly. +These operators sum more than two variables or dynamic parameters at a time +and may be created during ADFun :ref:`optimization ` . + 11-13 ===== The :ref:`for_hes_sparsity-name` and :ref:`ForSparseHes-name` routines diff --git a/include/cppad/local/sweep/for_hes.hpp b/include/cppad/local/sweep/for_hes.hpp index 7ab20bc7e..5e0bc1f1a 100644 --- a/include/cppad/local/sweep/for_hes.hpp +++ b/include/cppad/local/sweep/for_hes.hpp @@ -8,6 +8,7 @@ # include # include # include +# include /* {xrst_begin local_sweep_for_hes dev} @@ -182,19 +183,17 @@ void for_hes( // CPPAD_ASSERT_UNKNOWN( numvar > 0 ); // - // vecad_sparsity contains a sparsity pattern for each VecAD object. - // vecad_ind maps a VecAD index (beginning of the VecAD object) - // to the index for the corresponding set in vecad_sparsity. + // vecad_sparsity: forward Jacobian sparsity pattern for each VecAD object. + // vecad_ind: maps the VecAD index at beginning of the VecAD object + // to the index for the corresponding set in vecad_sparsity. size_t num_vecad_ind = play->num_var_vecad_ind_rec(); size_t num_vecad_vec = play->num_var_vecad_rec(); SetVector vecad_sparsity; pod_vector vecad_ind; - pod_vector vecad_jac; if( num_vecad_vec > 0 ) { size_t length; vecad_sparsity.resize(num_vecad_vec, np1); vecad_ind.extend(num_vecad_ind); - vecad_jac.extend(num_vecad_vec); size_t j = 0; for(size_t i = 0; i < num_vecad_vec; i++) { // length of this VecAD @@ -206,8 +205,6 @@ void for_hes( vecad_ind[j+k] = num_vecad_vec; // start of next VecAD j += length + 1; - // initialize this vector's reverse jacobian value - vecad_jac[i] = false; } CPPAD_ASSERT_UNKNOWN( j == play->num_var_vecad_ind_rec() ); } @@ -409,6 +406,7 @@ void for_hes( // ------------------------------------------------- case CSumOp: + var_op::csum_forward_hes(arg, i_var, n, for_hes_sparse); itr.correct_before_increment(); break; // ------------------------------------------------- diff --git a/include/cppad/local/var_op/csum_op.hpp b/include/cppad/local/var_op/csum_op.hpp index d0512a6db..e300f761a 100644 --- a/include/cppad/local/var_op/csum_op.hpp +++ b/include/cppad/local/var_op/csum_op.hpp @@ -705,18 +705,43 @@ Vector_set is the type used for vectors of sets. It must satisfy the :ref:`SetVector-name` concept. +n +* +is the number of independent variables on the tape. + +for_hes_sparse +************** +see :ref:`local_sweep_for_hes@for_hes_sparse` . + {xrst_end var_csum_forward_hes} */ // BEGIN_CSUM_FORWARD_HES template -inline void load_forward_hes( +inline void csum_forward_hes( const addr_t* arg , size_t i_z , - const Vector_set& var_sparsity , - const bool* var_rev_jac ) + size_t n , + Vector_set& for_hes_sparse ) // END_CSUM_FORWARD_HES -{ - +{ // + // np1 + size_t np1 = n + 1; + // + // for_hes_sparse + for_hes_sparse.clear(np1 + i_z); + // + // addition and subtraction variables + for(addr_t i = 5; i < arg[2]; ++i) + { CPPAD_ASSERT_UNKNOWN( size_t(arg[i]) < i_z ); + // + // linear functions only modify forward Jacobian sparsity + for_hes_sparse.binary_union( + np1 + i_z , // index in sparsity for result + np1 + i_z , // index in sparsity for left operand + np1 + size_t(arg[i]) , // index for right operand + for_hes_sparse // sparsity vector for right operand + ); + } } } } } // END namespace diff --git a/test_more/general/for_hes_sparsity.cpp b/test_more/general/for_hes_sparsity.cpp index 5ce59b77a..6c993495c 100644 --- a/test_more/general/for_hes_sparsity.cpp +++ b/test_more/general/for_hes_sparsity.cpp @@ -198,9 +198,7 @@ bool for_hes_sparsity(void) { bool ok = true; // ok &= test_old_interface(); - // - // 2DO: fix this case - ok &= ! test_csum(); + ok &= test_csum(); // return ok; }