Linear method eigensolver improvements

This issue is exploration into accelerating the eigensolver part of the linear method.
See #3633 for the umbrella issue with scaling the linear method with respect to the number of parameters.

The nonsymmetric generalized eigenvalue problem scales as N_p^3, so it will eventually dominate the time in each linear method step as the number of parameters increases.

This issue is to record some of the investigations and observations.

The problem is there are multiple orthogonal directions to accelerate: multicore, GPUs, and distributed memory. As far as I know, there are no good solutions that combine all of these.

Options
* LAPACK (invert overlap matrix + dgeev) - the current solution.
* Solve occurs on each MPI rank independently. Uses redundant computation to avoid communication.
* Vendor optimized Lapack (MKL, LibSci, etc). (mostly multicore, maybe some GPU)
* May be easy to use, but not a general solution. Might require some rework to enable threaded versions of the library.
* MAGMA (GPU acceleration) (https://icl.cs.utk.edu/magma/ )
* Provides a fairly easy boost in performance for runs that are already using a GPU. However, this is only a one-time boost.
* StarNEig (https://nlafet.github.io/StarNEig )
* In theory could use shared memory, GPU, and distributed memory. Unfortunately, eigenvectors are only implemented for shared memory.
* Scalapack (distributed memory)
* Distributed resources are likely plentiful as most runs are concerned with scaling the problem out to generate samples.
* Performs poorly on one rank (in one case, it's 5x slower than using MKL lapack directly)

The next step is to implement a Scalapack version of the generalized eigenvalue kernel and explore the scaling.