%%%%%%%%%%%%%%%%%%%
% XLiFE++ is an extended library of finite elements written in C++
%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%
 
\section{{\libtitle eigenCore} sub-library}

Most of the standard eigenvalue algorithms exploit projection processes in order to extract approximate eigenvectors from a given subspace. The basic idea of a projection method is to extract an approximate eigenvector from a specified low-dimensional subspace. If these approximation can be extracted, a small matrix eigenvalue problem is obtained. The numerical solution of the small $m\times m$ eigenvalue problem will be treated by the {\lib eigenCore} library.
The {\lib eigenCore} library is organized in multi-level, each of which utilizes subroutines provided by the lower one.
 \begin{itemize}
\item {\lib utils} : the lowest level provides various basic operations
\item {\lib decomposition} : the second level implements several decomposition algorithms
\item {\lib eigenSolver} : the highest level solves different eigen problem
\end{itemize}

\subsection{{\libtitle utils} part}

This part contains several classes which work on different basic operations.

\subsubsection{The {\classtitle VectorEigenDense} class}

Class {\class VectorEigenDense} is derived directly from {\class Vector}. Besides taking advantages of all funtionalities of {\class Vector}, this class supplements some particular functions to calculate Householder transformation. 
\vspace{.1cm}
\begin{lstlisting}
template<typename K>
class VectorEigenDense : public Vector<K>
{
public:
    typedef typename VectorEigenDense::type_t Scalar;
    typedef typename NumTraits<Scalar>::RealScalar RealScalar;
    typedef typename VectorEigenDense::it_vk it_vk;
    typedef typename VectorEigenDense::cit_vk cit_vk;

    VectorEigenDense() : Vector<K>(), acType_(_col) {}
    VectorEigenDense(const Dimen l) : Vector<K>(l, K()), acType_(_col) {}
    ...

\end{lstlisting}
\vspace{.1cm}
This class offers some additional functions to calculate House holder values.
\vspace{.1cm}
\begin{lstlisting}[]{}
void makeHouseHolderInPlace(K& tau, Real& beta);
void makeHouseHolderInPlace(K& tau, Real& beta, Number tail);
void makeHouseHolder(VectorEigenDense<K>& essential, K& tau, Real& beta) const;
\end{lstlisting}
\vspace{.2cm}
Note that the in the future, we may improve calculation performance of eigensolver for {\em dense} matrix by changing this class {\class VectorEigenDense}. By using the same block of allocated data, we can gain more efficiencies. However, for the moment, it's sufficient to have this class inherited from class {\class Vector} 
\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={Traits.hpp}}

\subsubsection{The {\classtitle MatrixEigenDense} class}

Class {\class MatrixEigenDense} inherits from the fundamental class {\class Matrix}. Basically, this class reuses all functionalities of its father in addition to some methods to support computing House holder, Cholesky composition and Jacobian rotation.
\vspace{.1cm}
\begin{lstlisting}
template<typename K>
class MatrixEigenDense : public Matrix<K>
{
public:
    typedef typename MatrixEigenDense<K>::type_t Scalar;
    typedef typename NumTraits<Scalar>::RealScalar RealScalar;
    typedef typename std::vector<K>::iterator it_vk;
    typedef typename std::vector<K>::const_iterator cit_vk;

    MatrixEigenDense() : Matrix<K>(), col_(1) {}
    MatrixEigenDense(const Dimen r) : Matrix<K>(r,r), col_(r) {}
    ...

\end{lstlisting}
\vspace{.1cm}
Householder computation can be invoked with
\vspace{.1cm}
\begin{lstlisting}[]{}
void applyHouseholderOnTheLeft(const VectorEigenDense<K>& essential, const K& tau);
void applyHouseholderOnTheRight(const VectorEigenDense<K>& essential, const K& tau);
\end{lstlisting}
\vspace{.1cm}
Cholesky decomposition can be calculated by
\vspace{.1cm}
\begin{lstlisting}[]{}
MatrixEigenDense<K> cholesky() const;
\end{lstlisting}
\vspace{.1cm}
or one can call Jacobian rotation with
\vspace{.1cm}
\begin{lstlisting}[]{}
template<typename OtherScalar>
inline void applyOnTheLeft(Index p, Index q, const JacobiRotation<OtherScalar>& j);
template<typename OtherScalar>
inline void applyOnTheRight(Index p, Index q, const JacobiRotation<OtherScalar>& j);
\end{lstlisting}
\vspace{.1cm}
Similar to class {\class VectorEigenDense},  in the future, this class may be improved to make {\lib eigenCore} more efficient.
\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={Traits.hpp}}

