# include <cppad/cppad.hpp>
bool sparse_sub_hes(void)
{     bool ok = true;
     using CppAD::AD;
     typedef CPPAD_TESTVECTOR(size_t)     SizeVector;
     typedef CPPAD_TESTVECTOR(double)     DoubleVector;
     typedef CppAD::sparse_rc<SizeVector> sparsity;
     //
     // domain space vector
     size_t n = 4;
     CPPAD_TESTVECTOR(AD<double>) ax(n);
     for(size_t j = 0; j < n; j++)
          ax[j] = double(j);

     // declare independent variables and start recording
     CppAD::Independent(ax);

     // range space vector
     size_t m = 1;
     CPPAD_TESTVECTOR(AD<double>) ay(m);
     ay[0] = 0.0;
     for(size_t j = 0; j < n; j++)
          ay[0] += double(j+1) * ax[0] * ax[j];

     // create f: x -> y and stop tape recording
     CppAD::ADFun<double> f(ax, ay);

     // sparsity pattern for the identity matrix
     size_t nr     = n;
     size_t nc     = n;
     size_t nnz_in = n;
     sparsity pattern_in(nr, nc, nnz_in);
     for(size_t k = 0; k < nnz_in; k++)
     {     size_t r = k;
          size_t c = k;
          pattern_in.set(k, r, c);
     }
     // compute sparsity pattern for J(x) = f'(x)
     bool transpose       = false;
     bool dependency      = false;
     bool internal_bool   = false;
     sparsity pattern_out;
     f.for_jac_sparsity(
          pattern_in, transpose, dependency, internal_bool, pattern_out
     );
     //
     // compute sparsity pattern for H(x) = f''(x)
     CPPAD_TESTVECTOR(bool) select_range(m);
     select_range[0]      = true;
     CppAD::sparse_hes_work work;
     f.rev_hes_sparsity(
          select_range, transpose, internal_bool, pattern_out
     );
     size_t nnz = pattern_out.nnz();
     ok        &= nnz == 7;
     ok        &= pattern_out.nr() == n;
     ok        &= pattern_out.nc() == n;
     {     // check results
          const SizeVector& row( pattern_out.row() );
          const SizeVector& col( pattern_out.col() );
          SizeVector row_major = pattern_out.row_major();
          //
          ok &= row[ row_major[0] ] ==  0  && col[ row_major[0] ] ==  0;
          ok &= row[ row_major[1] ] ==  0  && col[ row_major[1] ] ==  1;
          ok &= row[ row_major[2] ] ==  0  && col[ row_major[2] ] ==  2;
          ok &= row[ row_major[3] ] ==  0  && col[ row_major[3] ] ==  3;
          //
          ok &= row[ row_major[4] ] ==  1  && col[ row_major[4] ] ==  0;
          ok &= row[ row_major[5] ] ==  2  && col[ row_major[5] ] ==  0;
          ok &= row[ row_major[6] ] ==  3  && col[ row_major[6] ] ==  0;
     }
     //
     // Only interested in cross-terms. Since we are not computing rwo 0,
     // we do not need sparsity entries in row 0.
     CppAD::sparse_rc<SizeVector> subset_pattern(n, n, 3);
     for(size_t k = 0; k < 3; k++)
          subset_pattern.set(k, k+1, 0);
     CppAD::sparse_rcv<SizeVector, DoubleVector> subset( subset_pattern );
     //
     // argument and weight values for computation
     CPPAD_TESTVECTOR(double) x(n), w(m);
     for(size_t j = 0; j < n; j++)
          x[j] = double(n) / double(j+1);
     w[0] = 1.0;
     //
     std::string coloring = "cppad.general";
     size_t n_sweep = f.sparse_hes(
          x, w, subset, subset_pattern, coloring, work
     );
     ok &= n_sweep == 1;
     for(size_t k = 0; k < 3; k++)
     {     size_t i = k + 1;
          ok &= subset.val()[k] == double(i + 1);
     }
     //
     // convert subset from lower triangular to upper triangular
     for(size_t k = 0; k < 3; k++)
          subset_pattern.set(k, 0, k+1);
     subset = CppAD::sparse_rcv<SizeVector, DoubleVector>( subset_pattern );
     //
     // This will require more work because the Hessian is computed
     // column by column (not row by row).
     work.clear();
     n_sweep = f.sparse_hes(
          x, w, subset, subset_pattern, coloring, work
     );
     ok &= n_sweep == 3;
     //
     // but it will get the right answer
     for(size_t k = 0; k < 3; k++)
     {     size_t i = k + 1;
          ok &= subset.val()[k] == double(i + 1);
     }
     return ok;
}