In my application I need the real Schur form of a matrix, as well as its eigenvalues (as provided by Lapack's DGEES).
Eigen's RealSchur will give me the Schur matrices U and T, but does not offer an interface to subsequently extract the eigenvalues from the 1x1 or 2x2 diagonal blocks of T. Such interface would be highly appreciated.
It seems that this would merely require a small factorization of the member compute() of the EigenSolver template: the part below the comment:
// Compute eigenvalues from matT
in EigenSolver.h is exactly what I need --- but rather not duplicate.
IMHO, the cleanest solution is to move this code to a new public member RealSchur<>::eigenvalues(), that is also used by EigenSolver<>::compute().
I see 3 options:
1 - call it RealSchur::eigenvalues(), and compute them at the end of RealSchur::compute (), similar to EigenSolver. Perhaps, with an option to compute() to request them?
2 - call it RealSchur::eigenvalues(), and compute them in a cached vector at the first call to eigenvalues(). The cache would be invalidated at every call to compute().
3 - no cache, and thus we would call it extractEigenvalues() to make it clear that each call to this function involves some computations.
To be more consistent with the other classes, I would go for either 1 or 2.
Created attachment 843 [details]
This patch implements adds member eigenvalues() to the RealSchur template:
* new private member m_eigenvalues hosts the eigenvalues
of the matrix; these are calculated as part of the factorization.
This vector can be retrieved through accessor eigenvalues();
* new private member extractEigenvalues() resizes m_eigenvalues
and extracts the eigenvalues from a previously calculated m_T.
This member returns NumericalIssue if a !isfinite eigenvalue
is encountered, Success otherwise;
* computeFromHessenberg() calls extractEigenvalues() after doing the
factorization. If the factorization generated an error, that
error is kept for later inspection (Convergence failure, ...),
otherwise m_info is overwritten by the result of extractEigenvalues();
* compute() calls computeFromHessenberg() on a scaled matrix (as before).
After that function returns, the (inverse) scaling factor is applied to
m_eigenvalues (as it is to m_T itself);
* in the special case of a near-zero matrix m_eigenvalues if manually
set to zero (as is m_T).
Even if eigenvalues() is not called, there is one user-visible change:
m_info is set to NumericalIssue if non-finite eigenvalues are encountered;
before these went undetected.
This patch does not contain additional tests. This new functionality will be
indirectly tested by the EigenSolver tests once EigenSolver starts to use
RealSchur::eigenvalues(). (That is part of a patch 2/2).
Created attachment 844 [details]
This patch changes EigenSolver to use RealSchur::eigenvalues() to obtain the eigenvalues.
The two patches implement Gaels' proposition #1. I checked that all eigensolver_generic* and real_schur* tests pass after applying these patches.
Preparing the patches cost me a surprising amount of time, mostly because the intended error checking semantics of EigenSolver are not always clear (to me). As an example, let us assume that the computation of the Schur decomposition fails:
* the EigenSolver object is marked to be initialized, although not all matrices are set (or even resized in accordance with the size of the incoming matrix).
* m_eigenvectorsOk is set to computeEigenVectors, while m_matU is not calculated (or even resized) in case of a Schur failure.
I am happy to provide a documentation or code fix, but need some guidance here: what is the intended state of the RealSchur and EigenSolver objects after compute is called but does not return successfully? Should matrixT() and eigenvalues() (and maybe matrixU()) be resized at least? Or are no guarantees provided whatsoever in that case, and is it possible that the sizes are even inconsistent (with each other or with the matrix that was passed to compute)?