bool test_control_reduced_objective(void)
{    bool ok = true;
     double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
     typedef CppAD::AD<double> adouble;
     // ----------------------------------------------------------------------
     // control problem parameters
     size_t   J = 3;
     double   delta_t = 0.1;
     VECTOR(double) p(4), q(4);
     for(size_t i = 0; i < 4; i++)
     {    p[i] = double(5 - i);
          q[i] = double(i + 1);
     }
     // ----------------------------------------------------------------------
     // L_fun, F_fun
     CppAD::ADFun<double> F_fun, L_fun;
     control::rec_objective(F_fun, J,  delta_t, q);
     control::rec_constraint(L_fun, J, delta_t, p);
     size_t m = L_fun.Range();
     size_t n = L_fun.Domain() - m;
     ok &= F_fun.Range() == 1;
     ok &= F_fun.Domain() == n + m;
     // -----------------------------------------------------------------------
     // aL_fun
     CppAD::ADFun<adouble> aL_fun;
     adouble adelta_t = delta_t;
     VECTOR(adouble) ap(4);
     for(size_t i = 0; i < 4; i++)
          ap[i] = p[i];
     control::rec_constraint(aL_fun, J, adelta_t, ap);
     // -----------------------------------------------------------------------
     // control_solve
     double criteria = double(n + m) * eps99;
     control::implicit_solver control_solve(L_fun, aL_fun, criteria);
     // -----------------------------------------------------------------------
     // R_kedem, R_full, R_partial
     implicit_kedem<control::implicit_solver> R_kedem(
          L_fun, F_fun, control_solve
     );
     bool full_step  = true;
     size_t num_step = 2;
     implicit_newton<control::implicit_solver> R_full(
          full_step, num_step, aL_fun, F_fun, control_solve
     );
     full_step  = false;
     implicit_newton<control::implicit_solver> R_partial(
          full_step, num_step, aL_fun, F_fun, control_solve
     );
     // ----------------------------------------------------------------------
     // zero order forward
     VECTOR(double) x0(n);
     for(size_t i = 0; i < n; i++)
          x0[i] = 0.1;
     bool store = true;
     VECTOR(double) kedem_R0     = R_kedem.Forward(0, x0, store);
     VECTOR(double) full_R0      = R_full.Forward(0, x0);
     VECTOR(double) partial_R0   = R_partial.Forward(0, x0);
     ok &= CppAD::NearEqual(kedem_R0[0], partial_R0[0], eps99, eps99);
     ok &= CppAD::NearEqual(kedem_R0[0], full_R0[0],    eps99, eps99);
     // -----------------------------------------------------------------------
     // compute and check gradient w.r.t. x
     VECTOR(double) w(1);
     w[0] = 1.0;
     VECTOR(double) partial_dR = R_partial.Reverse(1, w);
     VECTOR(double) full_dR    = R_partial.Reverse(1, w);
     VECTOR(double) x1(n), kedem_dR(1);
     for(size_t i = 0; i < n; i++)
          x1[i] = 0.0;
     for(size_t j = 0; j < n; j++)
     {    x1[j] = 1.0;
          store = false;
          kedem_dR  = R_kedem.Forward(1, x1, store);
          ok &= CppAD::NearEqual(kedem_dR[0], partial_dR[j], eps99, eps99);
          ok &= CppAD::NearEqual(kedem_dR[0], full_dR[j],    eps99, eps99);
          x1[j] = 0.0;
     }
     // -----------------------------------------------------------------------
     // compute Hessian w.r.t. x corresponding to R_partial or R_full
     MATRIX(double) H(n, n);
     VECTOR(double) d_partial(n), d_full(n);
     w.resize(2);
     w[0] = 0.0;
     w[1] = 1.0;
     for(size_t i = 0; i < n; i++)
          x1[i] = 0.0;
     for(size_t i = 0; i < n; i++)
     {    x1[i] = 1.0;
          //
          // first order forward
          R_partial.Forward(1, x1);
          R_full.Forward(1, x1);
          //
          // second order reverse
          d_partial = R_partial.Reverse(2, w);
          d_full    = R_full.Reverse(2, w);
          for(size_t j = 0; j < n; j++)
          {    ok &= CppAD::NearEqual(
                    d_partial[2 * j + 0], d_full[2 * j + 0], eps99, eps99
               );
               ok &= CppAD::NearEqual(
                    d_partial[2 * j + 1], d_full[2 * j + 1], eps99, eps99
               );
               H(i, j) = d_partial[2 * j + 0];
          }
          x1[i] = 0.0;
     }
     // -----------------------------------------------------------------------
     // compare H with Hessian corresponding to kedem
     VECTOR(double) x2(n), ddR(1);
     for(size_t i = 0; i < n; i++)
          x2[i] = 0.0;
     for(size_t i = 0; i < n; i++)
     {    // check i-th diagonal element
          x1[i] = 1.0;
          //
          store = true;
          R_kedem.Forward(1, x1, true);
          store = false;
          ddR = R_kedem.Forward(2, x2, false);
          ok &= CppAD::NearEqual(ddR[0], H(i, i) / 2.0, eps99, eps99);
          //
          // check off diagonal elements
          for(size_t j = 0; j < n; j++) if( j != i )
          {    x1[j] = 1.0;
               //
               store = true;
               R_kedem.Forward(1, x1, store);
               store = false;
               ddR = R_kedem.Forward(2, x2, store);
               double Hij = ( 2.0 * ddR[0] - H(i, i) - H(j, j) ) / 2.0;
               ok &= CppAD::NearEqual(Hij, H(i, j), eps99, eps99);
               //
               x1[j] = 0.0;
          }
          //
          x1[i] = 0.0;
     }
     return ok;
}