\subsubsection{The {\classtitle JacobiRotation} class}

Class {\class JacobiRotation} represents a {\em Jacobian} or {\em Givens} rotation,  2D rotation in a plane defined by its cosine and sine. 
\vspace{.1cm}
\begin{lstlisting}
template<typename Scalar>
class JacobiRotation
{
  public:
    typedef typename NumTraits<Scalar>::RealScalar RealScalar;

    /** Default constructor without any initialization. */
    JacobiRotation() {}

    /** Construct a planar rotation from a cosine-sine pair (\a c, \c s). */
    JacobiRotation(const Scalar& c, const Scalar& s) : c_(c), s_(s) {}

    Scalar& c() { return c_; }
    Scalar c() const { return c_; }
    Scalar& s() { return s_; }
    Scalar s() const { return s_; }
    ...
\end{lstlisting}
\vspace{.1cm}
The class provide two basic operations: making a Jacobi rotation and applying this rotation on both the right and left side of a selfadjoint $2\times2$ matrix; making a Givens rotation and applying this rotation on the left side of  a vector.
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y, z}]
template<typename Scalar>
bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z);
template<typename Scalar>
void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* z);
\end{lstlisting}
\vspace{.1cm}

\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={Traits.hpp}}

\subsubsection{The {\classtitle NumTraits} class}

{\class NumTraits} greatly facilitate the management of the sort of extra parameters that come up during the implementation of eigensolver algorithm.
\begin{lstlisting}
template<typename K>
struct NumTraits {};

template<>
struct NumTraits<Real> {
    enum {
      IsComplex = 0
    };
    static inline bool isReal() { return true;}
    static inline bool isComplex() { return false;}
    static inline Real zero() { return Real(0.0);}
    static inline Real one() { return Real(1.0);}
    static inline Real imag(const Real& v) { return v;}
    static inline Real real(const Real& v) { return v;}
    static inline Real value(const Real& v) { return v;}
    static inline Real abs2(const Real& v) { return v*v;}
    static inline Real epsilon() { return std::numeric_limits<Real>::epsilon(); }
    typedef Real RealScalar;
    typedef Complex ComplexScalar;
};

template<>
struct NumTraits<Complex> {
    enum {
      IsComplex = 1
    };
    static inline bool isReal() {return false;}
    static inline bool isComplex() { return true;}
    static inline Complex zero() { return Complex(0.0,0.0);}
    static inline Complex one() { return Complex(1.0,0.0);}
    static inline Real imag(const Complex& v) { return v.imag();}
    static inline Real real(const Complex& v) { return v.real();}
    static inline Complex value(const Complex& v) { return v;}
    static inline Real abs2(const Complex& v) { return v.real()*v.real() + v.imag()*v.imag();}
    typedef Real RealScalar;
    typedef Complex ComplexScalar;
};
\end{lstlisting}
\vspace{.1cm}
\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={utils.h}}

\subsection{{\libtitle decomposition} part}

All classes of this section utilize functions provided from the lower level - {\class utils} - to implement various decomposition algorithms, which play an important role in solving different eigen problems.

\subsubsection{The {\classtitle Tridiagonalization} class}

Class {\class Tridiagonalization} performs a tridiagonal decomposition of a self-adjoint matrix $ A$ such that:
  $A = Q T Q^* $ where $Q$ is unitary and $ T $ a real symmetric tridiagonal matrix. A tridiagonal matrix is a matrix which has non-zero elements only on the main diagonal and the first diagonal below and above it. This class is used in SelfAdjointEigenSolver to compute the eigenvalues and eigenvectors of a selfadjoint matrix.
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] size}]
template<typename _MatrixType>
class Tridiagonalization
{
  public:

    /** \brief Synonym for the template parameter \p _MatrixType. */
    typedef _MatrixType MatrixType;

