Speedup symmetric eigendecomposition

Description

Symmetric eigendecomposition from library Eigen is implemented in a suboptimal way. It consists of two steps: tridiagonalization and tridiagonal eigensolver.

In Eigen, basic Householder tridiagonalization is implemented. This can be improved on using blocked Householder algorithm, which uses cache more efficiently and is therefore faster.

For symmetric eigensolver Eigen uses QR iteration. While very accurate, it is quite slow. For square symmetric matrix of size n it does O(n^3) operations. It can be replaced by, for example, the
MRRR algorithm that is significantly faster - O(n^2) and usually accurate enough. However, for cases where multiple eigenvalues are very close to each other - with relative difference (l1-l2)/l1 between
them close to machine precision - MRRR might produce more significant errors or in extreme cases even fail to converge.

Both steps can be further improved using GPUs.

Example

Function is added that computes symmetric eigendecomposition.

Expected Output

Symmetric eigendecomposition is faster than one implemented in Eigen.

Current Version:

v2.18.1

The issue does not identify files, tests, or an entry point. Start by locating the current Eigen symmetric eigendecomposition implementation and any relevant benchmarks, then assess blocked Householder tridiagonalization and an MRRR eigensolver as alternatives. Done means faster symmetric eigendecomposition while preserving acceptable accuracy and convergence.