Make a more generic `eigen_gaussian` function

In ITensor/ITensors.jl#1106, @b-kloss introduced a function eigen_gaussian which properly diagonalizes fermionic Gaussian Hamiltonians and correlation matrices with pairing terms using a Schur decomposition. This leaves a situation where users should use Julia's eigen function from LinearAlgebra for number conserving fermionic Gaussian Hamiltonians and states and eigen_gaussian for more general parity conserving fermionic Gaussian Hamiltonians and states. It would be nice to develop a common interface for both of those cases.

One option is that internally we have wrappers ConservesNf and ConservesNfParity which can wrap any object (i.e. a Hamiltonian, correlation matrix, Gaussian circuit, etc.) and mark that they are either number conserving or parity conserving, which helps with code organization so that we can have code that is mostly generic across those two cases but then specialized with dispatch as needed. We could expose those types and then dispatch on them in eigen_gaussian (or eigen). Functions like hopping_hamiltonian(::OpSum) and quadratic_hamiltonian(::OpSum) could output quadratic Hamiltonians wrapped in ConservesNf and ConservesNfParity (respectively) which could then automatically get used for dispatch, so some of these details could get hidden from users.

Start by reviewing eigen_gaussian and Julia's LinearAlgebra.eigen, then compare the number-conserving and parity-conserving paths described in the issue. Examine hopping_hamiltonian(::OpSum) and quadratic_hamiltonian(::OpSum), along with the proposed ConservesNf and ConservesNfParity wrappers. Done means both cases share a coherent interface while retaining specialized dispatch where needed.