    typedef typename MatrixType::type_t Scalar;
    typedef typename NumTraits<Scalar>::RealScalar RealScalar;
    typedef VectorEigenDense<Scalar> CoeffVectorType;
    typedef VectorEigenDense<RealScalar> DiagonalType;
    ...
    Tridiagonalization(Dimen size);
    Tridiagonalization(const MatrixType& matrix);
    ...	
\end{lstlisting}
\vspace{.1cm}
This class comes with default constructor with size of the matrix whose tridiagonal decomposition will be computed. In fact, this parameter is just a hint. After being created, object of this class can be used to calculate any square {\em dense} matrix with 
\vspace{.1cm}
\begin{lstlisting}[]{}
Tridiagonalization& compute(const MatrixType& matrix);
\end{lstlisting}
\vspace{.1cm}
The unitary matrix $Q$ can only be invoked after calling the {\itshape compute}:
\vspace{.1cm}
\begin{lstlisting}[]{}
MatrixType matrixQ();
\end{lstlisting}
\vspace{.1cm}
"Householder coefficients" can be retrieved by
\vspace{.1cm}
\begin{lstlisting}[]{}
inline CoeffVectorType householderCoefficients() const;
\end{lstlisting}
\vspace{.1cm}
We don't need the full-form of matrix $T$ in computation of eigenvalues; however its real diagonal and its real sub-diagonal are essential in calculating eigenvalues of the matrix. These two values can be retrieved with: 
\vspace{.1cm}
\begin{lstlisting}[]{}
DiagonalReturnType diagonal() const;
SubDiagonalReturnType subDiagonal() const;
\end{lstlisting}
\vspace{.1cm}
For example, the following $5\times 5$ real symmetric matrix 
\[
A=\left[
\begin{array}{cccccc}
1.36 &             -0.816  &             0.521  &              1.43 &             -0.144\\
-0.816 &             -0.659  &             0.794   &           -0.173  &            -0.406  \\
0.521  &            0.794 &             -0.541 &              0.461&              0.179  \\
1.43&              -0.173 &              0.461&               -1.43 &              0.822 \\
-0.144  &            -0.406 &              0.179 &              0.822&               -1.37
\end{array}
\right]
\]
after being tridiagonalized, has diagonal and sub-diagonal vectors\\
\mbox{}\hspace{10 mm} diagonal     = ( 1.36 -0.659 -0.541 -1.43 -1.37 ) \\
\mbox{}\hspace{10 mm} sub-diagonal     = ( -0.816 0.794 0.461 0.822 )

\displayInfos{library=eigenSolvers, test={test\_EigenSolverDecomposition.cpp}, header dep={Traits.hpp, VectorEigenDense.hpp, HouseHolderSequence.hpp}}

\subsubsection{The {\classtitle HessenbergDecomposition} class }

Class {\class HessenbergDecomposition} performs an Hessenberg decomposition of a matrix $ A $. In the real case, the Hessenberg decomposition consists of an orthogonal matrix $ Q $ and a Hessenberg matrix $ H $ such that $ A = Q H Q^T $. An orthogonal matrix is a matrix whose inverse equals its transpose ($ Q^{-1} = Q^T $). A Hessenberg matrix has zeros below the subdiagonal, so it is almost upper triangular. The Hessenberg decomposition of a complex matrix is $ A = Q H Q^* $ with $ Q $ unitary (that is, $ Q^{-1} = Q^* $). \\
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] size}]
template<typename _MatrixType>
class HessenbergDecomposition
{
  public:

    /** \brief Synonym for the template parameter \p _MatrixType. */
    typedef _MatrixType MatrixType;

    /** \brief Scalar type for matrices of type #MatrixType. */
    typedef typename MatrixType::type_t Scalar;

    /** \brief Type for vector of Householder coefficients.
      *
      * This is column vector with entries of type #Scalar. The length of the
      * vector is one less than the size of #MatrixType,
      */
    typedef VectorEigenDense<Scalar> CoeffVectorType;

    /** \brief Return type of matrixQ() */
    typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;

