master: for_hes.hpp: fix bug in CSumOp operator.

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 <optimize-name>` .
+
 11-13
 =====
 The :ref:`for_hes_sparsity-name` and :ref:`ForSparseHes-name` routines

include/cppad/local/sweep/for_hes.hpp

@@ -8,6 +8,7 @@
 # include <cppad/local/play/atom_op_info.hpp>
 # include <cppad/local/var_op/load_op.hpp>
 # include <cppad/local/var_op/store_op.hpp>
+# include <cppad/local/var_op/csum_op.hpp>
 
 /*
 {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<size_t> vecad_ind;
-   pod_vector<bool>   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;
          // -------------------------------------------------

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 <class Vector_set>
-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

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;
 }