Sparsity patterns for Jacobian and Hessian computations

[a-jp]

Based on clad discussion #780. We assume the current computations of a Jacobian and Hessian are returned as a dense matrix. We wish to augment the API to allow sparse representations of the Jacobian and Hessian to be computed. This places two requirements:

1. The calculation should return to the user two pieces of information a) in an as yet to be specified sparse data structure, the sparsity structure (the indices of the entries that are non-zero) and b) the values of the non-zero entries (see 2) in a compressed format.
2. The non-zero entries computed and returned should be structurally non-zero as a function of the differentiation, and not zero due to initial values or constants. Said another way, the zero entries that are excluded must be defined to be always zero, and not just zero due to specific initial values.

As a check on the results outputted from 1 and 2, one should be able to reconstruct a dense matrix from the two pieces of information provided by 1. with additional zeros where needed. For a sanity check, one would then expect the recombination of the sparsity structure, and the compressed non-zero entries into a dense matrix to be identical (within machine precision) to the currently computed dense matrices computed by calls to the Jacobian and Hessian routines.

It's for library developers to comment on the best mechanics for integration into clad, a small amount of the discussion is contained in the initial discussion #780. It's possible that something is to be learned from other works such as CasADi although I neither use or have any affiliation with that codebase. Likewise, sparse Hessians and Jacobians are well accommodated for, as noted in the original discussion, in CppAD.