    /** \brief Default constructor; the decomposition will be computed later.
      *
      * \param [in] size  The size of the matrix whose Hessenberg decomposition will be computed.
      *
      * The default constructor is useful in cases in which the user intends to
      * perform decompositions via compute().  The \p size parameter is only
      * used as a hint. It is not an error to give a wrong \p size, but it may
      * impair performance.
      *
      * \sa compute() for an example.
      */
    HessenbergDecomposition(Number size)
      : matrix_(size,size),
        isInitialized_(false)
    {
      if(size>1)
        hCoeffs_.resize(size-1);
    }
    ...	
\end{lstlisting}
\vspace{.1cm}
We can use the constructor  which computes the Hessenberg decomposition at construction time.
\vspace{.1cm}
\begin{lstlisting}[]{}
HessenbergDecomposition(const MatrixType&);
\end{lstlisting}
\vspace{.1cm}
Once the decomposition is computed, we can use the {\cmd matrixH()} and {\cmd matrixQ()} functions to construct the matrices $H$ and $Q$ in the decomposition.
\vspace{.1cm}
\begin{lstlisting}[]{}
MatrixType matrixQ();
MatrixType matrixH() const;
\end{lstlisting}
\vspace{.1cm}
Alternatively, calling the function {\cmd compute()} to compute the Hessenberg decomposition of a given matrix. then we can use these two functions above
\vspace{.1cm}
\begin{lstlisting}[]{}
HessenbergDecomposition& compute(const MatrixType& matrix);
\end{lstlisting}
\vspace{.1cm}
For instance, the following $5\times 5$ real symmetric matrix 
\[
A=\left[
\begin{array}{cccccc}
1.36 &             -0.816  &             0.521  &              1.43 &             -0.144\\
-0.816 &             -0.659  &             0.794   &           -0.173  &            -0.406  \\
0.521  &            0.794 &             -0.541 &              0.461&              0.179  \\
1.43&              -0.173 &              0.461&               -1.43 &              0.822 \\
-0.144  &            -0.406 &              0.179 &              0.822&               -1.37
\end{array}
\right]
\]
has Hessenberg decomposition
\footnotesize
\[
H=\left[
\begin{array}{cccccc}
1.36 &       1.7328972849  &                 0 &                  0   &                0 \\
1.7328972849&       -1.1933420269&      -0.96533126903&   -9.1846580094e-17&     2.021583842e-16  \\
0 &     -0.96533126903&       -1.2812029745&       0.21380585197&    1.9428902931e-16  \\
0 &                  0  &     0.21380585197 &      -1.6885671008&       0.34719629211 \\
0 &  0 &                  0&       0.34719629211&       0.16311210222
\end{array}
\right]
\]
\normalsize
and unitary matrix
\footnotesize
\[
Q=\left[
\begin{array}{cccccc}
1 &                  0  &                 0   &                0   &                0\\
0&      -0.47088769029&       0.12629876626&      -0.67227759101&      -0.55709626223  \\
0&       0.30065255716&      -0.19453390601&       0.43672469898 &     -0.82524913607 \\
0&       0.82520759451&      0.045097442878&      -0.56297063107&    -0.0079192904815 \\
0 &    -0.083097827699&      -0.97168482632&      -0.20094388887&      0.092438643767
\end{array}
\right]
\]
\normalsize

\displayInfos{library=eigenSolvers, test={test\_EigenSolverDecomposition.cpp}, header dep={Traits.hpp, VectorEigenDense.hpp}}

\subsubsection{The {\classtitle RealSchur} class }

Class {\class RealSchur} performs a real Schur decomposition of a square matrix. Given a real square matrix $A$, this class computes the real Schur decomposition: $ A = U T U^T $ where $U$ is a real orthogonal matrix and $T$ is a real quasi-triangular matrix. An orthogonal matrix is a matrix whose inverse is equal to its transpose, $ U^{-1} = U^T $. A quasi-triangular matrix is a block-triangular matrix whose diagonal consists of 1-by-1 blocks and 2-by-2 blocks with complex eigenvalues. The eigenvalues of the blocks on the diagonal of $T$ are the same as the eigenvalues of the matrix $A$, and thus the real Schur decomposition is used in {\class RealEigenSolver} to compute the eigendecomposition of a matrix.
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] size}]
template<typename _MatrixType>
class RealSchur
{
  public:
    typedef _MatrixType MatrixType;
    typedef typename MatrixType::type_t Scalar;
    typedef typename NumTraits<Scalar>::ComplexScalar ComplexScalar;

    typedef VectorEigenDense<ComplexScalar> EigenvalueType;
    typedef VectorEigenDense<Scalar> ColumnVectorType;

    /** \brief Default constructor.
      *
      * \param [in] size  Positive integer, size of the matrix whose Schur decomposition will be computed.
      *
      * The default constructor is useful in cases in which the user intends to
      * perform decompositions via compute().  The \p size parameter is only
      * used as a hint. It is not an error to give a wrong \p size, but it may
      * impair performance.
      *
      * \sa compute() for an example.
      */
    RealSchur(Index size)
            : matT_(size, size),
              matU_(size, size),
              hess_(size),
              maxIterations_(40),
              isInitialized_(false),
              matUisUptodate_(false),
              info_(_noConvergence)
    { }
    ...	
\end{lstlisting}
\vspace{.1cm}
We can use the constructor which computes the real Schur decomposition at construction time. 
\vspace{.1cm}
\begin{lstlisting}[]{}
RealSchur(const MatrixType& matrix, bool computeU = true);
\end{lstlisting}
\vspace{.1cm}
Once the decomposition is computed, we can use the {\cmd matrixU()} and {\cmd matrixT()} functions to retrieve the matrices $U$ and $T$ in the decomposition.
\vspace{.1cm}
\begin{lstlisting}[]{}
MatrixType matrixU() const;
MatrixType matrixT() const;
\end{lstlisting}
\vspace{.1cm}
Alternatively calling the function {\cmd compute()} to compute the real Schur decomposition of a given matrix, after that we are able to use these two functions above. \\
The implementation of Schur decomposition of real square matrix is adapted from \url{http://math.nist.gov/javanumerics/jama/}, whose code is based on \eispack. The Schur decomposition is computed by first reducing the matrix to $Hessenberg$ form using the class {\class HessenbergDecomposition}. The $Hessenberg$ matrix is then reduced to triangular form by performing $Francis QR$ iterations with implicit double shift.\\
There are several auxiliary functions in order to implement this algorithm :
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] norm}]
Scalar computeNormOfT();
Index findSmallSubdiagEntry(Index iu, Scalar norm);
void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector);
\end{lstlisting}
\vspace{.1cm}
For instance, the following $5\times 5$ non symmetric matrix 
\[
A=\left[
\begin{array}{cccccc}
 1.8779&             -0.5583&             -0.2099&              0.2696&              0.1097\\
0.9407&             -0.3114&             -1.6989&              0.4943&              1.1287  \\
0.7873&               -0.57&              0.6076&             -1.4831&               -0.29  \\
-0.8759&             -1.0257&             -0.1178&             -1.0203&              1.2616 \\
0.3199&             -0.9087&              0.6992&              -0.447&              0.4754
\end{array}
\right]
\]
has the quasi-triangular matrix in the Schur decomposition
\footnotesize
\[
H=\left[
\begin{array}{cccccc}
1.6552531332&      -0.30767059895 &     -0.45435806468  &      1.5602779282 &      0.68358491675\\ 
                    0 &       1.1891303949 &       1.9306398993 &      0.96785347442 &      0.57970343469 \\
                    0 &      -1.3969331077&       0.22387492903 &     -0.78380738142&          1.27402357 \\
                    0 &                  0  &                 0 &     0.092450608769&       -0.5563858236 \\
                    0  &                 0 &                  0 &                  0 &      -1.5315090659
\end{array}
\right]
\]
\normalsize
and the orthogonal matrix in the Schur decomposition.
\footnotesize
\[
Q=\left[
\begin{array}{cccccc}
0.72665944996&       0.28908942548&     -0.090578168382&       0.61658006012&     0.0042394246225 \\
        0.13365823926&       0.75252664471 &    -0.11720311204 &    -0.52512343394&      -0.35545454525 \\
        0.50462706182 &     -0.59105620177&      -0.20754951129&      -0.34476029525&       -0.4838053679 \\
        -0.1068784212 &     0.010284495663 &        0.747282602 &      0.23512571367&      -0.61217305049 \\
        0.43362431945&     -0.026033604674&       0.61363746549 &     -0.41222645012&      0.51459752796
\end{array}
\right]
\]
\normalsize
\displayInfos{library=eigenSolvers, test={test\_EigenSolverDecomposition.cpp}, header dep={Jacobi.hpp, VectorEigenDense.hpp, HessenbergDecomposition.hpp}}

\subsection{{\libtitle eigenSolver} main part}

In this section, there are solvers for three kinds of eigen problem :
\begin{itemize}
\item Hermitian eigenproblem 
\item Generalized Hermitian eigenproblem 
\item Non-Hermitian (real) eigenproblem
\end{itemize}
The following presents eigen solver classes corresponding to each kind of eigen problem.

\subsubsection{The {\classtitle SelfAdjointEigenSolver} class}

Class {\class SelfAdjointEigenSolver} computes eigenvalues and eigenvectors of self-adjoint matrices.
A matrix $A$ is selfadjoint if it equals its adjoint. For real matrices, this means that the matrix is symmetric: it equals its transpose. \\
This class computes the eigenvalues and eigenvectors of a self-adjoint matrix. These are the scalars $ \lambda $ and vectors $ v $ such that $ Av = \lambda v $.  The eigenvalues of a self-adjoint matrix are always real. If $ D $ is a diagonal matrix with the eigenvalues on the diagonal, and $ V $ is a matrix with the eigenvectors as its columns, then $ A = V D V^{-1} $ (for self-adjoint matrices, the matrix $ V $ is always invertible). This is called the eigendecomposition. The implemented algorithm exploits the fact that the matrix is self-adjoint, making it faster and more accurate than the general purpose eigenvalue algorithms implemented in RealEigenSolver.\\
\vspace{.1cm}
\begin{lstlisting}
template<typename _MatrixType>
class SelfAdjointEigenSolver
{
  public:

    typedef _MatrixType MatrixType;
    
    /** \brief Scalar type for matrices of type \p _MatrixType. */
    typedef typename MatrixType::type_t Scalar;

    /** \brief Real scalar type for \p _MatrixType.
      *
      * This is just \c Scalar if #Scalar is real (e.g., \c float or
      * \c double), and the type of the real part of \c Scalar if #Scalar is
      * complex.
      */
    typedef typename NumTraits<Scalar>::RealScalar RealScalar;
    
    /** \brief Type for vector of eigenvalues as returned by eigenvalues().
      *
      * This is a column vector with entries of type #RealScalar.
      * The length of the vector is the size of \p _MatrixType.
      */
    typedef VectorEigenDense<RealScalar> RealVectorType;
    typedef Tridiagonalization<MatrixType> TridiagonalizationType;

    /** \brief Default constructor for fixed-size matrices.
      *
      * The default constructor is useful in cases in which the user intends to
      * perform decompositions via compute().
      *
      */
    SelfAdjointEigenSolver()
        : eivec_(),
          eivalues_(),
          subdiag_(),
          maxIterations_(30),
          isInitialized_(false)
    { }
    ...	
\end{lstlisting}
\vspace{.1cm}
We can use the constructor which computes the eigenvalues and eigenvectors at construction time. 
\vspace{.1cm}
\begin{lstlisting}[]{}
SelfAdjointEigenSolver(const MatrixType& matrix, Int options = _computeEigenVector);
\end{lstlisting}
\vspace{.1cm}
Once the eigenvalue and eigenvectors are computed, they can be retrieved with 
\vspace{.1cm}
\begin{lstlisting}[]{}
const RealVectorType& eigenvalues() const;
const MatrixType& eigenvectors() const;
\end{lstlisting}
\vspace{.1cm}
Alternatively, we can call the function {\cmd compute()} to compute the eigenvalues and eigenvectors of a given matrix.
\vspace{.1cm}
\begin{lstlisting}[]{}
SelfAdjointEigenSolver& compute(const MatrixType& matrix, Int options = _computeEigenVector);
\end{lstlisting}
\vspace{.1cm}
This implementation uses a symmetric QR algorithm. The matrix is first reduced to tridiagonal form using the {\class Tridiagonalization} class. The tridiagonal matrix is then brought to diagonal form with implicit symmetric QR steps with Wilkinson shift.\\
For instance, the following $5\times 5$ real symmetric matrix 
\[
A=\left[
\begin{array}{cccccc}
1.36 &             -0.816  &             0.521  &              1.43 &             -0.144\\
-0.816 &             -0.659  &             0.794   &           -0.173  &            -0.406  \\
0.521  &            0.794 &             -0.541 &              0.461&              0.179  \\
1.43&              -0.173 &              0.461&               -1.43 &              0.822 \\
-0.144  &            -0.406 &              0.179 &              0.822&               -1.37
\end{array}
\right]
\]
has eigenvalues and corresponding (column) eigenvector\\
\mbox{}\hspace{10 mm} eigenvalues     = ( -2.64925 -1.76682 -0.742793 0.22723 2.29163 ) \\
\footnotesize
\[
\mbox{\normalsize eigenvectors =} \left[
\begin{array}{cccccc}
-0.32648814917&     -0.098193931162&       0.34709885718&     -0.010999922777&       0.87359305479 \\
-0.20806598152&     -0.64229394812&      0.22733262857 &     0.66233712946 &     -0.23194058392 \\
0.050094522003 &      0.62960023585&      -0.16376394055 &      0.73987891396&      0.16387387324 \\
0.72069515152 &     -0.39678270513  &    -0.40156098976&       0.11452959441&       0.38573789934 \\
-0.57289010572&      -0.15486595502&      -0.79985773898&     -0.025508690809 &     0.085967242812
\end{array}
\right]
\]
\normalsize

\displayInfos{library=eigenSolvers, test={test\_EigenSolverInternCore.cpp}, header dep={VectorEigenDense.hpp, Tridiagonalization.hpp}}

\subsubsection{The {\classtitle GeneralizedSelfAdjointEigenSolver} class}

Class {\class GeneralizedSelfAdjointEigenSolver} computes eigenvalues and eigenvectors of the generalized self-adjoint eigen problem. \\ This class solves the generalized eigenvalue problem $ Av = \lambda Bv $. In this case, the matrix $ A $ should be self-adjoint and the matrix $ B $ should be positive definite.
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] size}]
template<typename _MatrixType>
class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
{
    typedef SelfAdjointEigenSolver<_MatrixType> Base;
public:
    typedef _MatrixType MatrixType;

    /** \brief Default constructor for fixed-size matrices.
      *
      * The default constructor is useful in cases in which the user intends to
      * perform decompositions via compute().
      */
    GeneralizedSelfAdjointEigenSolver() : Base() {}

    /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
      *
      * \param [in]  size  Positive integer, size of the matrix whose
      * eigenvalues and eigenvectors will be computed.
      *
      * \sa compute() for an example
      */
    GeneralizedSelfAdjointEigenSolver(Index size)
        : Base(size)
    {}
    ...	
\end{lstlisting}
\vspace{.1cm}
We can use the constructor which computes the eigenvalues and eigenvectors at construction time. 
\vspace{.1cm}
\begin{lstlisting}[]{}
    GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
                                      int options = _computeEigenVector|_Ax_lBx)
      : Base(matA.numOfCols());
\end{lstlisting}
\vspace{.1cm}
Alternatively, we can call the function {\itshape compute()} to compute the eigenvalues and eigenvectors of a given matrix.
\vspace{.1cm}
\begin{lstlisting}[]{}
    GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
                                               int options = _computeEigenVector|_Ax_lBx);
\end{lstlisting}
\vspace{.1cm}
Eigenvalues and eigenvectors can be retrieved with same functions of {\class SelfAdjointEigenSolver}
\vspace{.1cm}
\begin{lstlisting}[]{}
const RealVectorType& eigenvalues() const;
const MatrixType& eigenvectors() const;
\end{lstlisting}
\vspace{.1cm}
For instance, the following $5\times 5$ real symmetric matrix $A$ and positive definite matrix $B$
\[
A=\left[
\begin{array}{cccccc}
1.36 &             -0.816  &             0.521  &              1.43 &             -0.144\\
-0.816 &             -0.659  &             0.794   &           -0.173  &            -0.406  \\
0.521  &            0.794 &             -0.541 &              0.461&              0.179  \\
1.43&              -0.173 &              0.461&               -1.43 &              0.822 \\
-0.144  &            -0.406 &              0.179 &              0.822&               -1.37
\end{array}
\right]
\]
\[
B=\left[
\begin{array}{cccccc}
1.1981& 1.4405& 1.5367& 1.0971& 1.1276\\
1.4405& 2.1207& 1.9115& 1.3836& 1.642\\
1.5367& 1.9115& 2.4828& 1.7702& 1.7543\\
1.0971& 1.3836& 1.7702& 1.4359& 1.2428\\
1.1276& 1.642& 1.7543& 1.2428& 1.424
\end{array}
\right]
\]
has eigenvalues and corresponding (column) eigenvector\\
\mbox{}\hspace{10 mm} eigenvalues     = ( -40.0667 -12.6587 -3.08602 0.0522151 40.3448 ) \\
\footnotesize
\[
\mbox{\normalsize eigenvectors =} \left[
\begin{array}{cccccc}
-2.5241477328&      -0.01326907506 &      0.83333293087 &    -0.040348533316 &       4.3939678718\\ 
3.7607674197  &      0.3488243247  &      1.1073133113   &    -0.3174355818   &    -4.3303408407 \\
3.0396247156   &      2.615639864  &    -0.34676637691   &   -0.36756457169   &    -3.9258019948 \\
0.45843793898  &     -2.0051996759  &    -0.65863554474 &    -0.039469282032   &     1.4762595546 \\
-6.6376554636   &    -1.9381010232  &    -0.97022016252   &   0.070809133211 &       4.9718019236
\end{array}
\right]
\]
\normalsize
\displayInfos{library=eigenSolvers, test={test\_EigenSolverInternCore.cpp}, header dep={Tridiagonalization.hpp}}

\subsubsection{The {\classtitle RealEigenSolver} class}

Class {\class RealEigenSolver} computes eigenvalues and eigenvectors of general real matrices. \\
The eigenvalues and eigenvectors of a matrix $ A $ are scalars $ \lambda $ and vectors $ v $ such that $ Av = \lambda v $.  If $ D $ is a diagonal matrix with the eigenvalues on the diagonal, and $ V $ is a matrix with the eigenvectors as its columns, then $ A V = V D $. The matrix $ V $ is almost always invertible, in which case we have $ A = V D V^{-1} $. \\
The eigenvalues and eigenvectors of the matrix may be complex, even when the matrix is real. However, we can choose real matrices $ V $ and $ D $ satisfying $ A V = V D $, just like the eigendecomposition, if the matrix $ D $ is not required to be diagonal, but if it is allowed to have blocks of the form\[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \] (where $ u $ and $ v $ are real numbers) on the diagonal.  These blocks correspond to complex eigenvalue pairs $ u \pm iv $. We call this variant of the eigendecomposition the pseudo-eigendecomposition.
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] size}]
template<typename _MatrixType>
class RealEigenSolver
{
  public:

    /** \brief Synonym for the template parameter \p _MatrixType. */
    typedef _MatrixType MatrixType;

    /** \brief Scalar type for matrices of type #MatrixType. */
    typedef typename MatrixType::type_t Scalar;
    typedef typename NumTraits<Scalar>::RealScalar RealScalar;
    typedef VectorEigenDense<Scalar> ColumnVectorType;
    typedef VectorEigenDense<Scalar> VectorType;

    /** \brief Complex scalar type for #MatrixType. 
      *
      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
      * \c float or \c double) and just \c Scalar if #Scalar is
      * complex.
      */
    typedef typename NumTraits<Scalar>::ComplexScalar ComplexScalar;
    typedef VectorEigenDense<ComplexScalar> ComplexVectorType;

    /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 
      *
      * This is a column vector with entries of type #ComplexScalar.
      * The length of the vector is the size of #MatrixType.
      */
    typedef VectorEigenDense<ComplexScalar> EigenvalueType;

    /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 
      *
      * This is a square matrix with entries of type #ComplexScalar. 
      * The size is the same as the size of #MatrixType.
      */
    typedef MatrixEigenDense<ComplexScalar> EigenvectorsType;

    /** \brief Default constructor.
      *
      * The default constructor is useful in cases in which the user intends to
      * perform decompositions via RealEigenSolver::compute(const MatrixType&, bool).
      *
      * \sa compute() for an example.
      */
 RealEigenSolver() : eivec_(), eivalues_(), isInitialized_(false), realSchur_(), matT_() {}

    /**
     * \brief Default constructor with memory preallocation
      *
      * Like the default constructor but with preallocation of the internal data
      * according to the specified problem \a size.
      * \sa RealEigenSolver()
      */
    RealEigenSolver(Index size)
      : eivec_(size, size),
        eivalues_(size),
        isInitialized_(false),
        eigenvectorsOk_(false),
        realSchur_(size),
        matT_(size, size)
    {}
    ...	
\end{lstlisting}
\vspace{.1cm}
Same as other two eigensolvers, we can use the constructor which computes the eigenvalues and eigenvectors at construction time. 
\vspace{.1cm}
\begin{lstlisting}[]{}
RealEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true);
\end{lstlisting}
\vspace{.1cm}
Once the eigenvalue and eigenvectors are computed, they can be retrieved with 
\vspace{.1cm}
\begin{lstlisting}[]{}
const EigenvalueType& eigenvalues() const;
EigenvectorsType eigenvectors() const;
\end{lstlisting}
\vspace{.1cm}
Alternatively, we can call the function {\cmd compute()} to compute the eigenvalues and eigenvectors of a given matrix.
\vspace{.1cm}
\begin{lstlisting}[]{}
RealEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
\end{lstlisting}
\vspace{.1cm}
This implementation first reduces the matrix to real Schur form using the {\class RealSchur} class. The Schur decomposition is then used to compute the eigenvalues and eigenvectors.\\
For instance, the following $5\times 5$ non symmetric matrix 
\[
A=\left[
\begin{array}{cccccc}
1.8779&             -0.5583&             -0.2099&              0.2696&              0.1097\\
0.9407&             -0.3114&             -1.6989&              0.4943&              1.1287  \\
0.7873&               -0.57&              0.6076&             -1.4831&               -0.29  \\
-0.8759&             -1.0257&             -0.1178&             -1.0203&              1.2616 \\
0.3199&             -0.9087&              0.6992&              -0.447&              0.4754
\end{array}
\right]
\]
has eigenvalues and corresponding (column) eigenvector\\
\mbox{}\hspace{10 mm} eigenvalues     = ( (1.65525,0) (0.706503,1.56973) (0.706503,-1.56973) (0.0924506,0) (-1.53151,0) ) \\
\[
\mbox{eigenvectors =} \left[
\begin{array}{ccccc}
(0.727,0)& (0.177,-0.001)& (0.177,0.001)& (-0.230,0)& (-0.041,0) \\
(0.133,0)& (0.534,-0.218)& (0.534,0.218)& (-0.664,0)& (-0.392,0) \\
(0.504,0)& (-0.444,0.0814) & (-0.444,-0.0814)& (-0.336,0)& (-0.587,0)\\ 
(-0.106,0)& (0.011,0.457)&  (0.0118,-0.457)& (0.135,0)& (-0.695,0) \\
(0.433,0)& (-0.036,0.469)& (-0.036,-0.469)& (-0.611,0)& (-0.121,0)
\end{array}
\right]
\]

\displayInfos{library=eigenSolvers, test={test\_EigenSolverInternCore.cpp}, header dep={VectorEigenDense.hpp, MatrixEigenDense.hpp, RealSchur.hpp}}