don't index Eigen

14 files changed, 0 insertions, 5484 deletions

diff --git a/eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h
deleted file mode 100644
index dc5fae0..0000000
--- a/eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ /dev/null
@@ -1,346 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2009 Claire Maurice
-// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
-#define EIGEN_COMPLEX_EIGEN_SOLVER_H
-
-#include "./ComplexSchur.h"
-
-namespace Eigen { 
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class ComplexEigenSolver
-  *
-  * \brief Computes eigenvalues and eigenvectors of general complex matrices
-  *
-  * \tparam _MatrixType the type of the matrix of which we are
-  * computing the eigendecomposition; this is expected to be an
-  * instantiation of the Matrix class template.
-  *
-  * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
-  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v
-  * \f$.  If \f$ D \f$ is a diagonal matrix with the eigenvalues on
-  * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as
-  * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is
-  * almost always invertible, in which case we have \f$ A = V D V^{-1}
-  * \f$. This is called the eigendecomposition.
-  *
-  * The main function in this class is compute(), which computes the
-  * eigenvalues and eigenvectors of a given function. The
-  * documentation for that function contains an example showing the
-  * main features of the class.
-  *
-  * \sa class EigenSolver, class SelfAdjointEigenSolver
-  */
-template<typename _MatrixType> class ComplexEigenSolver
-{
-  public:
-
-    /** \brief Synonym for the template parameter \p _MatrixType. */
-    typedef _MatrixType MatrixType;
-
-    enum {
-      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-      Options = MatrixType::Options,
-      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-    };
-
-    /** \brief Scalar type for matrices of type #MatrixType. */
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-    /** \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 std::complex<RealScalar> ComplexScalar;
-
-    /** \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 Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> 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 Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
-
-    /** \brief Default constructor.
-      *
-      * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via compute().
-      */
-    ComplexEigenSolver()
-            : m_eivec(),
-              m_eivalues(),
-              m_schur(),
-              m_isInitialized(false),
-              m_eigenvectorsOk(false),
-              m_matX()
-    {}
-
-    /** \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 ComplexEigenSolver()
-      */
-    explicit ComplexEigenSolver(Index size)
-            : m_eivec(size, size),
-              m_eivalues(size),
-              m_schur(size),
-              m_isInitialized(false),
-              m_eigenvectorsOk(false),
-              m_matX(size, size)
-    {}
-
-    /** \brief Constructor; computes eigendecomposition of given matrix.
-      *
-      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
-      *    eigenvalues are computed; if false, only the eigenvalues are
-      *    computed.
-      *
-      * This constructor calls compute() to compute the eigendecomposition.
-      */
-    template<typename InputType>
-    explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
-            : m_eivec(matrix.rows(),matrix.cols()),
-              m_eivalues(matrix.cols()),
-              m_schur(matrix.rows()),
-              m_isInitialized(false),
-              m_eigenvectorsOk(false),
-              m_matX(matrix.rows(),matrix.cols())
-    {
-      compute(matrix.derived(), computeEigenvectors);
-    }
-
-    /** \brief Returns the eigenvectors of given matrix.
-      *
-      * \returns  A const reference to the matrix whose columns are the eigenvectors.
-      *
-      * \pre Either the constructor
-      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
-      * function compute(const MatrixType& matrix, bool) has been called before
-      * to compute the eigendecomposition of a matrix, and
-      * \p computeEigenvectors was set to true (the default).
-      *
-      * This function returns a matrix whose columns are the eigenvectors. Column
-      * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
-      * \f$ as returned by eigenvalues().  The eigenvectors are normalized to
-      * have (Euclidean) norm equal to one. The matrix returned by this
-      * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
-      * V^{-1} \f$, if it exists.
-      *
-      * Example: \include ComplexEigenSolver_eigenvectors.cpp
-      * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
-      */
-    const EigenvectorType& eigenvectors() const
-    {
-      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
-      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
-      return m_eivec;
-    }
-
-    /** \brief Returns the eigenvalues of given matrix.
-      *
-      * \returns A const reference to the column vector containing the eigenvalues.
-      *
-      * \pre Either the constructor
-      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
-      * function compute(const MatrixType& matrix, bool) has been called before
-      * to compute the eigendecomposition of a matrix.
-      *
-      * This function returns a column vector containing the
-      * eigenvalues. Eigenvalues are repeated according to their
-      * algebraic multiplicity, so there are as many eigenvalues as
-      * rows in the matrix. The eigenvalues are not sorted in any particular
-      * order.
-      *
-      * Example: \include ComplexEigenSolver_eigenvalues.cpp
-      * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
-      */
-    const EigenvalueType& eigenvalues() const
-    {
-      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
-      return m_eivalues;
-    }
-
-    /** \brief Computes eigendecomposition of given matrix.
-      *
-      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
-      *    eigenvalues are computed; if false, only the eigenvalues are
-      *    computed.
-      * \returns    Reference to \c *this
-      *
-      * This function computes the eigenvalues of the complex matrix \p matrix.
-      * The eigenvalues() function can be used to retrieve them.  If
-      * \p computeEigenvectors is true, then the eigenvectors are also computed
-      * and can be retrieved by calling eigenvectors().
-      *
-      * The matrix is first reduced to Schur form using the
-      * ComplexSchur class. The Schur decomposition is then used to
-      * compute the eigenvalues and eigenvectors.
-      *
-      * The cost of the computation is dominated by the cost of the
-      * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
-      * is the size of the matrix.
-      *
-      * Example: \include ComplexEigenSolver_compute.cpp
-      * Output: \verbinclude ComplexEigenSolver_compute.out
-      */
-    template<typename InputType>
-    ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
-
-    /** \brief Reports whether previous computation was successful.
-      *
-      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
-      */
-    ComputationInfo info() const
-    {
-      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
-      return m_schur.info();
-    }
-
-    /** \brief Sets the maximum number of iterations allowed. */
-    ComplexEigenSolver& setMaxIterations(Index maxIters)
-    {
-      m_schur.setMaxIterations(maxIters);
-      return *this;
-    }
-
-    /** \brief Returns the maximum number of iterations. */
-    Index getMaxIterations()
-    {
-      return m_schur.getMaxIterations();
-    }
-
-  protected:
-    
-    static void check_template_parameters()
-    {
-      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
-    }
-    
-    EigenvectorType m_eivec;
-    EigenvalueType m_eivalues;
-    ComplexSchur<MatrixType> m_schur;
-    bool m_isInitialized;
-    bool m_eigenvectorsOk;
-    EigenvectorType m_matX;
-
-  private:
-    void doComputeEigenvectors(RealScalar matrixnorm);
-    void sortEigenvalues(bool computeEigenvectors);
-};
-
-
-template<typename MatrixType>
-template<typename InputType>
-ComplexEigenSolver<MatrixType>& 
-ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
-{
-  check_template_parameters();
-  
-  // this code is inspired from Jampack
-  eigen_assert(matrix.cols() == matrix.rows());
-
-  // Do a complex Schur decomposition, A = U T U^*
-  // The eigenvalues are on the diagonal of T.
-  m_schur.compute(matrix.derived(), computeEigenvectors);
-
-  if(m_schur.info() == Success)
-  {
-    m_eivalues = m_schur.matrixT().diagonal();
-    if(computeEigenvectors)
-      doComputeEigenvectors(m_schur.matrixT().norm());
-    sortEigenvalues(computeEigenvectors);
-  }
-
-  m_isInitialized = true;
-  m_eigenvectorsOk = computeEigenvectors;
-  return *this;
-}
-
-
-template<typename MatrixType>
-void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
-{
-  const Index n = m_eivalues.size();
-
-  matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)());
-
-  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
-  // The matrix X is unit triangular.
-  m_matX = EigenvectorType::Zero(n, n);
-  for(Index k=n-1 ; k>=0 ; k--)
-  {
-    m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
-    // Compute X(i,k) using the (i,k) entry of the equation X T = D X
-    for(Index i=k-1 ; i>=0 ; i--)
-    {
-      m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
-      if(k-i-1>0)
-        m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
-      ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
-      if(z==ComplexScalar(0))
-      {
-        // If the i-th and k-th eigenvalue are equal, then z equals 0.
-        // Use a small value instead, to prevent division by zero.
-        numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
-      }
-      m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
-    }
-  }
-
-  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
-  m_eivec.noalias() = m_schur.matrixU() * m_matX;
-  // .. and normalize the eigenvectors
-  for(Index k=0 ; k<n ; k++)
-  {
-    m_eivec.col(k).normalize();
-  }
-}
-
-
-template<typename MatrixType>
-void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
-{
-  const Index n =  m_eivalues.size();
-  for (Index i=0; i<n; i++)
-  {
-    Index k;
-    m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
-    if (k != 0)
-    {
-      k += i;
-      std::swap(m_eivalues[k],m_eivalues[i]);
-      if(computeEigenvectors)
-	m_eivec.col(i).swap(m_eivec.col(k));
-    }
-  }
-}
-
-} // end namespace Eigen
-
-#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
diff --git a/eigen/Eigen/src/Eigenvalues/ComplexSchur.h b/eigen/Eigen/src/Eigenvalues/ComplexSchur.h
deleted file mode 100644
index 7f38919..0000000
--- a/eigen/Eigen/src/Eigenvalues/ComplexSchur.h
+++ /dev/null
@@ -1,459 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2009 Claire Maurice
-// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_COMPLEX_SCHUR_H
-#define EIGEN_COMPLEX_SCHUR_H
-
-#include "./HessenbergDecomposition.h"
-
-namespace Eigen { 
-
-namespace internal {
-template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
-}
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class ComplexSchur
-  *
-  * \brief Performs a complex Schur decomposition of a real or complex square matrix
-  *
-  * \tparam _MatrixType the type of the matrix of which we are
-  * computing the Schur decomposition; this is expected to be an
-  * instantiation of the Matrix class template.
-  *
-  * Given a real or complex square matrix A, this class computes the
-  * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary
-  * complex matrix, and T is a complex upper triangular matrix.  The
-  * diagonal of the matrix T corresponds to the eigenvalues of the
-  * matrix A.
-  *
-  * Call the function compute() to compute the Schur decomposition of
-  * a given matrix. Alternatively, you can use the 
-  * ComplexSchur(const MatrixType&, bool) constructor which computes
-  * the Schur decomposition at construction time. Once the
-  * decomposition is computed, you can use the matrixU() and matrixT()
-  * functions to retrieve the matrices U and V in the decomposition.
-  *
-  * \note This code is inspired from Jampack
-  *
-  * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
-  */
-template<typename _MatrixType> class ComplexSchur
-{
-  public:
-    typedef _MatrixType MatrixType;
-    enum {
-      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-      Options = MatrixType::Options,
-      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-    };
-
-    /** \brief Scalar type for matrices of type \p _MatrixType. */
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-    /** \brief Complex scalar type for \p _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 std::complex<RealScalar> ComplexScalar;
-
-    /** \brief Type for the matrices in the Schur decomposition.
-      *
-      * This is a square matrix with entries of type #ComplexScalar. 
-      * The size is the same as the size of \p _MatrixType.
-      */
-    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
-
-    /** \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.
-      */
-    explicit ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
-      : m_matT(size,size),
-        m_matU(size,size),
-        m_hess(size),
-        m_isInitialized(false),
-        m_matUisUptodate(false),
-        m_maxIters(-1)
-    {}
-
-    /** \brief Constructor; computes Schur decomposition of given matrix. 
-      * 
-      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
-      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
-      *
-      * This constructor calls compute() to compute the Schur decomposition.
-      *
-      * \sa matrixT() and matrixU() for examples.
-      */
-    template<typename InputType>
-    explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
-      : m_matT(matrix.rows(),matrix.cols()),
-        m_matU(matrix.rows(),matrix.cols()),
-        m_hess(matrix.rows()),
-        m_isInitialized(false),
-        m_matUisUptodate(false),
-        m_maxIters(-1)
-    {
-      compute(matrix.derived(), computeU);
-    }
-
-    /** \brief Returns the unitary matrix in the Schur decomposition. 
-      *
-      * \returns A const reference to the matrix U.
-      *
-      * It is assumed that either the constructor
-      * ComplexSchur(const MatrixType& matrix, bool computeU) or the
-      * member function compute(const MatrixType& matrix, bool computeU)
-      * has been called before to compute the Schur decomposition of a
-      * matrix, and that \p computeU was set to true (the default
-      * value).
-      *
-      * Example: \include ComplexSchur_matrixU.cpp
-      * Output: \verbinclude ComplexSchur_matrixU.out
-      */
-    const ComplexMatrixType& matrixU() const
-    {
-      eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
-      eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
-      return m_matU;
-    }
-
-    /** \brief Returns the triangular matrix in the Schur decomposition. 
-      *
-      * \returns A const reference to the matrix T.
-      *
-      * It is assumed that either the constructor
-      * ComplexSchur(const MatrixType& matrix, bool computeU) or the
-      * member function compute(const MatrixType& matrix, bool computeU)
-      * has been called before to compute the Schur decomposition of a
-      * matrix.
-      *
-      * Note that this function returns a plain square matrix. If you want to reference
-      * only the upper triangular part, use:
-      * \code schur.matrixT().triangularView<Upper>() \endcode 
-      *
-      * Example: \include ComplexSchur_matrixT.cpp
-      * Output: \verbinclude ComplexSchur_matrixT.out
-      */
-    const ComplexMatrixType& matrixT() const
-    {
-      eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
-      return m_matT;
-    }
-
-    /** \brief Computes Schur decomposition of given matrix. 
-      * 
-      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
-      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
-
-      * \returns    Reference to \c *this
-      *
-      * The Schur decomposition is computed by first reducing the
-      * matrix to Hessenberg form using the class
-      * HessenbergDecomposition. The Hessenberg matrix is then reduced
-      * to triangular form by performing QR iterations with a single
-      * shift. The cost of computing the Schur decomposition depends
-      * on the number of iterations; as a rough guide, it may be taken
-      * on the number of iterations; as a rough guide, it may be taken
-      * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
-      * if \a computeU is false.
-      *
-      * Example: \include ComplexSchur_compute.cpp
-      * Output: \verbinclude ComplexSchur_compute.out
-      *
-      * \sa compute(const MatrixType&, bool, Index)
-      */
-    template<typename InputType>
-    ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
-    
-    /** \brief Compute Schur decomposition from a given Hessenberg matrix
-     *  \param[in] matrixH Matrix in Hessenberg form H
-     *  \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T
-     *  \param computeU Computes the matriX U of the Schur vectors
-     * \return Reference to \c *this
-     * 
-     *  This routine assumes that the matrix is already reduced in Hessenberg form matrixH
-     *  using either the class HessenbergDecomposition or another mean. 
-     *  It computes the upper quasi-triangular matrix T of the Schur decomposition of H
-     *  When computeU is true, this routine computes the matrix U such that 
-     *  A = U T U^T =  (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix
-     * 
-     * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix
-     * is not available, the user should give an identity matrix (Q.setIdentity())
-     * 
-     * \sa compute(const MatrixType&, bool)
-     */
-    template<typename HessMatrixType, typename OrthMatrixType>
-    ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU=true);
-
-    /** \brief Reports whether previous computation was successful.
-      *
-      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
-      */
-    ComputationInfo info() const
-    {
-      eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
-      return m_info;
-    }
-
-    /** \brief Sets the maximum number of iterations allowed. 
-      *
-      * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size
-      * of the matrix.
-      */
-    ComplexSchur& setMaxIterations(Index maxIters)
-    {
-      m_maxIters = maxIters;
-      return *this;
-    }
-
-    /** \brief Returns the maximum number of iterations. */
-    Index getMaxIterations()
-    {
-      return m_maxIters;
-    }
-
-    /** \brief Maximum number of iterations per row.
-      *
-      * If not otherwise specified, the maximum number of iterations is this number times the size of the
-      * matrix. It is currently set to 30.
-      */
-    static const int m_maxIterationsPerRow = 30;
-
-  protected:
-    ComplexMatrixType m_matT, m_matU;
-    HessenbergDecomposition<MatrixType> m_hess;
-    ComputationInfo m_info;
-    bool m_isInitialized;
-    bool m_matUisUptodate;
-    Index m_maxIters;
-
-  private:  
-    bool subdiagonalEntryIsNeglegible(Index i);
-    ComplexScalar computeShift(Index iu, Index iter);
-    void reduceToTriangularForm(bool computeU);
-    friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
-};
-
-/** If m_matT(i+1,i) is neglegible in floating point arithmetic
-  * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and
-  * return true, else return false. */
-template<typename MatrixType>
-inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
-{
-  RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
-  RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
-  if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
-  {
-    m_matT.coeffRef(i+1,i) = ComplexScalar(0);
-    return true;
-  }
-  return false;
-}
-
-
-/** Compute the shift in the current QR iteration. */
-template<typename MatrixType>
-typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
-{
-  using std::abs;
-  if (iter == 10 || iter == 20) 
-  {
-    // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
-    return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
-  }
-
-  // compute the shift as one of the eigenvalues of t, the 2x2
-  // diagonal block on the bottom of the active submatrix
-  Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
-  RealScalar normt = t.cwiseAbs().sum();
-  t /= normt;     // the normalization by sf is to avoid under/overflow
-
-  ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
-  ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
-  ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
-  ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
-  ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
-  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
-  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
-
-  if(numext::norm1(eival1) > numext::norm1(eival2))
-    eival2 = det / eival1;
-  else
-    eival1 = det / eival2;
-
-  // choose the eigenvalue closest to the bottom entry of the diagonal
-  if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
-    return normt * eival1;
-  else
-    return normt * eival2;
-}
-
-
-template<typename MatrixType>
-template<typename InputType>
-ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
-{
-  m_matUisUptodate = false;
-  eigen_assert(matrix.cols() == matrix.rows());
-
-  if(matrix.cols() == 1)
-  {
-    m_matT = matrix.derived().template cast<ComplexScalar>();
-    if(computeU)  m_matU = ComplexMatrixType::Identity(1,1);
-    m_info = Success;
-    m_isInitialized = true;
-    m_matUisUptodate = computeU;
-    return *this;
-  }
-
-  internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU);
-  computeFromHessenberg(m_matT, m_matU, computeU);
-  return *this;
-}
-
-template<typename MatrixType>
-template<typename HessMatrixType, typename OrthMatrixType>
-ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
-{
-  m_matT = matrixH;
-  if(computeU)
-    m_matU = matrixQ;
-  reduceToTriangularForm(computeU);
-  return *this;
-}
-namespace internal {
-
-/* Reduce given matrix to Hessenberg form */
-template<typename MatrixType, bool IsComplex>
-struct complex_schur_reduce_to_hessenberg
-{
-  // this is the implementation for the case IsComplex = true
-  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
-  {
-    _this.m_hess.compute(matrix);
-    _this.m_matT = _this.m_hess.matrixH();
-    if(computeU)  _this.m_matU = _this.m_hess.matrixQ();
-  }
-};
-
-template<typename MatrixType>
-struct complex_schur_reduce_to_hessenberg<MatrixType, false>
-{
-  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
-  {
-    typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
-
-    // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
-    _this.m_hess.compute(matrix);
-    _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
-    if(computeU)  
-    {
-      // This may cause an allocation which seems to be avoidable
-      MatrixType Q = _this.m_hess.matrixQ(); 
-      _this.m_matU = Q.template cast<ComplexScalar>();
-    }
-  }
-};
-
-} // end namespace internal
-
-// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
-template<typename MatrixType>
-void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
-{  
-  Index maxIters = m_maxIters;
-  if (maxIters == -1)
-    maxIters = m_maxIterationsPerRow * m_matT.rows();
-
-  // The matrix m_matT is divided in three parts. 
-  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
-  // Rows il,...,iu is the part we are working on (the active submatrix).
-  // Rows iu+1,...,end are already brought in triangular form.
-  Index iu = m_matT.cols() - 1;
-  Index il;
-  Index iter = 0; // number of iterations we are working on the (iu,iu) element
-  Index totalIter = 0; // number of iterations for whole matrix
-
-  while(true)
-  {
-    // find iu, the bottom row of the active submatrix
-    while(iu > 0)
-    {
-      if(!subdiagonalEntryIsNeglegible(iu-1)) break;
-      iter = 0;
-      --iu;
-    }
-
-    // if iu is zero then we are done; the whole matrix is triangularized
-    if(iu==0) break;
-
-    // if we spent too many iterations, we give up
-    iter++;
-    totalIter++;
-    if(totalIter > maxIters) break;
-
-    // find il, the top row of the active submatrix
-    il = iu-1;
-    while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
-    {
-      --il;
-    }
-
-    /* perform the QR step using Givens rotations. The first rotation
-       creates a bulge; the (il+2,il) element becomes nonzero. This
-       bulge is chased down to the bottom of the active submatrix. */
-
-    ComplexScalar shift = computeShift(iu, iter);
-    JacobiRotation<ComplexScalar> rot;
-    rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
-    m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
-    m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
-    if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
-
-    for(Index i=il+1 ; i<iu ; i++)
-    {
-      rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
-      m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
-      m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
-      m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
-      if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
-    }
-  }
-
-  if(totalIter <= maxIters)
-    m_info = Success;
-  else
-    m_info = NoConvergence;
-
-  m_isInitialized = true;
-  m_matUisUptodate = computeU;
-}
-
-} // end namespace Eigen
-
-#endif // EIGEN_COMPLEX_SCHUR_H
diff --git a/eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h b/eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h
deleted file mode 100644
index 4980a3e..0000000
--- a/eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h
+++ /dev/null
@@ -1,91 +0,0 @@
-/*
- Copyright (c) 2011, Intel Corporation. All rights reserved.
-
- Redistribution and use in source and binary forms, with or without modification,
- are permitted provided that the following conditions are met:
-
- * Redistributions of source code must retain the above copyright notice, this
-   list of conditions and the following disclaimer.
- * Redistributions in binary form must reproduce the above copyright notice,
-   this list of conditions and the following disclaimer in the documentation
-   and/or other materials provided with the distribution.
- * Neither the name of Intel Corporation nor the names of its contributors may
-   be used to endorse or promote products derived from this software without
-   specific prior written permission.
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
- WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
- ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
- ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
- ********************************************************************************
- *   Content : Eigen bindings to LAPACKe
- *    Complex Schur needed to complex unsymmetrical eigenvalues/eigenvectors.
- ********************************************************************************
-*/
-
-#ifndef EIGEN_COMPLEX_SCHUR_LAPACKE_H
-#define EIGEN_COMPLEX_SCHUR_LAPACKE_H
-
-namespace Eigen { 
-
-/** \internal Specialization for the data types supported by LAPACKe */
-
-#define EIGEN_LAPACKE_SCHUR_COMPLEX(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \
-template<> template<typename InputType> inline \
-ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
-ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, bool computeU) \
-{ \
-  typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> MatrixType; \
-  typedef MatrixType::RealScalar RealScalar; \
-  typedef std::complex<RealScalar> ComplexScalar; \
-\
-  eigen_assert(matrix.cols() == matrix.rows()); \
-\
-  m_matUisUptodate = false; \
-  if(matrix.cols() == 1) \
-  { \
-    m_matT = matrix.derived().template cast<ComplexScalar>(); \
-    if(computeU)  m_matU = ComplexMatrixType::Identity(1,1); \
-      m_info = Success; \
-      m_isInitialized = true; \
-      m_matUisUptodate = computeU; \
-      return *this; \
-  } \
-  lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), sdim, info; \
-  lapack_int matrix_order = LAPACKE_COLROW; \
-  char jobvs, sort='N'; \
-  LAPACK_##LAPACKE_PREFIX_U##_SELECT1 select = 0; \
-  jobvs = (computeU) ? 'V' : 'N'; \
-  m_matU.resize(n, n); \
-  lapack_int ldvs  = internal::convert_index<lapack_int>(m_matU.outerStride()); \
-  m_matT = matrix; \
-  lapack_int lda = internal::convert_index<lapack_int>(m_matT.outerStride()); \
-  Matrix<EIGTYPE, Dynamic, Dynamic> w; \
-  w.resize(n, 1);\
-  info = LAPACKE_##LAPACKE_PREFIX##gees( matrix_order, jobvs, sort, select, n, (LAPACKE_TYPE*)m_matT.data(), lda, &sdim, (LAPACKE_TYPE*)w.data(), (LAPACKE_TYPE*)m_matU.data(), ldvs ); \
-  if(info == 0) \
-    m_info = Success; \
-  else \
-    m_info = NoConvergence; \
-\
-  m_isInitialized = true; \
-  m_matUisUptodate = computeU; \
-  return *this; \
-\
-}
-
-EIGEN_LAPACKE_SCHUR_COMPLEX(dcomplex, lapack_complex_double, z, Z, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_LAPACKE_SCHUR_COMPLEX(scomplex, lapack_complex_float,  c, C, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_LAPACKE_SCHUR_COMPLEX(dcomplex, lapack_complex_double, z, Z, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_LAPACKE_SCHUR_COMPLEX(scomplex, lapack_complex_float,  c, C, RowMajor, LAPACK_ROW_MAJOR)
-
-} // end namespace Eigen
-
-#endif // EIGEN_COMPLEX_SCHUR_LAPACKE_H
diff --git a/eigen/Eigen/src/Eigenvalues/EigenSolver.h b/eigen/Eigen/src/Eigenvalues/EigenSolver.h
deleted file mode 100644
index f205b18..0000000
--- a/eigen/Eigen/src/Eigenvalues/EigenSolver.h
+++ /dev/null
@@ -1,622 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_EIGENSOLVER_H
-#define EIGEN_EIGENSOLVER_H
-
-#include "./RealSchur.h"
-
-namespace Eigen { 
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class EigenSolver
-  *
-  * \brief Computes eigenvalues and eigenvectors of general matrices
-  *
-  * \tparam _MatrixType the type of the matrix of which we are computing the
-  * eigendecomposition; this is expected to be an instantiation of the Matrix
-  * class template. Currently, only real matrices are supported.
-  *
-  * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
-  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$.  If
-  * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
-  * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
-  * V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
-  * have \f$ A = V D V^{-1} \f$. This is called the eigendecomposition.
-  *
-  * The eigenvalues and eigenvectors of a matrix may be complex, even when the
-  * matrix is real. However, we can choose real matrices \f$ V \f$ and \f$ D
-  * \f$ satisfying \f$ A V = V D \f$, just like the eigendecomposition, if the
-  * matrix \f$ D \f$ is not required to be diagonal, but if it is allowed to
-  * have blocks of the form
-  * \f[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f]
-  * (where \f$ u \f$ and \f$ v \f$ are real numbers) on the diagonal.  These
-  * blocks correspond to complex eigenvalue pairs \f$ u \pm iv \f$. We call
-  * this variant of the eigendecomposition the pseudo-eigendecomposition.
-  *
-  * Call the function compute() to compute the eigenvalues and eigenvectors of
-  * a given matrix. Alternatively, you can use the 
-  * EigenSolver(const MatrixType&, bool) constructor which computes the
-  * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
-  * eigenvectors are computed, they can be retrieved with the eigenvalues() and
-  * eigenvectors() functions. The pseudoEigenvalueMatrix() and
-  * pseudoEigenvectors() methods allow the construction of the
-  * pseudo-eigendecomposition.
-  *
-  * The documentation for EigenSolver(const MatrixType&, bool) contains an
-  * example of the typical use of this class.
-  *
-  * \note The implementation is adapted from
-  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
-  * Their code is based on EISPACK.
-  *
-  * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
-  */
-template<typename _MatrixType> class EigenSolver
-{
-  public:
-
-    /** \brief Synonym for the template parameter \p _MatrixType. */
-    typedef _MatrixType MatrixType;
-
-    enum {
-      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-      Options = MatrixType::Options,
-      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-    };
-
-    /** \brief Scalar type for matrices of type #MatrixType. */
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-    /** \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 std::complex<RealScalar> ComplexScalar;
-
-    /** \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 Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> 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 Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
-
-    /** \brief Default constructor.
-      *
-      * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
-      *
-      * \sa compute() for an example.
-      */
-    EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
-
-    /** \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 EigenSolver()
-      */
-    explicit EigenSolver(Index size)
-      : m_eivec(size, size),
-        m_eivalues(size),
-        m_isInitialized(false),
-        m_eigenvectorsOk(false),
-        m_realSchur(size),
-        m_matT(size, size), 
-        m_tmp(size)
-    {}
-
-    /** \brief Constructor; computes eigendecomposition of given matrix. 
-      * 
-      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
-      *    eigenvalues are computed; if false, only the eigenvalues are
-      *    computed. 
-      *
-      * This constructor calls compute() to compute the eigenvalues
-      * and eigenvectors.
-      *
-      * Example: \include EigenSolver_EigenSolver_MatrixType.cpp
-      * Output: \verbinclude EigenSolver_EigenSolver_MatrixType.out
-      *
-      * \sa compute()
-      */
-    template<typename InputType>
-    explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
-      : m_eivec(matrix.rows(), matrix.cols()),
-        m_eivalues(matrix.cols()),
-        m_isInitialized(false),
-        m_eigenvectorsOk(false),
-        m_realSchur(matrix.cols()),
-        m_matT(matrix.rows(), matrix.cols()), 
-        m_tmp(matrix.cols())
-    {
-      compute(matrix.derived(), computeEigenvectors);
-    }
-
-    /** \brief Returns the eigenvectors of given matrix. 
-      *
-      * \returns  %Matrix whose columns are the (possibly complex) eigenvectors.
-      *
-      * \pre Either the constructor 
-      * EigenSolver(const MatrixType&,bool) or the member function
-      * compute(const MatrixType&, bool) has been called before, and
-      * \p computeEigenvectors was set to true (the default).
-      *
-      * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
-      * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
-      * eigenvectors are normalized to have (Euclidean) norm equal to one. The
-      * matrix returned by this function is the matrix \f$ V \f$ in the
-      * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists.
-      *
-      * Example: \include EigenSolver_eigenvectors.cpp
-      * Output: \verbinclude EigenSolver_eigenvectors.out
-      *
-      * \sa eigenvalues(), pseudoEigenvectors()
-      */
-    EigenvectorsType eigenvectors() const;
-
-    /** \brief Returns the pseudo-eigenvectors of given matrix. 
-      *
-      * \returns  Const reference to matrix whose columns are the pseudo-eigenvectors.
-      *
-      * \pre Either the constructor 
-      * EigenSolver(const MatrixType&,bool) or the member function
-      * compute(const MatrixType&, bool) has been called before, and
-      * \p computeEigenvectors was set to true (the default).
-      *
-      * The real matrix \f$ V \f$ returned by this function and the
-      * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix()
-      * satisfy \f$ AV = VD \f$.
-      *
-      * Example: \include EigenSolver_pseudoEigenvectors.cpp
-      * Output: \verbinclude EigenSolver_pseudoEigenvectors.out
-      *
-      * \sa pseudoEigenvalueMatrix(), eigenvectors()
-      */
-    const MatrixType& pseudoEigenvectors() const
-    {
-      eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
-      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
-      return m_eivec;
-    }
-
-    /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition.
-      *
-      * \returns  A block-diagonal matrix.
-      *
-      * \pre Either the constructor 
-      * EigenSolver(const MatrixType&,bool) or the member function
-      * compute(const MatrixType&, bool) has been called before.
-      *
-      * The matrix \f$ D \f$ returned by this function is real and
-      * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2
-      * blocks of the form
-      * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$.
-      * These blocks are not sorted in any particular order.
-      * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by
-      * pseudoEigenvectors() satisfy \f$ AV = VD \f$.
-      *
-      * \sa pseudoEigenvectors() for an example, eigenvalues()
-      */
-    MatrixType pseudoEigenvalueMatrix() const;
-
-    /** \brief Returns the eigenvalues of given matrix. 
-      *
-      * \returns A const reference to the column vector containing the eigenvalues.
-      *
-      * \pre Either the constructor 
-      * EigenSolver(const MatrixType&,bool) or the member function
-      * compute(const MatrixType&, bool) has been called before.
-      *
-      * The eigenvalues are repeated according to their algebraic multiplicity,
-      * so there are as many eigenvalues as rows in the matrix. The eigenvalues 
-      * are not sorted in any particular order.
-      *
-      * Example: \include EigenSolver_eigenvalues.cpp
-      * Output: \verbinclude EigenSolver_eigenvalues.out
-      *
-      * \sa eigenvectors(), pseudoEigenvalueMatrix(),
-      *     MatrixBase::eigenvalues()
-      */
-    const EigenvalueType& eigenvalues() const
-    {
-      eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
-      return m_eivalues;
-    }
-
-    /** \brief Computes eigendecomposition of given matrix. 
-      * 
-      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
-      *    eigenvalues are computed; if false, only the eigenvalues are
-      *    computed. 
-      * \returns    Reference to \c *this
-      *
-      * This function computes the eigenvalues of the real matrix \p matrix.
-      * The eigenvalues() function can be used to retrieve them.  If 
-      * \p computeEigenvectors is true, then the eigenvectors are also computed
-      * and can be retrieved by calling eigenvectors().
-      *
-      * The matrix is first reduced to real Schur form using the RealSchur
-      * class. The Schur decomposition is then used to compute the eigenvalues
-      * and eigenvectors.
-      *
-      * The cost of the computation is dominated by the cost of the
-      * Schur decomposition, which is very approximately \f$ 25n^3 \f$
-      * (where \f$ n \f$ is the size of the matrix) if \p computeEigenvectors 
-      * is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false.
-      *
-      * This method reuses of the allocated data in the EigenSolver object.
-      *
-      * Example: \include EigenSolver_compute.cpp
-      * Output: \verbinclude EigenSolver_compute.out
-      */
-    template<typename InputType>
-    EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
-
-    /** \returns NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise. */
-    ComputationInfo info() const
-    {
-      eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
-      return m_info;
-    }
-
-    /** \brief Sets the maximum number of iterations allowed. */
-    EigenSolver& setMaxIterations(Index maxIters)
-    {
-      m_realSchur.setMaxIterations(maxIters);
-      return *this;
-    }
-
-    /** \brief Returns the maximum number of iterations. */
-    Index getMaxIterations()
-    {
-      return m_realSchur.getMaxIterations();
-    }
-
-  private:
-    void doComputeEigenvectors();
-
-  protected:
-    
-    static void check_template_parameters()
-    {
-      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
-      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
-    }
-    
-    MatrixType m_eivec;
-    EigenvalueType m_eivalues;
-    bool m_isInitialized;
-    bool m_eigenvectorsOk;
-    ComputationInfo m_info;
-    RealSchur<MatrixType> m_realSchur;
-    MatrixType m_matT;
-
-    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
-    ColumnVectorType m_tmp;
-};
-
-template<typename MatrixType>
-MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
-{
-  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
-  const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
-  Index n = m_eivalues.rows();
-  MatrixType matD = MatrixType::Zero(n,n);
-  for (Index i=0; i<n; ++i)
-  {
-    if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)), precision))
-      matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
-    else
-    {
-      matD.template block<2,2>(i,i) <<  numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
-                                       -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
-      ++i;
-    }
-  }
-  return matD;
-}
-
-template<typename MatrixType>
-typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
-{
-  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
-  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
-  const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
-  Index n = m_eivec.cols();
-  EigenvectorsType matV(n,n);
-  for (Index j=0; j<n; ++j)
-  {
-    if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n)
-    {
-      // we have a real eigen value
-      matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
-      matV.col(j).normalize();
-    }
-    else
-    {
-      // we have a pair of complex eigen values
-      for (Index i=0; i<n; ++i)
-      {
-        matV.coeffRef(i,j)   = ComplexScalar(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
-        matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
-      }
-      matV.col(j).normalize();
-      matV.col(j+1).normalize();
-      ++j;
-    }
-  }
-  return matV;
-}
-
-template<typename MatrixType>
-template<typename InputType>
-EigenSolver<MatrixType>& 
-EigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
-{
-  check_template_parameters();
-  
-  using std::sqrt;
-  using std::abs;
-  using numext::isfinite;
-  eigen_assert(matrix.cols() == matrix.rows());
-
-  // Reduce to real Schur form.
-  m_realSchur.compute(matrix.derived(), computeEigenvectors);
-  
-  m_info = m_realSchur.info();
-
-  if (m_info == Success)
-  {
-    m_matT = m_realSchur.matrixT();
-    if (computeEigenvectors)
-      m_eivec = m_realSchur.matrixU();
-  
-    // Compute eigenvalues from matT
-    m_eivalues.resize(matrix.cols());
-    Index i = 0;
-    while (i < matrix.cols()) 
-    {
-      if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
-      {
-        m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
-        if(!(isfinite)(m_eivalues.coeffRef(i)))
-        {
-          m_isInitialized = true;
-          m_eigenvectorsOk = false;
-          m_info = NumericalIssue;
-          return *this;
-        }
-        ++i;
-      }
-      else
-      {
-        Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
-        Scalar z;
-        // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
-        // without overflow
-        {
-          Scalar t0 = m_matT.coeff(i+1, i);
-          Scalar t1 = m_matT.coeff(i, i+1);
-          Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
-          t0 /= maxval;
-          t1 /= maxval;
-          Scalar p0 = p/maxval;
-          z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
-        }
-        
-        m_eivalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
-        m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
-        if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1))))
-        {
-          m_isInitialized = true;
-          m_eigenvectorsOk = false;
-          m_info = NumericalIssue;
-          return *this;
-        }
-        i += 2;
-      }
-    }
-    
-    // Compute eigenvectors.
-    if (computeEigenvectors)
-      doComputeEigenvectors();
-  }
-
-  m_isInitialized = true;
-  m_eigenvectorsOk = computeEigenvectors;
-
-  return *this;
-}
-
-
-template<typename MatrixType>
-void EigenSolver<MatrixType>::doComputeEigenvectors()
-{
-  using std::abs;
-  const Index size = m_eivec.cols();
-  const Scalar eps = NumTraits<Scalar>::epsilon();
-
-  // inefficient! this is already computed in RealSchur
-  Scalar norm(0);
-  for (Index j = 0; j < size; ++j)
-  {
-    norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
-  }
-  
-  // Backsubstitute to find vectors of upper triangular form
-  if (norm == Scalar(0))
-  {
-    return;
-  }
-
-  for (Index n = size-1; n >= 0; n--)
-  {
-    Scalar p = m_eivalues.coeff(n).real();
-    Scalar q = m_eivalues.coeff(n).imag();
-
-    // Scalar vector
-    if (q == Scalar(0))
-    {
-      Scalar lastr(0), lastw(0);
-      Index l = n;
-
-      m_matT.coeffRef(n,n) = Scalar(1);
-      for (Index i = n-1; i >= 0; i--)
-      {
-        Scalar w = m_matT.coeff(i,i) - p;
-        Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
-
-        if (m_eivalues.coeff(i).imag() < Scalar(0))
-        {
-          lastw = w;
-          lastr = r;
-        }
-        else
-        {
-          l = i;
-          if (m_eivalues.coeff(i).imag() == Scalar(0))
-          {
-            if (w != Scalar(0))
-              m_matT.coeffRef(i,n) = -r / w;
-            else
-              m_matT.coeffRef(i,n) = -r / (eps * norm);
-          }
-          else // Solve real equations
-          {
-            Scalar x = m_matT.coeff(i,i+1);
-            Scalar y = m_matT.coeff(i+1,i);
-            Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
-            Scalar t = (x * lastr - lastw * r) / denom;
-            m_matT.coeffRef(i,n) = t;
-            if (abs(x) > abs(lastw))
-              m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
-            else
-              m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
-          }
-
-          // Overflow control
-          Scalar t = abs(m_matT.coeff(i,n));
-          if ((eps * t) * t > Scalar(1))
-            m_matT.col(n).tail(size-i) /= t;
-        }
-      }
-    }
-    else if (q < Scalar(0) && n > 0) // Complex vector
-    {
-      Scalar lastra(0), lastsa(0), lastw(0);
-      Index l = n-1;
-
-      // Last vector component imaginary so matrix is triangular
-      if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
-      {
-        m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
-        m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
-      }
-      else
-      {
-        ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
-        m_matT.coeffRef(n-1,n-1) = numext::real(cc);
-        m_matT.coeffRef(n-1,n) = numext::imag(cc);
-      }
-      m_matT.coeffRef(n,n-1) = Scalar(0);
-      m_matT.coeffRef(n,n) = Scalar(1);
-      for (Index i = n-2; i >= 0; i--)
-      {
-        Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
-        Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
-        Scalar w = m_matT.coeff(i,i) - p;
-
-        if (m_eivalues.coeff(i).imag() < Scalar(0))
-        {
-          lastw = w;
-          lastra = ra;
-          lastsa = sa;
-        }
-        else
-        {
-          l = i;
-          if (m_eivalues.coeff(i).imag() == RealScalar(0))
-          {
-            ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q);
-            m_matT.coeffRef(i,n-1) = numext::real(cc);
-            m_matT.coeffRef(i,n) = numext::imag(cc);
-          }
-          else
-          {
-            // Solve complex equations
-            Scalar x = m_matT.coeff(i,i+1);
-            Scalar y = m_matT.coeff(i+1,i);
-            Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
-            Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
-            if ((vr == Scalar(0)) && (vi == Scalar(0)))
-              vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
-
-            ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi);
-            m_matT.coeffRef(i,n-1) = numext::real(cc);
-            m_matT.coeffRef(i,n) = numext::imag(cc);
-            if (abs(x) > (abs(lastw) + abs(q)))
-            {
-              m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
-              m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
-            }
-            else
-            {
-              cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q);
-              m_matT.coeffRef(i+1,n-1) = numext::real(cc);
-              m_matT.coeffRef(i+1,n) = numext::imag(cc);
-            }
-          }
-
-          // Overflow control
-          Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
-          if ((eps * t) * t > Scalar(1))
-            m_matT.block(i, n-1, size-i, 2) /= t;
-
-        }
-      }
-      
-      // We handled a pair of complex conjugate eigenvalues, so need to skip them both
-      n--;
-    }
-    else
-    {
-      eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen
-    }
-  }
-
-  // Back transformation to get eigenvectors of original matrix
-  for (Index j = size-1; j >= 0; j--)
-  {
-    m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
-    m_eivec.col(j) = m_tmp;
-  }
-}
-
-} // end namespace Eigen
-
-#endif // EIGEN_EIGENSOLVER_H
diff --git a/eigen/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h b/eigen/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
deleted file mode 100644
index 87d789b..0000000
--- a/eigen/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
+++ /dev/null
@@ -1,418 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
-// Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
-#define EIGEN_GENERALIZEDEIGENSOLVER_H
-
-#include "./RealQZ.h"
-
-namespace Eigen { 
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class GeneralizedEigenSolver
-  *
-  * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices
-  *
-  * \tparam _MatrixType the type of the matrices of which we are computing the
-  * eigen-decomposition; this is expected to be an instantiation of the Matrix
-  * class template. Currently, only real matrices are supported.
-  *
-  * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars
-  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$.  If
-  * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
-  * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
-  * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
-  * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition.
-  *
-  * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the
-  * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is
-  * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$
-  * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero,
-  * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that:
-  * \f$ \mu_i A v_i = B v_i \f$, or even \f$ \mu_i u_i^T A  = u_i^T B \f$ where \f$ u_i \f$ is
-  * called the left eigenvector.
-  *
-  * Call the function compute() to compute the generalized eigenvalues and eigenvectors of
-  * a given matrix pair. Alternatively, you can use the
-  * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the
-  * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
-  * eigenvectors are computed, they can be retrieved with the eigenvalues() and
-  * eigenvectors() functions.
-  *
-  * Here is an usage example of this class:
-  * Example: \include GeneralizedEigenSolver.cpp
-  * Output: \verbinclude GeneralizedEigenSolver.out
-  *
-  * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
-  */
-template<typename _MatrixType> class GeneralizedEigenSolver
-{
-  public:
-
-    /** \brief Synonym for the template parameter \p _MatrixType. */
-    typedef _MatrixType MatrixType;
-
-    enum {
-      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-      Options = MatrixType::Options,
-      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-    };
-
-    /** \brief Scalar type for matrices of type #MatrixType. */
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-    /** \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 std::complex<RealScalar> ComplexScalar;
-
-    /** \brief Type for vector of real scalar values eigenvalues as returned by betas().
-      *
-      * This is a column vector with entries of type #Scalar.
-      * The length of the vector is the size of #MatrixType.
-      */
-    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
-
-    /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas().
-      *
-      * This is a column vector with entries of type #ComplexScalar.
-      * The length of the vector is the size of #MatrixType.
-      */
-    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType;
-
-    /** \brief Expression type for the eigenvalues as returned by eigenvalues().
-      */
-    typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> 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 Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
-
-    /** \brief Default constructor.
-      *
-      * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
-      *
-      * \sa compute() for an example.
-      */
-    GeneralizedEigenSolver()
-      : m_eivec(),
-        m_alphas(),
-        m_betas(),
-        m_valuesOkay(false),
-        m_vectorsOkay(false),
-        m_realQZ()
-    {}
-
-    /** \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 GeneralizedEigenSolver()
-      */
-    explicit GeneralizedEigenSolver(Index size)
-      : m_eivec(size, size),
-        m_alphas(size),
-        m_betas(size),
-        m_valuesOkay(false),
-        m_vectorsOkay(false),
-        m_realQZ(size),
-        m_tmp(size)
-    {}
-
-    /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
-      * 
-      * \param[in]  A  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  B  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
-      *    eigenvalues are computed; if false, only the eigenvalues are computed.
-      *
-      * This constructor calls compute() to compute the generalized eigenvalues
-      * and eigenvectors.
-      *
-      * \sa compute()
-      */
-    GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
-      : m_eivec(A.rows(), A.cols()),
-        m_alphas(A.cols()),
-        m_betas(A.cols()),
-        m_valuesOkay(false),
-        m_vectorsOkay(false),
-        m_realQZ(A.cols()),
-        m_tmp(A.cols())
-    {
-      compute(A, B, computeEigenvectors);
-    }
-
-    /* \brief Returns the computed generalized eigenvectors.
-      *
-      * \returns  %Matrix whose columns are the (possibly complex) right eigenvectors.
-      * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
-      *
-      * \pre Either the constructor 
-      * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
-      * compute(const MatrixType&, const MatrixType& bool) has been called before, and
-      * \p computeEigenvectors was set to true (the default).
-      *
-      * \sa eigenvalues()
-      */
-    EigenvectorsType eigenvectors() const {
-      eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
-      return m_eivec;
-    }
-
-    /** \brief Returns an expression of the computed generalized eigenvalues.
-      *
-      * \returns An expression of the column vector containing the eigenvalues.
-      *
-      * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode
-      * Not that betas might contain zeros. It is therefore not recommended to use this function,
-      * but rather directly deal with the alphas and betas vectors.
-      *
-      * \pre Either the constructor 
-      * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function
-      * compute(const MatrixType&,const MatrixType&,bool) has been called before.
-      *
-      * The eigenvalues are repeated according to their algebraic multiplicity,
-      * so there are as many eigenvalues as rows in the matrix. The eigenvalues 
-      * are not sorted in any particular order.
-      *
-      * \sa alphas(), betas(), eigenvectors()
-      */
-    EigenvalueType eigenvalues() const
-    {
-      eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
-      return EigenvalueType(m_alphas,m_betas);
-    }
-
-    /** \returns A const reference to the vectors containing the alpha values
-      *
-      * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
-      *
-      * \sa betas(), eigenvalues() */
-    ComplexVectorType alphas() const
-    {
-      eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
-      return m_alphas;
-    }
-
-    /** \returns A const reference to the vectors containing the beta values
-      *
-      * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
-      *
-      * \sa alphas(), eigenvalues() */
-    VectorType betas() const
-    {
-      eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
-      return m_betas;
-    }
-
-    /** \brief Computes generalized eigendecomposition of given matrix.
-      * 
-      * \param[in]  A  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  B  Square matrix whose eigendecomposition is to be computed.
-      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
-      *    eigenvalues are computed; if false, only the eigenvalues are
-      *    computed. 
-      * \returns    Reference to \c *this
-      *
-      * This function computes the eigenvalues of the real matrix \p matrix.
-      * The eigenvalues() function can be used to retrieve them.  If 
-      * \p computeEigenvectors is true, then the eigenvectors are also computed
-      * and can be retrieved by calling eigenvectors().
-      *
-      * The matrix is first reduced to real generalized Schur form using the RealQZ
-      * class. The generalized Schur decomposition is then used to compute the eigenvalues
-      * and eigenvectors.
-      *
-      * The cost of the computation is dominated by the cost of the
-      * generalized Schur decomposition.
-      *
-      * This method reuses of the allocated data in the GeneralizedEigenSolver object.
-      */
-    GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
-
-    ComputationInfo info() const
-    {
-      eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
-      return m_realQZ.info();
-    }
-
-    /** Sets the maximal number of iterations allowed.
-    */
-    GeneralizedEigenSolver& setMaxIterations(Index maxIters)
-    {
-      m_realQZ.setMaxIterations(maxIters);
-      return *this;
-    }
-
-  protected:
-    
-    static void check_template_parameters()
-    {
-      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
-      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
-    }
-    
-    EigenvectorsType m_eivec;
-    ComplexVectorType m_alphas;
-    VectorType m_betas;
-    bool m_valuesOkay, m_vectorsOkay;
-    RealQZ<MatrixType> m_realQZ;
-    ComplexVectorType m_tmp;
-};
-
-template<typename MatrixType>
-GeneralizedEigenSolver<MatrixType>&
-GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
-{
-  check_template_parameters();
-  
-  using std::sqrt;
-  using std::abs;
-  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
-  Index size = A.cols();
-  m_valuesOkay = false;
-  m_vectorsOkay = false;
-  // Reduce to generalized real Schur form:
-  // A = Q S Z and B = Q T Z
-  m_realQZ.compute(A, B, computeEigenvectors);
-  if (m_realQZ.info() == Success)
-  {
-    // Resize storage
-    m_alphas.resize(size);
-    m_betas.resize(size);
-    if (computeEigenvectors)
-    {
-      m_eivec.resize(size,size);
-      m_tmp.resize(size);
-    }
-
-    // Aliases:
-    Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
-    ComplexVectorType &cv = m_tmp;
-    const MatrixType &mS = m_realQZ.matrixS();
-    const MatrixType &mT = m_realQZ.matrixT();
-
-    Index i = 0;
-    while (i < size)
-    {
-      if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
-      {
-        // Real eigenvalue
-        m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
-        m_betas.coeffRef(i)  = mT.diagonal().coeff(i);
-        if (computeEigenvectors)
-        {
-          v.setConstant(Scalar(0.0));
-          v.coeffRef(i) = Scalar(1.0);
-          // For singular eigenvalues do nothing more
-          if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
-          {
-            // Non-singular eigenvalue
-            const Scalar alpha = real(m_alphas.coeffRef(i));
-            const Scalar beta = m_betas.coeffRef(i);
-            for (Index j = i-1; j >= 0; j--)
-            {
-              const Index st = j+1;
-              const Index sz = i-j;
-              if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
-              {
-                // 2x2 block
-                Matrix<Scalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( v.segment(st,sz) );
-                Matrix<Scalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
-                v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
-                j--;
-              }
-              else
-              {
-                v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
-              }
-            }
-          }
-          m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v;
-          m_eivec.col(i).real().normalize();
-          m_eivec.col(i).imag().setConstant(0);
-        }
-        ++i;
-      }
-      else
-      {
-        // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
-        // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
-
-        // T =  [a 0]
-        //      [0 b]
-        RealScalar a = mT.diagonal().coeff(i),
-                   b = mT.diagonal().coeff(i+1);
-        const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
-
-        // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
-        Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
-
-        Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
-        Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
-        const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
-        m_alphas.coeffRef(i)   = conj(alpha);
-        m_alphas.coeffRef(i+1) = alpha;
-
-        if (computeEigenvectors) {
-          // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
-          cv.setZero();
-          cv.coeffRef(i+1) = Scalar(1.0);
-          // here, the "static_cast" workaound expression template issues.
-          cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
-                          / (static_cast<Scalar>(beta*mS.coeffRef(i,i))   - alpha*mT.coeffRef(i,i));
-          for (Index j = i-1; j >= 0; j--)
-          {
-            const Index st = j+1;
-            const Index sz = i+1-j;
-            if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
-            {
-              // 2x2 block
-              Matrix<ComplexScalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( cv.segment(st,sz) );
-              Matrix<ComplexScalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
-              cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
-              j--;
-            } else {
-              cv.coeffRef(j) =  cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
-                              / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
-            }
-          }
-          m_eivec.col(i+1).noalias() = (m_realQZ.matrixZ().transpose() * cv);
-          m_eivec.col(i+1).normalize();
-          m_eivec.col(i) = m_eivec.col(i+1).conjugate();
-        }
-        i += 2;
-      }
-    }
-
-    m_valuesOkay = true;
-    m_vectorsOkay = computeEigenvectors;
-  }
-  return *this;
-}
-
-} // end namespace Eigen
-
-#endif // EIGEN_GENERALIZEDEIGENSOLVER_H
diff --git a/eigen/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h b/eigen/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h
deleted file mode 100644
index 5f6bb82..0000000
--- a/eigen/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h
+++ /dev/null
@@ -1,226 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
-#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
-
-#include "./Tridiagonalization.h"
-
-namespace Eigen { 
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class GeneralizedSelfAdjointEigenSolver
-  *
-  * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
-  *
-  * \tparam _MatrixType the type of the matrix of which we are computing the
-  * eigendecomposition; this is expected to be an instantiation of the Matrix
-  * class template.
-  *
-  * This class solves the generalized eigenvalue problem
-  * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
-  * selfadjoint and the matrix \f$ B \f$ should be positive definite.
-  *
-  * Only the \b lower \b triangular \b part of the input matrix is referenced.
-  *
-  * Call the function compute() to compute the eigenvalues and eigenvectors of
-  * a given matrix. Alternatively, you can use the
-  * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
-  * constructor which computes the eigenvalues and eigenvectors at construction time.
-  * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
-  * and eigenvectors() functions.
-  *
-  * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
-  * contains an example of the typical use of this class.
-  *
-  * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
-  */
-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(). This constructor
-      * can only be used if \p _MatrixType is a fixed-size matrix; use
-      * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
-      */
-    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.
-      *
-      * This constructor is useful for dynamic-size matrices, when 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
-      */
-    explicit GeneralizedSelfAdjointEigenSolver(Index size)
-        : Base(size)
-    {}
-
-    /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
-      *
-      * \param[in]  matA  Selfadjoint matrix in matrix pencil.
-      *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  matB  Positive-definite matrix in matrix pencil.
-      *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
-      *                     Default is #ComputeEigenvectors|#Ax_lBx.
-      *
-      * This constructor calls compute(const MatrixType&, const MatrixType&, int)
-      * to compute the eigenvalues and (if requested) the eigenvectors of the
-      * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
-      * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
-      * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
-      * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
-      * \a options contains ComputeEigenvectors.
-      *
-      * In addition, the two following variants can be solved via \p options:
-      * - \c ABx_lx: \f$ ABx = \lambda x \f$
-      * - \c BAx_lx: \f$ BAx = \lambda x \f$
-      *
-      * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
-      *
-      * \sa compute(const MatrixType&, const MatrixType&, int)
-      */
-    GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
-                                      int options = ComputeEigenvectors|Ax_lBx)
-      : Base(matA.cols())
-    {
-      compute(matA, matB, options);
-    }
-
-    /** \brief Computes generalized eigendecomposition of given matrix pencil.
-      *
-      * \param[in]  matA  Selfadjoint matrix in matrix pencil.
-      *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  matB  Positive-definite matrix in matrix pencil.
-      *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
-      *                     Default is #ComputeEigenvectors|#Ax_lBx.
-      *
-      * \returns    Reference to \c *this
-      *
-      * Accoring to \p options, this function computes eigenvalues and (if requested)
-      * the eigenvectors of one of the following three generalized eigenproblems:
-      * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
-      * - \c ABx_lx: \f$ ABx = \lambda x \f$
-      * - \c BAx_lx: \f$ BAx = \lambda x \f$
-      * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
-      * matrix \f$ B \f$.
-      * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
-      *
-      * The eigenvalues() function can be used to retrieve
-      * the eigenvalues. If \p options contains ComputeEigenvectors, then the
-      * eigenvectors are also computed and can be retrieved by calling
-      * eigenvectors().
-      *
-      * The implementation uses LLT to compute the Cholesky decomposition
-      * \f$ B = LL^* \f$ and computes the classical eigendecomposition
-      * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
-      * and of \f$ L^{*} A L \f$ otherwise. This solves the
-      * generalized eigenproblem, because any solution of the generalized
-      * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
-      * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
-      * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
-      * can be made for the two other variants.
-      *
-      * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
-      *
-      * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
-      */
-    GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
-                                               int options = ComputeEigenvectors|Ax_lBx);
-
-  protected:
-
-};
-
-
-template<typename MatrixType>
-GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
-compute(const MatrixType& matA, const MatrixType& matB, int options)
-{
-  eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
-  eigen_assert((options&~(EigVecMask|GenEigMask))==0
-          && (options&EigVecMask)!=EigVecMask
-          && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
-           || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
-          && "invalid option parameter");
-
-  bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
-
-  // Compute the cholesky decomposition of matB = L L' = U'U
-  LLT<MatrixType> cholB(matB);
-
-  int type = (options&GenEigMask);
-  if(type==0)
-    type = Ax_lBx;
-
-  if(type==Ax_lBx)
-  {
-    // compute C = inv(L) A inv(L')
-    MatrixType matC = matA.template selfadjointView<Lower>();
-    cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
-    cholB.matrixU().template solveInPlace<OnTheRight>(matC);
-
-    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
-
-    // transform back the eigen vectors: evecs = inv(U) * evecs
-    if(computeEigVecs)
-      cholB.matrixU().solveInPlace(Base::m_eivec);
-  }
-  else if(type==ABx_lx)
-  {
-    // compute C = L' A L
-    MatrixType matC = matA.template selfadjointView<Lower>();
-    matC = matC * cholB.matrixL();
-    matC = cholB.matrixU() * matC;
-
-    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
-
-    // transform back the eigen vectors: evecs = inv(U) * evecs
-    if(computeEigVecs)
-      cholB.matrixU().solveInPlace(Base::m_eivec);
-  }
-  else if(type==BAx_lx)
-  {
-    // compute C = L' A L
-    MatrixType matC = matA.template selfadjointView<Lower>();
-    matC = matC * cholB.matrixL();
-    matC = cholB.matrixU() * matC;
-
-    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
-
-    // transform back the eigen vectors: evecs = L * evecs
-    if(computeEigVecs)
-      Base::m_eivec = cholB.matrixL() * Base::m_eivec;
-  }
-
-  return *this;
-}
-
-} // end namespace Eigen
-
-#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
diff --git a/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h
deleted file mode 100644
index f647f69..0000000
--- a/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ /dev/null
@@ -1,374 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_HESSENBERGDECOMPOSITION_H
-#define EIGEN_HESSENBERGDECOMPOSITION_H
-
-namespace Eigen { 
-
-namespace internal {
-  
-template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
-template<typename MatrixType>
-struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
-{
-  typedef MatrixType ReturnType;
-};
-
-}
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class HessenbergDecomposition
-  *
-  * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
-  *
-  * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
-  *
-  * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
-  * the real case, the Hessenberg decomposition consists of an orthogonal
-  * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
-  * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
-  * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
-  * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
-  * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
-  * \f$ Q^{-1} = Q^* \f$).
-  *
-  * Call the function compute() to compute the Hessenberg decomposition of a
-  * given matrix. Alternatively, you can use the
-  * HessenbergDecomposition(const MatrixType&) constructor which computes the
-  * Hessenberg decomposition at construction time. Once the decomposition is
-  * computed, you can use the matrixH() and matrixQ() functions to construct
-  * the matrices H and Q in the decomposition.
-  *
-  * The documentation for matrixH() contains an example of the typical use of
-  * this class.
-  *
-  * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
-  */
-template<typename _MatrixType> class HessenbergDecomposition
-{
-  public:
-
-    /** \brief Synonym for the template parameter \p _MatrixType. */
-    typedef _MatrixType MatrixType;
-
-    enum {
-      Size = MatrixType::RowsAtCompileTime,
-      SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
-      Options = MatrixType::Options,
-      MaxSize = MatrixType::MaxRowsAtCompileTime,
-      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
-    };
-
-    /** \brief Scalar type for matrices of type #MatrixType. */
-    typedef typename MatrixType::Scalar Scalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-    /** \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, if it is a fixed-side
-      * type.
-      */
-    typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
-
-    /** \brief Return type of matrixQ() */
-    typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
-    
-    typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
-
-    /** \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.
-      */
-    explicit HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
-      : m_matrix(size,size),
-        m_temp(size),
-        m_isInitialized(false)
-    {
-      if(size>1)
-        m_hCoeffs.resize(size-1);
-    }
-
-    /** \brief Constructor; computes Hessenberg decomposition of given matrix.
-      *
-      * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
-      *
-      * This constructor calls compute() to compute the Hessenberg
-      * decomposition.
-      *
-      * \sa matrixH() for an example.
-      */
-    template<typename InputType>
-    explicit HessenbergDecomposition(const EigenBase<InputType>& matrix)
-      : m_matrix(matrix.derived()),
-        m_temp(matrix.rows()),
-        m_isInitialized(false)
-    {
-      if(matrix.rows()<2)
-      {
-        m_isInitialized = true;
-        return;
-      }
-      m_hCoeffs.resize(matrix.rows()-1,1);
-      _compute(m_matrix, m_hCoeffs, m_temp);
-      m_isInitialized = true;
-    }
-
-    /** \brief Computes Hessenberg decomposition of given matrix.
-      *
-      * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
-      * \returns    Reference to \c *this
-      *
-      * The Hessenberg decomposition is computed by bringing the columns of the
-      * matrix successively in the required form using Householder reflections
-      * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
-      * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
-      * denotes the size of the given matrix.
-      *
-      * This method reuses of the allocated data in the HessenbergDecomposition
-      * object.
-      *
-      * Example: \include HessenbergDecomposition_compute.cpp
-      * Output: \verbinclude HessenbergDecomposition_compute.out
-      */
-    template<typename InputType>
-    HessenbergDecomposition& compute(const EigenBase<InputType>& matrix)
-    {
-      m_matrix = matrix.derived();
-      if(matrix.rows()<2)
-      {
-        m_isInitialized = true;
-        return *this;
-      }
-      m_hCoeffs.resize(matrix.rows()-1,1);
-      _compute(m_matrix, m_hCoeffs, m_temp);
-      m_isInitialized = true;
-      return *this;
-    }
-
-    /** \brief Returns the Householder coefficients.
-      *
-      * \returns a const reference to the vector of Householder coefficients
-      *
-      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
-      * or the member function compute(const MatrixType&) has been called
-      * before to compute the Hessenberg decomposition of a matrix.
-      *
-      * The Householder coefficients allow the reconstruction of the matrix
-      * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
-      *
-      * \sa packedMatrix(), \ref Householder_Module "Householder module"
-      */
-    const CoeffVectorType& householderCoefficients() const
-    {
-      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
-      return m_hCoeffs;
-    }
-
-    /** \brief Returns the internal representation of the decomposition
-      *
-      *	\returns a const reference to a matrix with the internal representation
-      *	         of the decomposition.
-      *
-      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
-      * or the member function compute(const MatrixType&) has been called
-      * before to compute the Hessenberg decomposition of a matrix.
-      *
-      * The returned matrix contains the following information:
-      *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
-      *  - the rest of the lower part contains the Householder vectors that, combined with
-      *    Householder coefficients returned by householderCoefficients(),
-      *    allows to reconstruct the matrix Q as
-      *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
-      *    Here, the matrices \f$ H_i \f$ are the Householder transformations
-      *       \f$ H_i = (I - h_i v_i v_i^T) \f$
-      *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
-      *    \f$ v_i \f$ is the Householder vector defined by
-      *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
-      *    with M the matrix returned by this function.
-      *
-      * See LAPACK for further details on this packed storage.
-      *
-      * Example: \include HessenbergDecomposition_packedMatrix.cpp
-      * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
-      *
-      * \sa householderCoefficients()
-      */
-    const MatrixType& packedMatrix() const
-    {
-      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
-      return m_matrix;
-    }
-
-    /** \brief Reconstructs the orthogonal matrix Q in the decomposition
-      *
-      * \returns object representing the matrix Q
-      *
-      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
-      * or the member function compute(const MatrixType&) has been called
-      * before to compute the Hessenberg decomposition of a matrix.
-      *
-      * This function returns a light-weight object of template class
-      * HouseholderSequence. You can either apply it directly to a matrix or
-      * you can convert it to a matrix of type #MatrixType.
-      *
-      * \sa matrixH() for an example, class HouseholderSequence
-      */
-    HouseholderSequenceType matrixQ() const
-    {
-      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
-      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
-             .setLength(m_matrix.rows() - 1)
-             .setShift(1);
-    }
-
-    /** \brief Constructs the Hessenberg matrix H in the decomposition
-      *
-      * \returns expression object representing the matrix H
-      *
-      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
-      * or the member function compute(const MatrixType&) has been called
-      * before to compute the Hessenberg decomposition of a matrix.
-      *
-      * The object returned by this function constructs the Hessenberg matrix H
-      * when it is assigned to a matrix or otherwise evaluated. The matrix H is
-      * constructed from the packed matrix as returned by packedMatrix(): The
-      * upper part (including the subdiagonal) of the packed matrix contains
-      * the matrix H. It may sometimes be better to directly use the packed
-      * matrix instead of constructing the matrix H.
-      *
-      * Example: \include HessenbergDecomposition_matrixH.cpp
-      * Output: \verbinclude HessenbergDecomposition_matrixH.out
-      *
-      * \sa matrixQ(), packedMatrix()
-      */
-    MatrixHReturnType matrixH() const
-    {
-      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
-      return MatrixHReturnType(*this);
-    }
-
-  private:
-
-    typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
-    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
-
-  protected:
-    MatrixType m_matrix;
-    CoeffVectorType m_hCoeffs;
-    VectorType m_temp;
-    bool m_isInitialized;
-};
-
-/** \internal
-  * Performs a tridiagonal decomposition of \a matA in place.
-  *
-  * \param matA the input selfadjoint matrix
-  * \param hCoeffs returned Householder coefficients
-  *
-  * The result is written in the lower triangular part of \a matA.
-  *
-  * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
-  *
-  * \sa packedMatrix()
-  */
-template<typename MatrixType>
-void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
-{
-  eigen_assert(matA.rows()==matA.cols());
-  Index n = matA.rows();
-  temp.resize(n);
-  for (Index i = 0; i<n-1; ++i)
-  {
-    // let's consider the vector v = i-th column starting at position i+1
-    Index remainingSize = n-i-1;
-    RealScalar beta;
-    Scalar h;
-    matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
-    matA.col(i).coeffRef(i+1) = beta;
-    hCoeffs.coeffRef(i) = h;
-
-    // Apply similarity transformation to remaining columns,
-    // i.e., compute A = H A H'
-
-    // A = H A
-    matA.bottomRightCorner(remainingSize, remainingSize)
-        .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
-
-    // A = A H'
-    matA.rightCols(remainingSize)
-        .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0));
-  }
-}
-
-namespace internal {
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \brief Expression type for return value of HessenbergDecomposition::matrixH()
-  *
-  * \tparam MatrixType type of matrix in the Hessenberg decomposition
-  *
-  * Objects of this type represent the Hessenberg matrix in the Hessenberg
-  * decomposition of some matrix. The object holds a reference to the
-  * HessenbergDecomposition class until the it is assigned or evaluated for
-  * some other reason (the reference should remain valid during the life time
-  * of this object). This class is the return type of
-  * HessenbergDecomposition::matrixH(); there is probably no other use for this
-  * class.
-  */
-template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
-: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
-{
-  public:
-    /** \brief Constructor.
-      *
-      * \param[in] hess  Hessenberg decomposition
-      */
-    HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
-
-    /** \brief Hessenberg matrix in decomposition.
-      *
-      * \param[out] result  Hessenberg matrix in decomposition \p hess which
-      *                     was passed to the constructor
-      */
-    template <typename ResultType>
-    inline void evalTo(ResultType& result) const
-    {
-      result = m_hess.packedMatrix();
-      Index n = result.rows();
-      if (n>2)
-        result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
-    }
-
-    Index rows() const { return m_hess.packedMatrix().rows(); }
-    Index cols() const { return m_hess.packedMatrix().cols(); }
-
-  protected:
-    const HessenbergDecomposition<MatrixType>& m_hess;
-};
-
-} // end namespace internal
-
-} // end namespace Eigen
-
-#endif // EIGEN_HESSENBERGDECOMPOSITION_H
diff --git a/eigen/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h b/eigen/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
deleted file mode 100644
index e4e4260..0000000
--- a/eigen/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
+++ /dev/null
@@ -1,158 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_MATRIXBASEEIGENVALUES_H
-#define EIGEN_MATRIXBASEEIGENVALUES_H
-
-namespace Eigen { 
-
-namespace internal {
-
-template<typename Derived, bool IsComplex>
-struct eigenvalues_selector
-{
-  // this is the implementation for the case IsComplex = true
-  static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
-  run(const MatrixBase<Derived>& m)
-  {
-    typedef typename Derived::PlainObject PlainObject;
-    PlainObject m_eval(m);
-    return ComplexEigenSolver<PlainObject>(m_eval, false).eigenvalues();
-  }
-};
-
-template<typename Derived>
-struct eigenvalues_selector<Derived, false>
-{
-  static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
-  run(const MatrixBase<Derived>& m)
-  {
-    typedef typename Derived::PlainObject PlainObject;
-    PlainObject m_eval(m);
-    return EigenSolver<PlainObject>(m_eval, false).eigenvalues();
-  }
-};
-
-} // end namespace internal
-
-/** \brief Computes the eigenvalues of a matrix 
-  * \returns Column vector containing the eigenvalues.
-  *
-  * \eigenvalues_module
-  * This function computes the eigenvalues with the help of the EigenSolver
-  * class (for real matrices) or the ComplexEigenSolver class (for complex
-  * matrices). 
-  *
-  * The eigenvalues are repeated according to their algebraic multiplicity,
-  * so there are as many eigenvalues as rows in the matrix.
-  *
-  * The SelfAdjointView class provides a better algorithm for selfadjoint
-  * matrices.
-  *
-  * Example: \include MatrixBase_eigenvalues.cpp
-  * Output: \verbinclude MatrixBase_eigenvalues.out
-  *
-  * \sa EigenSolver::eigenvalues(), ComplexEigenSolver::eigenvalues(),
-  *     SelfAdjointView::eigenvalues()
-  */
-template<typename Derived>
-inline typename MatrixBase<Derived>::EigenvaluesReturnType
-MatrixBase<Derived>::eigenvalues() const
-{
-  return internal::eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
-}
-
-/** \brief Computes the eigenvalues of a matrix
-  * \returns Column vector containing the eigenvalues.
-  *
-  * \eigenvalues_module
-  * This function computes the eigenvalues with the help of the
-  * SelfAdjointEigenSolver class.  The eigenvalues are repeated according to
-  * their algebraic multiplicity, so there are as many eigenvalues as rows in
-  * the matrix.
-  *
-  * Example: \include SelfAdjointView_eigenvalues.cpp
-  * Output: \verbinclude SelfAdjointView_eigenvalues.out
-  *
-  * \sa SelfAdjointEigenSolver::eigenvalues(), MatrixBase::eigenvalues()
-  */
-template<typename MatrixType, unsigned int UpLo> 
-inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType
-SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
-{
-  PlainObject thisAsMatrix(*this);
-  return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix, false).eigenvalues();
-}
-
-
-
-/** \brief Computes the L2 operator norm
-  * \returns Operator norm of the matrix.
-  *
-  * \eigenvalues_module
-  * This function computes the L2 operator norm of a matrix, which is also
-  * known as the spectral norm. The norm of a matrix \f$ A \f$ is defined to be
-  * \f[ \|A\|_2 = \max_x \frac{\|Ax\|_2}{\|x\|_2} \f]
-  * where the maximum is over all vectors and the norm on the right is the
-  * Euclidean vector norm. The norm equals the largest singular value, which is
-  * the square root of the largest eigenvalue of the positive semi-definite
-  * matrix \f$ A^*A \f$.
-  *
-  * The current implementation uses the eigenvalues of \f$ A^*A \f$, as computed
-  * by SelfAdjointView::eigenvalues(), to compute the operator norm of a
-  * matrix.  The SelfAdjointView class provides a better algorithm for
-  * selfadjoint matrices.
-  *
-  * Example: \include MatrixBase_operatorNorm.cpp
-  * Output: \verbinclude MatrixBase_operatorNorm.out
-  *
-  * \sa SelfAdjointView::eigenvalues(), SelfAdjointView::operatorNorm()
-  */
-template<typename Derived>
-inline typename MatrixBase<Derived>::RealScalar
-MatrixBase<Derived>::operatorNorm() const
-{
-  using std::sqrt;
-  typename Derived::PlainObject m_eval(derived());
-  // FIXME if it is really guaranteed that the eigenvalues are already sorted,
-  // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
-  return sqrt((m_eval*m_eval.adjoint())
-                 .eval()
-		 .template selfadjointView<Lower>()
-		 .eigenvalues()
-		 .maxCoeff()
-		 );
-}
-
-/** \brief Computes the L2 operator norm
-  * \returns Operator norm of the matrix.
-  *
-  * \eigenvalues_module
-  * This function computes the L2 operator norm of a self-adjoint matrix. For a
-  * self-adjoint matrix, the operator norm is the largest eigenvalue.
-  *
-  * The current implementation uses the eigenvalues of the matrix, as computed
-  * by eigenvalues(), to compute the operator norm of the matrix.
-  *
-  * Example: \include SelfAdjointView_operatorNorm.cpp
-  * Output: \verbinclude SelfAdjointView_operatorNorm.out
-  *
-  * \sa eigenvalues(), MatrixBase::operatorNorm()
-  */
-template<typename MatrixType, unsigned int UpLo>
-inline typename SelfAdjointView<MatrixType, UpLo>::RealScalar
-SelfAdjointView<MatrixType, UpLo>::operatorNorm() const
-{
-  return eigenvalues().cwiseAbs().maxCoeff();
-}
-
-} // end namespace Eigen
-
-#endif
diff --git a/eigen/Eigen/src/Eigenvalues/RealQZ.h b/eigen/Eigen/src/Eigenvalues/RealQZ.h
deleted file mode 100644
index b3a910d..0000000
--- a/eigen/Eigen/src/Eigenvalues/RealQZ.h
+++ /dev/null
@@ -1,654 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_REAL_QZ_H
-#define EIGEN_REAL_QZ_H
-
-namespace Eigen {
-
-  /** \eigenvalues_module \ingroup Eigenvalues_Module
-   *
-   *
-   * \class RealQZ
-   *
-   * \brief Performs a real QZ decomposition of a pair of square matrices
-   *
-   * \tparam _MatrixType the type of the matrix of which we are computing the
-   * real QZ decomposition; this is expected to be an instantiation of the
-   * Matrix class template.
-   *
-   * Given a real square matrices A and B, this class computes the real QZ
-   * decomposition: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ where Q and Z are
-   * real orthogonal matrixes, T is upper-triangular matrix, and S is upper
-   * quasi-triangular matrix. An orthogonal matrix is a matrix whose
-   * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular
-   * matrix is a block-triangular matrix whose diagonal consists of 1-by-1
-   * blocks and 2-by-2 blocks where further reduction is impossible due to
-   * complex eigenvalues. 
-   *
-   * The eigenvalues of the pencil \f$ A - z B \f$ can be obtained from
-   * 1x1 and 2x2 blocks on the diagonals of S and T.
-   *
-   * Call the function compute() to compute the real QZ decomposition of a
-   * given pair of matrices. Alternatively, you can use the 
-   * RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ)
-   * constructor which computes the real QZ decomposition at construction
-   * time. Once the decomposition is computed, you can use the matrixS(),
-   * matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices
-   * S, T, Q and Z in the decomposition. If computeQZ==false, some time
-   * is saved by not computing matrices Q and Z.
-   *
-   * Example: \include RealQZ_compute.cpp
-   * Output: \include RealQZ_compute.out
-   *
-   * \note The implementation is based on the algorithm in "Matrix Computations"
-   * by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for
-   * generalized eigenvalue problems" by C.B.Moler and G.W.Stewart.
-   *
-   * \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver
-   */
-
-  template<typename _MatrixType> class RealQZ
-  {
-    public:
-      typedef _MatrixType MatrixType;
-      enum {
-        RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-        ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-        Options = MatrixType::Options,
-        MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-        MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-      };
-      typedef typename MatrixType::Scalar Scalar;
-      typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
-      typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-      typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
-      typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
-
-      /** \brief Default constructor.
-       *
-       * \param [in] size  Positive integer, size of the matrix whose QZ 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.
-       */
-      explicit RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) :
-        m_S(size, size),
-        m_T(size, size),
-        m_Q(size, size),
-        m_Z(size, size),
-        m_workspace(size*2),
-        m_maxIters(400),
-        m_isInitialized(false)
-        { }
-
-      /** \brief Constructor; computes real QZ decomposition of given matrices
-       * 
-       * \param[in]  A          Matrix A.
-       * \param[in]  B          Matrix B.
-       * \param[in]  computeQZ  If false, A and Z are not computed.
-       *
-       * This constructor calls compute() to compute the QZ decomposition.
-       */
-      RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true) :
-        m_S(A.rows(),A.cols()),
-        m_T(A.rows(),A.cols()),
-        m_Q(A.rows(),A.cols()),
-        m_Z(A.rows(),A.cols()),
-        m_workspace(A.rows()*2),
-        m_maxIters(400),
-        m_isInitialized(false) {
-          compute(A, B, computeQZ);
-        }
-
-      /** \brief Returns matrix Q in the QZ decomposition. 
-       *
-       * \returns A const reference to the matrix Q.
-       */
-      const MatrixType& matrixQ() const {
-        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
-        eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
-        return m_Q;
-      }
-
-      /** \brief Returns matrix Z in the QZ decomposition. 
-       *
-       * \returns A const reference to the matrix Z.
-       */
-      const MatrixType& matrixZ() const {
-        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
-        eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
-        return m_Z;
-      }
-
-      /** \brief Returns matrix S in the QZ decomposition. 
-       *
-       * \returns A const reference to the matrix S.
-       */
-      const MatrixType& matrixS() const {
-        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
-        return m_S;
-      }
-
-      /** \brief Returns matrix S in the QZ decomposition. 
-       *
-       * \returns A const reference to the matrix S.
-       */
-      const MatrixType& matrixT() const {
-        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
-        return m_T;
-      }
-
-      /** \brief Computes QZ decomposition of given matrix. 
-       * 
-       * \param[in]  A          Matrix A.
-       * \param[in]  B          Matrix B.
-       * \param[in]  computeQZ  If false, A and Z are not computed.
-       * \returns    Reference to \c *this
-       */
-      RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true);
-
-      /** \brief Reports whether previous computation was successful.
-       *
-       * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
-       */
-      ComputationInfo info() const
-      {
-        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
-        return m_info;
-      }
-
-      /** \brief Returns number of performed QR-like iterations.
-      */
-      Index iterations() const
-      {
-        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
-        return m_global_iter;
-      }
-
-      /** Sets the maximal number of iterations allowed to converge to one eigenvalue
-       * or decouple the problem.
-      */
-      RealQZ& setMaxIterations(Index maxIters)
-      {
-        m_maxIters = maxIters;
-        return *this;
-      }
-
-    private:
-
-      MatrixType m_S, m_T, m_Q, m_Z;
-      Matrix<Scalar,Dynamic,1> m_workspace;
-      ComputationInfo m_info;
-      Index m_maxIters;
-      bool m_isInitialized;
-      bool m_computeQZ;
-      Scalar m_normOfT, m_normOfS;
-      Index m_global_iter;
-
-      typedef Matrix<Scalar,3,1> Vector3s;
-      typedef Matrix<Scalar,2,1> Vector2s;
-      typedef Matrix<Scalar,2,2> Matrix2s;
-      typedef JacobiRotation<Scalar> JRs;
-
-      void hessenbergTriangular();
-      void computeNorms();
-      Index findSmallSubdiagEntry(Index iu);
-      Index findSmallDiagEntry(Index f, Index l);
-      void splitOffTwoRows(Index i);
-      void pushDownZero(Index z, Index f, Index l);
-      void step(Index f, Index l, Index iter);
-
-  }; // RealQZ
-
-  /** \internal Reduces S and T to upper Hessenberg - triangular form */
-  template<typename MatrixType>
-    void RealQZ<MatrixType>::hessenbergTriangular()
-    {
-
-      const Index dim = m_S.cols();
-
-      // perform QR decomposition of T, overwrite T with R, save Q
-      HouseholderQR<MatrixType> qrT(m_T);
-      m_T = qrT.matrixQR();
-      m_T.template triangularView<StrictlyLower>().setZero();
-      m_Q = qrT.householderQ();
-      // overwrite S with Q* S
-      m_S.applyOnTheLeft(m_Q.adjoint());
-      // init Z as Identity
-      if (m_computeQZ)
-        m_Z = MatrixType::Identity(dim,dim);
-      // reduce S to upper Hessenberg with Givens rotations
-      for (Index j=0; j<=dim-3; j++) {
-        for (Index i=dim-1; i>=j+2; i--) {
-          JRs G;
-          // kill S(i,j)
-          if(m_S.coeff(i,j) != 0)
-          {
-            G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j));
-            m_S.coeffRef(i,j) = Scalar(0.0);
-            m_S.rightCols(dim-j-1).applyOnTheLeft(i-1,i,G.adjoint());
-            m_T.rightCols(dim-i+1).applyOnTheLeft(i-1,i,G.adjoint());
-            // update Q
-            if (m_computeQZ)
-              m_Q.applyOnTheRight(i-1,i,G);
-          }
-          // kill T(i,i-1)
-          if(m_T.coeff(i,i-1)!=Scalar(0))
-          {
-            G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i));
-            m_T.coeffRef(i,i-1) = Scalar(0.0);
-            m_S.applyOnTheRight(i,i-1,G);
-            m_T.topRows(i).applyOnTheRight(i,i-1,G);
-            // update Z
-            if (m_computeQZ)
-              m_Z.applyOnTheLeft(i,i-1,G.adjoint());
-          }
-        }
-      }
-    }
-
-  /** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */
-  template<typename MatrixType>
-    inline void RealQZ<MatrixType>::computeNorms()
-    {
-      const Index size = m_S.cols();
-      m_normOfS = Scalar(0.0);
-      m_normOfT = Scalar(0.0);
-      for (Index j = 0; j < size; ++j)
-      {
-        m_normOfS += m_S.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
-        m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum();
-      }
-    }
-
-
-  /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */
-  template<typename MatrixType>
-    inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu)
-    {
-      using std::abs;
-      Index res = iu;
-      while (res > 0)
-      {
-        Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res));
-        if (s == Scalar(0.0))
-          s = m_normOfS;
-        if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
-          break;
-        res--;
-      }
-      return res;
-    }
-
-  /** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1)  */
-  template<typename MatrixType>
-    inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l)
-    {
-      using std::abs;
-      Index res = l;
-      while (res >= f) {
-        if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
-          break;
-        res--;
-      }
-      return res;
-    }
-
-  /** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */
-  template<typename MatrixType>
-    inline void RealQZ<MatrixType>::splitOffTwoRows(Index i)
-    {
-      using std::abs;
-      using std::sqrt;
-      const Index dim=m_S.cols();
-      if (abs(m_S.coeff(i+1,i))==Scalar(0))
-        return;
-      Index j = findSmallDiagEntry(i,i+1);
-      if (j==i-1)
-      {
-        // block of (S T^{-1})
-        Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>().
-          template solve<OnTheRight>(m_S.template block<2,2>(i,i));
-        Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1));
-        Scalar q = p*p + STi(1,0)*STi(0,1);
-        if (q>=0) {
-          Scalar z = sqrt(q);
-          // one QR-like iteration for ABi - lambda I
-          // is enough - when we know exact eigenvalue in advance,
-          // convergence is immediate
-          JRs G;
-          if (p>=0)
-            G.makeGivens(p + z, STi(1,0));
-          else
-            G.makeGivens(p - z, STi(1,0));
-          m_S.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint());
-          m_T.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint());
-          // update Q
-          if (m_computeQZ)
-            m_Q.applyOnTheRight(i,i+1,G);
-
-          G.makeGivens(m_T.coeff(i+1,i+1), m_T.coeff(i+1,i));
-          m_S.topRows(i+2).applyOnTheRight(i+1,i,G);
-          m_T.topRows(i+2).applyOnTheRight(i+1,i,G);
-          // update Z
-          if (m_computeQZ)
-            m_Z.applyOnTheLeft(i+1,i,G.adjoint());
-
-          m_S.coeffRef(i+1,i) = Scalar(0.0);
-          m_T.coeffRef(i+1,i) = Scalar(0.0);
-        }
-      }
-      else
-      {
-        pushDownZero(j,i,i+1);
-      }
-    }
-
-  /** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */
-  template<typename MatrixType>
-    inline void RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l)
-    {
-      JRs G;
-      const Index dim = m_S.cols();
-      for (Index zz=z; zz<l; zz++)
-      {
-        // push 0 down
-        Index firstColS = zz>f ? (zz-1) : zz;
-        G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1));
-        m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,G.adjoint());
-        m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,G.adjoint());
-        m_T.coeffRef(zz+1,zz+1) = Scalar(0.0);
-        // update Q
-        if (m_computeQZ)
-          m_Q.applyOnTheRight(zz,zz+1,G);
-        // kill S(zz+1, zz-1)
-        if (zz>f)
-        {
-          G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1));
-          m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,G);
-          m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,G);
-          m_S.coeffRef(zz+1,zz-1) = Scalar(0.0);
-          // update Z
-          if (m_computeQZ)
-            m_Z.applyOnTheLeft(zz,zz-1,G.adjoint());
-        }
-      }
-      // finally kill S(l,l-1)
-      G.makeGivens(m_S.coeff(l,l), m_S.coeff(l,l-1));
-      m_S.applyOnTheRight(l,l-1,G);
-      m_T.applyOnTheRight(l,l-1,G);
-      m_S.coeffRef(l,l-1)=Scalar(0.0);
-      // update Z
-      if (m_computeQZ)
-        m_Z.applyOnTheLeft(l,l-1,G.adjoint());
-    }
-
-  /** \internal QR-like iterative step for block f..l */
-  template<typename MatrixType>
-    inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter)
-    {
-      using std::abs;
-      const Index dim = m_S.cols();
-
-      // x, y, z
-      Scalar x, y, z;
-      if (iter==10)
-      {
-        // Wilkinson ad hoc shift
-        const Scalar
-          a11=m_S.coeff(f+0,f+0), a12=m_S.coeff(f+0,f+1),
-          a21=m_S.coeff(f+1,f+0), a22=m_S.coeff(f+1,f+1), a32=m_S.coeff(f+2,f+1),
-          b12=m_T.coeff(f+0,f+1),
-          b11i=Scalar(1.0)/m_T.coeff(f+0,f+0),
-          b22i=Scalar(1.0)/m_T.coeff(f+1,f+1),
-          a87=m_S.coeff(l-1,l-2),
-          a98=m_S.coeff(l-0,l-1),
-          b77i=Scalar(1.0)/m_T.coeff(l-2,l-2),
-          b88i=Scalar(1.0)/m_T.coeff(l-1,l-1);
-        Scalar ss = abs(a87*b77i) + abs(a98*b88i),
-               lpl = Scalar(1.5)*ss,
-               ll = ss*ss;
-        x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i
-          - a11*a21*b12*b11i*b11i*b22i;
-        y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i 
-          - a21*a21*b12*b11i*b11i*b22i;
-        z = a21*a32*b11i*b22i;
-      }
-      else if (iter==16)
-      {
-        // another exceptional shift
-        x = m_S.coeff(f,f)/m_T.coeff(f,f)-m_S.coeff(l,l)/m_T.coeff(l,l) + m_S.coeff(l,l-1)*m_T.coeff(l-1,l) /
-          (m_T.coeff(l-1,l-1)*m_T.coeff(l,l));
-        y = m_S.coeff(f+1,f)/m_T.coeff(f,f);
-        z = 0;
-      }
-      else if (iter>23 && !(iter%8))
-      {
-        // extremely exceptional shift
-        x = internal::random<Scalar>(-1.0,1.0);
-        y = internal::random<Scalar>(-1.0,1.0);
-        z = internal::random<Scalar>(-1.0,1.0);
-      }
-      else
-      {
-        // Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1
-        // where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where
-        // U and V are 2x2 bottom right sub matrices of A and B. Thus:
-        //  = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1)
-        //  = AB^-1AB^-1 + det(M) - tr(M)(AB^-1)
-        // Since we are only interested in having x, y, z with a correct ratio, we have:
-        const Scalar
-          a11 = m_S.coeff(f,f),     a12 = m_S.coeff(f,f+1),
-          a21 = m_S.coeff(f+1,f),   a22 = m_S.coeff(f+1,f+1),
-                                    a32 = m_S.coeff(f+2,f+1),
-
-          a88 = m_S.coeff(l-1,l-1), a89 = m_S.coeff(l-1,l),
-          a98 = m_S.coeff(l,l-1),   a99 = m_S.coeff(l,l),
-
-          b11 = m_T.coeff(f,f),     b12 = m_T.coeff(f,f+1),
-                                    b22 = m_T.coeff(f+1,f+1),
-
-          b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l),
-                                    b99 = m_T.coeff(l,l);
-
-        x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21)
-          + a12/b22 - (a11/b11)*(b12/b22);
-        y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99);
-        z = a32/b22;
-      }
-
-      JRs G;
-
-      for (Index k=f; k<=l-2; k++)
-      {
-        // variables for Householder reflections
-        Vector2s essential2;
-        Scalar tau, beta;
-
-        Vector3s hr(x,y,z);
-
-        // Q_k to annihilate S(k+1,k-1) and S(k+2,k-1)
-        hr.makeHouseholderInPlace(tau, beta);
-        essential2 = hr.template bottomRows<2>();
-        Index fc=(std::max)(k-1,Index(0));  // first col to update
-        m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data());
-        m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data());
-        if (m_computeQZ)
-          m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data());
-        if (k>f)
-          m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0);
-
-        // Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k)
-        hr << m_T.coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1);
-        hr.makeHouseholderInPlace(tau, beta);
-        essential2 = hr.template bottomRows<2>();
-        {
-          Index lr = (std::min)(k+4,dim); // last row to update
-          Map<Matrix<Scalar,Dynamic,1> > tmp(m_workspace.data(),lr);
-          // S
-          tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2;
-          tmp += m_S.col(k+2).head(lr);
-          m_S.col(k+2).head(lr) -= tau*tmp;
-          m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
-          // T
-          tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2;
-          tmp += m_T.col(k+2).head(lr);
-          m_T.col(k+2).head(lr) -= tau*tmp;
-          m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
-        }
-        if (m_computeQZ)
-        {
-          // Z
-          Map<Matrix<Scalar,1,Dynamic> > tmp(m_workspace.data(),dim);
-          tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k));
-          tmp += m_Z.row(k+2);
-          m_Z.row(k+2) -= tau*tmp;
-          m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp);
-        }
-        m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0);
-
-        // Z_{k2} to annihilate T(k+1,k)
-        G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k));
-        m_S.applyOnTheRight(k+1,k,G);
-        m_T.applyOnTheRight(k+1,k,G);
-        // update Z
-        if (m_computeQZ)
-          m_Z.applyOnTheLeft(k+1,k,G.adjoint());
-        m_T.coeffRef(k+1,k) = Scalar(0.0);
-
-        // update x,y,z
-        x = m_S.coeff(k+1,k);
-        y = m_S.coeff(k+2,k);
-        if (k < l-2)
-          z = m_S.coeff(k+3,k);
-      } // loop over k
-
-      // Q_{n-1} to annihilate y = S(l,l-2)
-      G.makeGivens(x,y);
-      m_S.applyOnTheLeft(l-1,l,G.adjoint());
-      m_T.applyOnTheLeft(l-1,l,G.adjoint());
-      if (m_computeQZ)
-        m_Q.applyOnTheRight(l-1,l,G);
-      m_S.coeffRef(l,l-2) = Scalar(0.0);
-
-      // Z_{n-1} to annihilate T(l,l-1)
-      G.makeGivens(m_T.coeff(l,l),m_T.coeff(l,l-1));
-      m_S.applyOnTheRight(l,l-1,G);
-      m_T.applyOnTheRight(l,l-1,G);
-      if (m_computeQZ)
-        m_Z.applyOnTheLeft(l,l-1,G.adjoint());
-      m_T.coeffRef(l,l-1) = Scalar(0.0);
-    }
-
-  template<typename MatrixType>
-    RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ)
-    {
-
-      const Index dim = A_in.cols();
-
-      eigen_assert (A_in.rows()==dim && A_in.cols()==dim 
-          && B_in.rows()==dim && B_in.cols()==dim 
-          && "Need square matrices of the same dimension");
-
-      m_isInitialized = true;
-      m_computeQZ = computeQZ;
-      m_S = A_in; m_T = B_in;
-      m_workspace.resize(dim*2);
-      m_global_iter = 0;
-
-      // entrance point: hessenberg triangular decomposition
-      hessenbergTriangular();
-      // compute L1 vector norms of T, S into m_normOfS, m_normOfT
-      computeNorms();
-
-      Index l = dim-1, 
-            f, 
-            local_iter = 0;
-
-      while (l>0 && local_iter<m_maxIters)
-      {
-        f = findSmallSubdiagEntry(l);
-        // now rows and columns f..l (including) decouple from the rest of the problem
-        if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0);
-        if (f == l) // One root found
-        {
-          l--;
-          local_iter = 0;
-        }
-        else if (f == l-1) // Two roots found
-        {
-          splitOffTwoRows(f);
-          l -= 2;
-          local_iter = 0;
-        }
-        else // No convergence yet
-        {
-          // if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations
-          Index z = findSmallDiagEntry(f,l);
-          if (z>=f)
-          {
-            // zero found
-            pushDownZero(z,f,l);
-          }
-          else
-          {
-            // We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg 
-            // and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to
-            // apply a QR-like iteration to rows and columns f..l.
-            step(f,l, local_iter);
-            local_iter++;
-            m_global_iter++;
-          }
-        }
-      }
-      // check if we converged before reaching iterations limit
-      m_info = (local_iter<m_maxIters) ? Success : NoConvergence;
-
-      // For each non triangular 2x2 diagonal block of S,
-      //    reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD.
-      // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors,
-      // and is in par with Lapack/Matlab QZ.
-      if(m_info==Success)
-      {
-        for(Index i=0; i<dim-1; ++i)
-        {
-          if(m_S.coeff(i+1, i) != Scalar(0))
-          {
-            JacobiRotation<Scalar> j_left, j_right;
-            internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right);
-
-            // Apply resulting Jacobi rotations
-            m_S.applyOnTheLeft(i,i+1,j_left);
-            m_S.applyOnTheRight(i,i+1,j_right);
-            m_T.applyOnTheLeft(i,i+1,j_left);
-            m_T.applyOnTheRight(i,i+1,j_right);
-            m_T(i+1,i) = m_T(i,i+1) = Scalar(0);
-
-            if(m_computeQZ) {
-              m_Q.applyOnTheRight(i,i+1,j_left.transpose());
-              m_Z.applyOnTheLeft(i,i+1,j_right.transpose());
-            }
-
-            i++;
-          }
-        }
-      }
-
-      return *this;
-    } // end compute
-
-} // end namespace Eigen
-
-#endif //EIGEN_REAL_QZ
diff --git a/eigen/Eigen/src/Eigenvalues/RealSchur.h b/eigen/Eigen/src/Eigenvalues/RealSchur.h
deleted file mode 100644
index 17ea903..0000000
--- a/eigen/Eigen/src/Eigenvalues/RealSchur.h
+++ /dev/null
@@ -1,546 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_REAL_SCHUR_H
-#define EIGEN_REAL_SCHUR_H
-
-#include "./HessenbergDecomposition.h"
-
-namespace Eigen { 
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class RealSchur
-  *
-  * \brief Performs a real Schur decomposition of a square matrix
-  *
-  * \tparam _MatrixType the type of the matrix of which we are computing the
-  * real Schur decomposition; this is expected to be an instantiation of the
-  * Matrix class template.
-  *
-  * Given a real square matrix A, this class computes the real Schur
-  * decomposition: \f$ A = U T U^T \f$ 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, \f$ U^{-1} = U^T \f$. 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 EigenSolver to compute
-  * the eigendecomposition of a matrix.
-  *
-  * Call the function compute() to compute the real Schur decomposition of a
-  * given matrix. Alternatively, you can use the RealSchur(const MatrixType&, bool)
-  * constructor which computes the real Schur decomposition at construction
-  * time. Once the decomposition is computed, you can use the matrixU() and
-  * matrixT() functions to retrieve the matrices U and T in the decomposition.
-  *
-  * The documentation of RealSchur(const MatrixType&, bool) contains an example
-  * of the typical use of this class.
-  *
-  * \note The implementation is adapted from
-  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
-  * Their code is based on EISPACK.
-  *
-  * \sa class ComplexSchur, class EigenSolver, class ComplexEigenSolver
-  */
-template<typename _MatrixType> class RealSchur
-{
-  public:
-    typedef _MatrixType MatrixType;
-    enum {
-      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-      Options = MatrixType::Options,
-      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-    };
-    typedef typename MatrixType::Scalar Scalar;
-    typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
-    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> 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.
-      */
-    explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
-            : m_matT(size, size),
-              m_matU(size, size),
-              m_workspaceVector(size),
-              m_hess(size),
-              m_isInitialized(false),
-              m_matUisUptodate(false),
-              m_maxIters(-1)
-    { }
-
-    /** \brief Constructor; computes real Schur decomposition of given matrix. 
-      * 
-      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
-      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
-      *
-      * This constructor calls compute() to compute the Schur decomposition.
-      *
-      * Example: \include RealSchur_RealSchur_MatrixType.cpp
-      * Output: \verbinclude RealSchur_RealSchur_MatrixType.out
-      */
-    template<typename InputType>
-    explicit RealSchur(const EigenBase<InputType>& matrix, bool computeU = true)
-            : m_matT(matrix.rows(),matrix.cols()),
-              m_matU(matrix.rows(),matrix.cols()),
-              m_workspaceVector(matrix.rows()),
-              m_hess(matrix.rows()),
-              m_isInitialized(false),
-              m_matUisUptodate(false),
-              m_maxIters(-1)
-    {
-      compute(matrix.derived(), computeU);
-    }
-
-    /** \brief Returns the orthogonal matrix in the Schur decomposition. 
-      *
-      * \returns A const reference to the matrix U.
-      *
-      * \pre Either the constructor RealSchur(const MatrixType&, bool) or the
-      * member function compute(const MatrixType&, bool) has been called before
-      * to compute the Schur decomposition of a matrix, and \p computeU was set
-      * to true (the default value).
-      *
-      * \sa RealSchur(const MatrixType&, bool) for an example
-      */
-    const MatrixType& matrixU() const
-    {
-      eigen_assert(m_isInitialized && "RealSchur is not initialized.");
-      eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
-      return m_matU;
-    }
-
-    /** \brief Returns the quasi-triangular matrix in the Schur decomposition. 
-      *
-      * \returns A const reference to the matrix T.
-      *
-      * \pre Either the constructor RealSchur(const MatrixType&, bool) or the
-      * member function compute(const MatrixType&, bool) has been called before
-      * to compute the Schur decomposition of a matrix.
-      *
-      * \sa RealSchur(const MatrixType&, bool) for an example
-      */
-    const MatrixType& matrixT() const
-    {
-      eigen_assert(m_isInitialized && "RealSchur is not initialized.");
-      return m_matT;
-    }
-  
-    /** \brief Computes Schur decomposition of given matrix. 
-      * 
-      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
-      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
-      * \returns    Reference to \c *this
-      *
-      * The Schur decomposition is computed by first reducing the matrix to
-      * Hessenberg form using the class HessenbergDecomposition. The Hessenberg
-      * matrix is then reduced to triangular form by performing Francis QR
-      * iterations with implicit double shift. The cost of computing the Schur
-      * decomposition depends on the number of iterations; as a rough guide, it
-      * may be taken to be \f$25n^3\f$ flops if \a computeU is true and
-      * \f$10n^3\f$ flops if \a computeU is false.
-      *
-      * Example: \include RealSchur_compute.cpp
-      * Output: \verbinclude RealSchur_compute.out
-      *
-      * \sa compute(const MatrixType&, bool, Index)
-      */
-    template<typename InputType>
-    RealSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
-
-    /** \brief Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T
-     *  \param[in] matrixH Matrix in Hessenberg form H
-     *  \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T
-     *  \param computeU Computes the matriX U of the Schur vectors
-     * \return Reference to \c *this
-     * 
-     *  This routine assumes that the matrix is already reduced in Hessenberg form matrixH
-     *  using either the class HessenbergDecomposition or another mean. 
-     *  It computes the upper quasi-triangular matrix T of the Schur decomposition of H
-     *  When computeU is true, this routine computes the matrix U such that 
-     *  A = U T U^T =  (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix
-     * 
-     * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix
-     * is not available, the user should give an identity matrix (Q.setIdentity())
-     * 
-     * \sa compute(const MatrixType&, bool)
-     */
-    template<typename HessMatrixType, typename OrthMatrixType>
-    RealSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU);
-    /** \brief Reports whether previous computation was successful.
-      *
-      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
-      */
-    ComputationInfo info() const
-    {
-      eigen_assert(m_isInitialized && "RealSchur is not initialized.");
-      return m_info;
-    }
-
-    /** \brief Sets the maximum number of iterations allowed. 
-      *
-      * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size
-      * of the matrix.
-      */
-    RealSchur& setMaxIterations(Index maxIters)
-    {
-      m_maxIters = maxIters;
-      return *this;
-    }
-
-    /** \brief Returns the maximum number of iterations. */
-    Index getMaxIterations()
-    {
-      return m_maxIters;
-    }
-
-    /** \brief Maximum number of iterations per row.
-      *
-      * If not otherwise specified, the maximum number of iterations is this number times the size of the
-      * matrix. It is currently set to 40.
-      */
-    static const int m_maxIterationsPerRow = 40;
-
-  private:
-    
-    MatrixType m_matT;
-    MatrixType m_matU;
-    ColumnVectorType m_workspaceVector;
-    HessenbergDecomposition<MatrixType> m_hess;
-    ComputationInfo m_info;
-    bool m_isInitialized;
-    bool m_matUisUptodate;
-    Index m_maxIters;
-
-    typedef Matrix<Scalar,3,1> Vector3s;
-
-    Scalar computeNormOfT();
-    Index findSmallSubdiagEntry(Index iu);
-    void splitOffTwoRows(Index iu, bool computeU, const 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, Scalar* workspace);
-};
-
-
-template<typename MatrixType>
-template<typename InputType>
-RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
-{
-  const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)();
-
-  eigen_assert(matrix.cols() == matrix.rows());
-  Index maxIters = m_maxIters;
-  if (maxIters == -1)
-    maxIters = m_maxIterationsPerRow * matrix.rows();
-
-  Scalar scale = matrix.derived().cwiseAbs().maxCoeff();
-  if(scale<considerAsZero)
-  {
-    m_matT.setZero(matrix.rows(),matrix.cols());
-    if(computeU)
-      m_matU.setIdentity(matrix.rows(),matrix.cols());
-    m_info = Success;
-    m_isInitialized = true;
-    m_matUisUptodate = computeU;
-    return *this;
-  }
-
-  // Step 1. Reduce to Hessenberg form
-  m_hess.compute(matrix.derived()/scale);
-
-  // Step 2. Reduce to real Schur form  
-  computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU);
-
-  m_matT *= scale;
-  
-  return *this;
-}
-template<typename MatrixType>
-template<typename HessMatrixType, typename OrthMatrixType>
-RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU)
-{
-  using std::abs;
-
-  m_matT = matrixH;
-  if(computeU)
-    m_matU = matrixQ;
-  
-  Index maxIters = m_maxIters;
-  if (maxIters == -1)
-    maxIters = m_maxIterationsPerRow * matrixH.rows();
-  m_workspaceVector.resize(m_matT.cols());
-  Scalar* workspace = &m_workspaceVector.coeffRef(0);
-
-  // The matrix m_matT is divided in three parts. 
-  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
-  // Rows il,...,iu is the part we are working on (the active window).
-  // Rows iu+1,...,end are already brought in triangular form.
-  Index iu = m_matT.cols() - 1;
-  Index iter = 0;      // iteration count for current eigenvalue
-  Index totalIter = 0; // iteration count for whole matrix
-  Scalar exshift(0);   // sum of exceptional shifts
-  Scalar norm = computeNormOfT();
-
-  if(norm!=Scalar(0))
-  {
-    while (iu >= 0)
-    {
-      Index il = findSmallSubdiagEntry(iu);
-
-      // Check for convergence
-      if (il == iu) // One root found
-      {
-        m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
-        if (iu > 0)
-          m_matT.coeffRef(iu, iu-1) = Scalar(0);
-        iu--;
-        iter = 0;
-      }
-      else if (il == iu-1) // Two roots found
-      {
-        splitOffTwoRows(iu, computeU, exshift);
-        iu -= 2;
-        iter = 0;
-      }
-      else // No convergence yet
-      {
-        // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
-        Vector3s firstHouseholderVector = Vector3s::Zero(), shiftInfo;
-        computeShift(iu, iter, exshift, shiftInfo);
-        iter = iter + 1;
-        totalIter = totalIter + 1;
-        if (totalIter > maxIters) break;
-        Index im;
-        initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
-        performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
-      }
-    }
-  }
-  if(totalIter <= maxIters)
-    m_info = Success;
-  else
-    m_info = NoConvergence;
-
-  m_isInitialized = true;
-  m_matUisUptodate = computeU;
-  return *this;
-}
-
-/** \internal Computes and returns vector L1 norm of T */
-template<typename MatrixType>
-inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
-{
-  const Index size = m_matT.cols();
-  // FIXME to be efficient the following would requires a triangular reduxion code
-  // Scalar norm = m_matT.upper().cwiseAbs().sum() 
-  //               + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
-  Scalar norm(0);
-  for (Index j = 0; j < size; ++j)
-    norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
-  return norm;
-}
-
-/** \internal Look for single small sub-diagonal element and returns its index */
-template<typename MatrixType>
-inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu)
-{
-  using std::abs;
-  Index res = iu;
-  while (res > 0)
-  {
-    Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
-    if (abs(m_matT.coeff(res,res-1)) <= NumTraits<Scalar>::epsilon() * s)
-      break;
-    res--;
-  }
-  return res;
-}
-
-/** \internal Update T given that rows iu-1 and iu decouple from the rest. */
-template<typename MatrixType>
-inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift)
-{
-  using std::sqrt;
-  using std::abs;
-  const Index size = m_matT.cols();
-
-  // The eigenvalues of the 2x2 matrix [a b; c d] are 
-  // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
-  Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
-  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);   // q = tr^2 / 4 - det = discr/4
-  m_matT.coeffRef(iu,iu) += exshift;
-  m_matT.coeffRef(iu-1,iu-1) += exshift;
-
-  if (q >= Scalar(0)) // Two real eigenvalues
-  {
-    Scalar z = sqrt(abs(q));
-    JacobiRotation<Scalar> rot;
-    if (p >= Scalar(0))
-      rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
-    else
-      rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
-
-    m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
-    m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
-    m_matT.coeffRef(iu, iu-1) = Scalar(0); 
-    if (computeU)
-      m_matU.applyOnTheRight(iu-1, iu, rot);
-  }
-
-  if (iu > 1) 
-    m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
-}
-
-/** \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */
-template<typename MatrixType>
-inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
-{
-  using std::sqrt;
-  using std::abs;
-  shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
-  shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
-  shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
-
-  // Wilkinson's original ad hoc shift
-  if (iter == 10)
-  {
-    exshift += shiftInfo.coeff(0);
-    for (Index i = 0; i <= iu; ++i)
-      m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
-    Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
-    shiftInfo.coeffRef(0) = Scalar(0.75) * s;
-    shiftInfo.coeffRef(1) = Scalar(0.75) * s;
-    shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
-  }
-
-  // MATLAB's new ad hoc shift
-  if (iter == 30)
-  {
-    Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
-    s = s * s + shiftInfo.coeff(2);
-    if (s > Scalar(0))
-    {
-      s = sqrt(s);
-      if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
-        s = -s;
-      s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
-      s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
-      exshift += s;
-      for (Index i = 0; i <= iu; ++i)
-        m_matT.coeffRef(i,i) -= s;
-      shiftInfo.setConstant(Scalar(0.964));
-    }
-  }
-}
-
-/** \internal Compute index im at which Francis QR step starts and the first Householder vector. */
-template<typename MatrixType>
-inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
-{
-  using std::abs;
-  Vector3s& v = firstHouseholderVector; // alias to save typing
-
-  for (im = iu-2; im >= il; --im)
-  {
-    const Scalar Tmm = m_matT.coeff(im,im);
-    const Scalar r = shiftInfo.coeff(0) - Tmm;
-    const Scalar s = shiftInfo.coeff(1) - Tmm;
-    v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
-    v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
-    v.coeffRef(2) = m_matT.coeff(im+2,im+1);
-    if (im == il) {
-      break;
-    }
-    const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
-    const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
-    if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
-      break;
-  }
-}
-
-/** \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */
-template<typename MatrixType>
-inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
-{
-  eigen_assert(im >= il);
-  eigen_assert(im <= iu-2);
-
-  const Index size = m_matT.cols();
-
-  for (Index k = im; k <= iu-2; ++k)
-  {
-    bool firstIteration = (k == im);
-
-    Vector3s v;
-    if (firstIteration)
-      v = firstHouseholderVector;
-    else
-      v = m_matT.template block<3,1>(k,k-1);
-
-    Scalar tau, beta;
-    Matrix<Scalar, 2, 1> ess;
-    v.makeHouseholder(ess, tau, beta);
-    
-    if (beta != Scalar(0)) // if v is not zero
-    {
-      if (firstIteration && k > il)
-        m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
-      else if (!firstIteration)
-        m_matT.coeffRef(k,k-1) = beta;
-
-      // These Householder transformations form the O(n^3) part of the algorithm
-      m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
-      m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
-      if (computeU)
-        m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
-    }
-  }
-
-  Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
-  Scalar tau, beta;
-  Matrix<Scalar, 1, 1> ess;
-  v.makeHouseholder(ess, tau, beta);
-
-  if (beta != Scalar(0)) // if v is not zero
-  {
-    m_matT.coeffRef(iu-1, iu-2) = beta;
-    m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
-    m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
-    if (computeU)
-      m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
-  }
-
-  // clean up pollution due to round-off errors
-  for (Index i = im+2; i <= iu; ++i)
-  {
-    m_matT.coeffRef(i,i-2) = Scalar(0);
-    if (i > im+2)
-      m_matT.coeffRef(i,i-3) = Scalar(0);
-  }
-}
-
-} // end namespace Eigen
-
-#endif // EIGEN_REAL_SCHUR_H
diff --git a/eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h b/eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h
deleted file mode 100644
index 2c22517..0000000
--- a/eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h
+++ /dev/null
@@ -1,77 +0,0 @@
-/*
- Copyright (c) 2011, Intel Corporation. All rights reserved.
-
- Redistribution and use in source and binary forms, with or without modification,
- are permitted provided that the following conditions are met:
-
- * Redistributions of source code must retain the above copyright notice, this
-   list of conditions and the following disclaimer.
- * Redistributions in binary form must reproduce the above copyright notice,
-   this list of conditions and the following disclaimer in the documentation
-   and/or other materials provided with the distribution.
- * Neither the name of Intel Corporation nor the names of its contributors may
-   be used to endorse or promote products derived from this software without
-   specific prior written permission.
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
- WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
- ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
- ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
- ********************************************************************************
- *   Content : Eigen bindings to LAPACKe
- *    Real Schur needed to real unsymmetrical eigenvalues/eigenvectors.
- ********************************************************************************
-*/
-
-#ifndef EIGEN_REAL_SCHUR_LAPACKE_H
-#define EIGEN_REAL_SCHUR_LAPACKE_H
-
-namespace Eigen { 
-
-/** \internal Specialization for the data types supported by LAPACKe */
-
-#define EIGEN_LAPACKE_SCHUR_REAL(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \
-template<> template<typename InputType> inline \
-RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
-RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, bool computeU) \
-{ \
-  eigen_assert(matrix.cols() == matrix.rows()); \
-\
-  lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), sdim, info; \
-  lapack_int matrix_order = LAPACKE_COLROW; \
-  char jobvs, sort='N'; \
-  LAPACK_##LAPACKE_PREFIX_U##_SELECT2 select = 0; \
-  jobvs = (computeU) ? 'V' : 'N'; \
-  m_matU.resize(n, n); \
-  lapack_int ldvs  = internal::convert_index<lapack_int>(m_matU.outerStride()); \
-  m_matT = matrix; \
-  lapack_int lda = internal::convert_index<lapack_int>(m_matT.outerStride()); \
-  Matrix<EIGTYPE, Dynamic, Dynamic> wr, wi; \
-  wr.resize(n, 1); wi.resize(n, 1); \
-  info = LAPACKE_##LAPACKE_PREFIX##gees( matrix_order, jobvs, sort, select, n, (LAPACKE_TYPE*)m_matT.data(), lda, &sdim, (LAPACKE_TYPE*)wr.data(), (LAPACKE_TYPE*)wi.data(), (LAPACKE_TYPE*)m_matU.data(), ldvs ); \
-  if(info == 0) \
-    m_info = Success; \
-  else \
-    m_info = NoConvergence; \
-\
-  m_isInitialized = true; \
-  m_matUisUptodate = computeU; \
-  return *this; \
-\
-}
-
-EIGEN_LAPACKE_SCHUR_REAL(double,   double, d, D, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_LAPACKE_SCHUR_REAL(float,    float,  s, S, ColMajor, LAPACK_COL_MAJOR)
-EIGEN_LAPACKE_SCHUR_REAL(double,   double, d, D, RowMajor, LAPACK_ROW_MAJOR)
-EIGEN_LAPACKE_SCHUR_REAL(float,    float,  s, S, RowMajor, LAPACK_ROW_MAJOR)
-
-} // end namespace Eigen
-
-#endif // EIGEN_REAL_SCHUR_LAPACKE_H
diff --git a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
deleted file mode 100644
index 9ddd553..0000000
--- a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ /dev/null
@@ -1,870 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
-#define EIGEN_SELFADJOINTEIGENSOLVER_H
-
-#include "./Tridiagonalization.h"
-
-namespace Eigen { 
-
-template<typename _MatrixType>
-class GeneralizedSelfAdjointEigenSolver;
-
-namespace internal {
-template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues;
-template<typename MatrixType, typename DiagType, typename SubDiagType>
-ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec);
-}
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class SelfAdjointEigenSolver
-  *
-  * \brief Computes eigenvalues and eigenvectors of selfadjoint matrices
-  *
-  * \tparam _MatrixType the type of the matrix of which we are computing the
-  * eigendecomposition; this is expected to be an instantiation of the Matrix
-  * class template.
-  *
-  * A matrix \f$ A \f$ 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
-  * selfadjoint matrix. These are the scalars \f$ \lambda \f$ and vectors
-  * \f$ v \f$ such that \f$ Av = \lambda v \f$.  The eigenvalues of a
-  * selfadjoint matrix are always real. If \f$ D \f$ is a diagonal matrix with
-  * the eigenvalues on the diagonal, and \f$ V \f$ is a matrix with the
-  * eigenvectors as its columns, then \f$ A = V D V^{-1} \f$ (for selfadjoint
-  * matrices, the matrix \f$ V \f$ is always invertible). This is called the
-  * eigendecomposition.
-  *
-  * The algorithm exploits the fact that the matrix is selfadjoint, making it
-  * faster and more accurate than the general purpose eigenvalue algorithms
-  * implemented in EigenSolver and ComplexEigenSolver.
-  *
-  * Only the \b lower \b triangular \b part of the input matrix is referenced.
-  *
-  * Call the function compute() to compute the eigenvalues and eigenvectors of
-  * a given matrix. Alternatively, you can use the
-  * SelfAdjointEigenSolver(const MatrixType&, int) constructor which computes
-  * the eigenvalues and eigenvectors at construction time. Once the eigenvalue
-  * and eigenvectors are computed, they can be retrieved with the eigenvalues()
-  * and eigenvectors() functions.
-  *
-  * The documentation for SelfAdjointEigenSolver(const MatrixType&, int)
-  * contains an example of the typical use of this class.
-  *
-  * To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and
-  * the likes, see the class GeneralizedSelfAdjointEigenSolver.
-  *
-  * \sa MatrixBase::eigenvalues(), class EigenSolver, class ComplexEigenSolver
-  */
-template<typename _MatrixType> class SelfAdjointEigenSolver
-{
-  public:
-
-    typedef _MatrixType MatrixType;
-    enum {
-      Size = MatrixType::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-      Options = MatrixType::Options,
-      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-    };
-    
-    /** \brief Scalar type for matrices of type \p _MatrixType. */
-    typedef typename MatrixType::Scalar Scalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-    
-    typedef Matrix<Scalar,Size,Size,ColMajor,MaxColsAtCompileTime,MaxColsAtCompileTime> EigenvectorsType;
-
-    /** \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>::Real RealScalar;
-    
-    friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>;
-
-    /** \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 typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
-    typedef Tridiagonalization<MatrixType> TridiagonalizationType;
-    typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType;
-
-    /** \brief Default constructor for fixed-size matrices.
-      *
-      * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via compute(). This constructor
-      * can only be used if \p _MatrixType is a fixed-size matrix; use
-      * SelfAdjointEigenSolver(Index) for dynamic-size matrices.
-      *
-      * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
-      */
-    EIGEN_DEVICE_FUNC
-    SelfAdjointEigenSolver()
-        : m_eivec(),
-          m_eivalues(),
-          m_subdiag(),
-          m_isInitialized(false)
-    { }
-
-    /** \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.
-      *
-      * This constructor is useful for dynamic-size matrices, when 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
-      */
-    EIGEN_DEVICE_FUNC
-    explicit SelfAdjointEigenSolver(Index size)
-        : m_eivec(size, size),
-          m_eivalues(size),
-          m_subdiag(size > 1 ? size - 1 : 1),
-          m_isInitialized(false)
-    {}
-
-    /** \brief Constructor; computes eigendecomposition of given matrix.
-      *
-      * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
-      *    be computed. Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
-      *
-      * This constructor calls compute(const MatrixType&, int) to compute the
-      * eigenvalues of the matrix \p matrix. The eigenvectors are computed if
-      * \p options equals #ComputeEigenvectors.
-      *
-      * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out
-      *
-      * \sa compute(const MatrixType&, int)
-      */
-    template<typename InputType>
-    EIGEN_DEVICE_FUNC
-    explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors)
-      : m_eivec(matrix.rows(), matrix.cols()),
-        m_eivalues(matrix.cols()),
-        m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
-        m_isInitialized(false)
-    {
-      compute(matrix.derived(), options);
-    }
-
-    /** \brief Computes eigendecomposition of given matrix.
-      *
-      * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
-      *    be computed. Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
-      * \returns    Reference to \c *this
-      *
-      * This function computes the eigenvalues of \p matrix.  The eigenvalues()
-      * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
-      * then the eigenvectors are also computed and can be retrieved by
-      * calling eigenvectors().
-      *
-      * This implementation uses a symmetric QR algorithm. The matrix is first
-      * reduced to tridiagonal form using the Tridiagonalization class. The
-      * tridiagonal matrix is then brought to diagonal form with implicit
-      * symmetric QR steps with Wilkinson shift. Details can be found in
-      * Section 8.3 of Golub \& Van Loan, <i>%Matrix Computations</i>.
-      *
-      * The cost of the computation is about \f$ 9n^3 \f$ if the eigenvectors
-      * are required and \f$ 4n^3/3 \f$ if they are not required.
-      *
-      * This method reuses the memory in the SelfAdjointEigenSolver object that
-      * was allocated when the object was constructed, if the size of the
-      * matrix does not change.
-      *
-      * Example: \include SelfAdjointEigenSolver_compute_MatrixType.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType.out
-      *
-      * \sa SelfAdjointEigenSolver(const MatrixType&, int)
-      */
-    template<typename InputType>
-    EIGEN_DEVICE_FUNC
-    SelfAdjointEigenSolver& compute(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors);
-    
-    /** \brief Computes eigendecomposition of given matrix using a closed-form algorithm
-      *
-      * This is a variant of compute(const MatrixType&, int options) which
-      * directly solves the underlying polynomial equation.
-      * 
-      * Currently only 2x2 and 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d).
-      * 
-      * This method is usually significantly faster than the QR iterative algorithm
-      * but it might also be less accurate. It is also worth noting that
-      * for 3x3 matrices it involves trigonometric operations which are
-      * not necessarily available for all scalar types.
-      * 
-      * For the 3x3 case, we observed the following worst case relative error regarding the eigenvalues:
-      *   - double: 1e-8
-      *   - float:  1e-3
-      *
-      * \sa compute(const MatrixType&, int options)
-      */
-    EIGEN_DEVICE_FUNC
-    SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors);
-
-    /**
-      *\brief Computes the eigen decomposition from a tridiagonal symmetric matrix
-      *
-      * \param[in] diag The vector containing the diagonal of the matrix.
-      * \param[in] subdiag The subdiagonal of the matrix.
-      * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
-      * \returns Reference to \c *this
-      *
-      * This function assumes that the matrix has been reduced to tridiagonal form.
-      *
-      * \sa compute(const MatrixType&, int) for more information
-      */
-    SelfAdjointEigenSolver& computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options=ComputeEigenvectors);
-
-    /** \brief Returns the eigenvectors of given matrix.
-      *
-      * \returns  A const reference to the matrix whose columns are the eigenvectors.
-      *
-      * \pre The eigenvectors have been computed before.
-      *
-      * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
-      * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
-      * eigenvectors are normalized to have (Euclidean) norm equal to one. If
-      * this object was used to solve the eigenproblem for the selfadjoint
-      * matrix \f$ A \f$, then the matrix returned by this function is the
-      * matrix \f$ V \f$ in the eigendecomposition \f$ A = V D V^{-1} \f$.
-      *
-      * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out
-      *
-      * \sa eigenvalues()
-      */
-    EIGEN_DEVICE_FUNC
-    const EigenvectorsType& eigenvectors() const
-    {
-      eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
-      return m_eivec;
-    }
-
-    /** \brief Returns the eigenvalues of given matrix.
-      *
-      * \returns A const reference to the column vector containing the eigenvalues.
-      *
-      * \pre The eigenvalues have been computed before.
-      *
-      * The eigenvalues are repeated according to their algebraic multiplicity,
-      * so there are as many eigenvalues as rows in the matrix. The eigenvalues
-      * are sorted in increasing order.
-      *
-      * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out
-      *
-      * \sa eigenvectors(), MatrixBase::eigenvalues()
-      */
-    EIGEN_DEVICE_FUNC
-    const RealVectorType& eigenvalues() const
-    {
-      eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      return m_eivalues;
-    }
-
-    /** \brief Computes the positive-definite square root of the matrix.
-      *
-      * \returns the positive-definite square root of the matrix
-      *
-      * \pre The eigenvalues and eigenvectors of a positive-definite matrix
-      * have been computed before.
-      *
-      * The square root of a positive-definite matrix \f$ A \f$ is the
-      * positive-definite matrix whose square equals \f$ A \f$. This function
-      * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the
-      * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$.
-      *
-      * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out
-      *
-      * \sa operatorInverseSqrt(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a>
-      */
-    EIGEN_DEVICE_FUNC
-    MatrixType operatorSqrt() const
-    {
-      eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
-      return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
-    }
-
-    /** \brief Computes the inverse square root of the matrix.
-      *
-      * \returns the inverse positive-definite square root of the matrix
-      *
-      * \pre The eigenvalues and eigenvectors of a positive-definite matrix
-      * have been computed before.
-      *
-      * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to
-      * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is
-      * cheaper than first computing the square root with operatorSqrt() and
-      * then its inverse with MatrixBase::inverse().
-      *
-      * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out
-      *
-      * \sa operatorSqrt(), MatrixBase::inverse(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a>
-      */
-    EIGEN_DEVICE_FUNC
-    MatrixType operatorInverseSqrt() const
-    {
-      eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
-      return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
-    }
-
-    /** \brief Reports whether previous computation was successful.
-      *
-      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
-      */
-    EIGEN_DEVICE_FUNC
-    ComputationInfo info() const
-    {
-      eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      return m_info;
-    }
-
-    /** \brief Maximum number of iterations.
-      *
-      * The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n
-      * denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK).
-      */
-    static const int m_maxIterations = 30;
-
-  protected:
-    static void check_template_parameters()
-    {
-      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
-    }
-    
-    EigenvectorsType m_eivec;
-    RealVectorType m_eivalues;
-    typename TridiagonalizationType::SubDiagonalType m_subdiag;
-    ComputationInfo m_info;
-    bool m_isInitialized;
-    bool m_eigenvectorsOk;
-};
-
-namespace internal {
-/** \internal
-  *
-  * \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  * Performs a QR step on a tridiagonal symmetric matrix represented as a
-  * pair of two vectors \a diag and \a subdiag.
-  *
-  * \param diag the diagonal part of the input selfadjoint tridiagonal matrix
-  * \param subdiag the sub-diagonal part of the input selfadjoint tridiagonal matrix
-  * \param start starting index of the submatrix to work on
-  * \param end last+1 index of the submatrix to work on
-  * \param matrixQ pointer to the column-major matrix holding the eigenvectors, can be 0
-  * \param n size of the input matrix
-  *
-  * For compilation efficiency reasons, this procedure does not use eigen expression
-  * for its arguments.
-  *
-  * Implemented from Golub's "Matrix Computations", algorithm 8.3.2:
-  * "implicit symmetric QR step with Wilkinson shift"
-  */
-template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
-EIGEN_DEVICE_FUNC
-static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
-}
-
-template<typename MatrixType>
-template<typename InputType>
-EIGEN_DEVICE_FUNC
-SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
-::compute(const EigenBase<InputType>& a_matrix, int options)
-{
-  check_template_parameters();
-  
-  const InputType &matrix(a_matrix.derived());
-  
-  using std::abs;
-  eigen_assert(matrix.cols() == matrix.rows());
-  eigen_assert((options&~(EigVecMask|GenEigMask))==0
-          && (options&EigVecMask)!=EigVecMask
-          && "invalid option parameter");
-  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
-  Index n = matrix.cols();
-  m_eivalues.resize(n,1);
-
-  if(n==1)
-  {
-    m_eivec = matrix;
-    m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0));
-    if(computeEigenvectors)
-      m_eivec.setOnes(n,n);
-    m_info = Success;
-    m_isInitialized = true;
-    m_eigenvectorsOk = computeEigenvectors;
-    return *this;
-  }
-
-  // declare some aliases
-  RealVectorType& diag = m_eivalues;
-  EigenvectorsType& mat = m_eivec;
-
-  // map the matrix coefficients to [-1:1] to avoid over- and underflow.
-  mat = matrix.template triangularView<Lower>();
-  RealScalar scale = mat.cwiseAbs().maxCoeff();
-  if(scale==RealScalar(0)) scale = RealScalar(1);
-  mat.template triangularView<Lower>() /= scale;
-  m_subdiag.resize(n-1);
-  internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
-
-  m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
-  
-  // scale back the eigen values
-  m_eivalues *= scale;
-
-  m_isInitialized = true;
-  m_eigenvectorsOk = computeEigenvectors;
-  return *this;
-}
-
-template<typename MatrixType>
-SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
-::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options)
-{
-  //TODO : Add an option to scale the values beforehand
-  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
-
-  m_eivalues = diag;
-  m_subdiag = subdiag;
-  if (computeEigenvectors)
-  {
-    m_eivec.setIdentity(diag.size(), diag.size());
-  }
-  m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
-
-  m_isInitialized = true;
-  m_eigenvectorsOk = computeEigenvectors;
-  return *this;
-}
-
-namespace internal {
-/**
-  * \internal
-  * \brief Compute the eigendecomposition from a tridiagonal matrix
-  *
-  * \param[in,out] diag : On input, the diagonal of the matrix, on output the eigenvalues
-  * \param[in,out] subdiag : The subdiagonal part of the matrix (entries are modified during the decomposition)
-  * \param[in] maxIterations : the maximum number of iterations
-  * \param[in] computeEigenvectors : whether the eigenvectors have to be computed or not
-  * \param[out] eivec : The matrix to store the eigenvectors if computeEigenvectors==true. Must be allocated on input.
-  * \returns \c Success or \c NoConvergence
-  */
-template<typename MatrixType, typename DiagType, typename SubDiagType>
-ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec)
-{
-  using std::abs;
-
-  ComputationInfo info;
-  typedef typename MatrixType::Scalar Scalar;
-
-  Index n = diag.size();
-  Index end = n-1;
-  Index start = 0;
-  Index iter = 0; // total number of iterations
-  
-  typedef typename DiagType::RealScalar RealScalar;
-  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
-  const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
-  
-  while (end>0)
-  {
-    for (Index i = start; i<end; ++i)
-      if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1])),precision) || abs(subdiag[i]) <= considerAsZero)
-        subdiag[i] = 0;
-
-    // find the largest unreduced block
-    while (end>0 && subdiag[end-1]==RealScalar(0))
-    {
-      end--;
-    }
-    if (end<=0)
-      break;
-
-    // if we spent too many iterations, we give up
-    iter++;
-    if(iter > maxIterations * n) break;
-
-    start = end - 1;
-    while (start>0 && subdiag[start-1]!=0)
-      start--;
-
-    internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
-  }
-  if (iter <= maxIterations * n)
-    info = Success;
-  else
-    info = NoConvergence;
-
-  // Sort eigenvalues and corresponding vectors.
-  // TODO make the sort optional ?
-  // TODO use a better sort algorithm !!
-  if (info == Success)
-  {
-    for (Index i = 0; i < n-1; ++i)
-    {
-      Index k;
-      diag.segment(i,n-i).minCoeff(&k);
-      if (k > 0)
-      {
-        std::swap(diag[i], diag[k+i]);
-        if(computeEigenvectors)
-          eivec.col(i).swap(eivec.col(k+i));
-      }
-    }
-  }
-  return info;
-}
-  
-template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
-{
-  EIGEN_DEVICE_FUNC
-  static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
-  { eig.compute(A,options); }
-};
-
-template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
-{
-  typedef typename SolverType::MatrixType MatrixType;
-  typedef typename SolverType::RealVectorType VectorType;
-  typedef typename SolverType::Scalar Scalar;
-  typedef typename SolverType::EigenvectorsType EigenvectorsType;
-  
-
-  /** \internal
-   * Computes the roots of the characteristic polynomial of \a m.
-   * For numerical stability m.trace() should be near zero and to avoid over- or underflow m should be normalized.
-   */
-  EIGEN_DEVICE_FUNC
-  static inline void computeRoots(const MatrixType& m, VectorType& roots)
-  {
-    EIGEN_USING_STD_MATH(sqrt)
-    EIGEN_USING_STD_MATH(atan2)
-    EIGEN_USING_STD_MATH(cos)
-    EIGEN_USING_STD_MATH(sin)
-    const Scalar s_inv3 = Scalar(1)/Scalar(3);
-    const Scalar s_sqrt3 = sqrt(Scalar(3));
-
-    // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0.  The
-    // eigenvalues are the roots to this equation, all guaranteed to be
-    // real-valued, because the matrix is symmetric.
-    Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(1,0)*m(2,0)*m(2,1) - m(0,0)*m(2,1)*m(2,1) - m(1,1)*m(2,0)*m(2,0) - m(2,2)*m(1,0)*m(1,0);
-    Scalar c1 = m(0,0)*m(1,1) - m(1,0)*m(1,0) + m(0,0)*m(2,2) - m(2,0)*m(2,0) + m(1,1)*m(2,2) - m(2,1)*m(2,1);
-    Scalar c2 = m(0,0) + m(1,1) + m(2,2);
-
-    // Construct the parameters used in classifying the roots of the equation
-    // and in solving the equation for the roots in closed form.
-    Scalar c2_over_3 = c2*s_inv3;
-    Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3;
-    a_over_3 = numext::maxi(a_over_3, Scalar(0));
-
-    Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
-
-    Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
-    q = numext::maxi(q, Scalar(0));
-
-    // Compute the eigenvalues by solving for the roots of the polynomial.
-    Scalar rho = sqrt(a_over_3);
-    Scalar theta = atan2(sqrt(q),half_b)*s_inv3;  // since sqrt(q) > 0, atan2 is in [0, pi] and theta is in [0, pi/3]
-    Scalar cos_theta = cos(theta);
-    Scalar sin_theta = sin(theta);
-    // roots are already sorted, since cos is monotonically decreasing on [0, pi]
-    roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); // == 2*rho*cos(theta+2pi/3)
-    roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); // == 2*rho*cos(theta+ pi/3)
-    roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
-  }
-
-  EIGEN_DEVICE_FUNC
-  static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative)
-  {
-    using std::abs;
-    Index i0;
-    // Find non-zero column i0 (by construction, there must exist a non zero coefficient on the diagonal):
-    mat.diagonal().cwiseAbs().maxCoeff(&i0);
-    // mat.col(i0) is a good candidate for an orthogonal vector to the current eigenvector,
-    // so let's save it:
-    representative = mat.col(i0);
-    Scalar n0, n1;
-    VectorType c0, c1;
-    n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm();
-    n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm();
-    if(n0>n1) res = c0/std::sqrt(n0);
-    else      res = c1/std::sqrt(n1);
-
-    return true;
-  }
-
-  EIGEN_DEVICE_FUNC
-  static inline void run(SolverType& solver, const MatrixType& mat, int options)
-  {
-    eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
-    eigen_assert((options&~(EigVecMask|GenEigMask))==0
-            && (options&EigVecMask)!=EigVecMask
-            && "invalid option parameter");
-    bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
-    
-    EigenvectorsType& eivecs = solver.m_eivec;
-    VectorType& eivals = solver.m_eivalues;
-  
-    // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
-    Scalar shift = mat.trace() / Scalar(3);
-    // TODO Avoid this copy. Currently it is necessary to suppress bogus values when determining maxCoeff and for computing the eigenvectors later
-    MatrixType scaledMat = mat.template selfadjointView<Lower>();
-    scaledMat.diagonal().array() -= shift;
-    Scalar scale = scaledMat.cwiseAbs().maxCoeff();
-    if(scale > 0) scaledMat /= scale;   // TODO for scale==0 we could save the remaining operations
-
-    // compute the eigenvalues
-    computeRoots(scaledMat,eivals);
-
-    // compute the eigenvectors
-    if(computeEigenvectors)
-    {
-      if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
-      {
-        // All three eigenvalues are numerically the same
-        eivecs.setIdentity();
-      }
-      else
-      {
-        MatrixType tmp;
-        tmp = scaledMat;
-
-        // Compute the eigenvector of the most distinct eigenvalue
-        Scalar d0 = eivals(2) - eivals(1);
-        Scalar d1 = eivals(1) - eivals(0);
-        Index k(0), l(2);
-        if(d0 > d1)
-        {
-          numext::swap(k,l);
-          d0 = d1;
-        }
-
-        // Compute the eigenvector of index k
-        {
-          tmp.diagonal().array () -= eivals(k);
-          // By construction, 'tmp' is of rank 2, and its kernel corresponds to the respective eigenvector.
-          extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
-        }
-
-        // Compute eigenvector of index l
-        if(d0<=2*Eigen::NumTraits<Scalar>::epsilon()*d1)
-        {
-          // If d0 is too small, then the two other eigenvalues are numerically the same,
-          // and thus we only have to ortho-normalize the near orthogonal vector we saved above.
-          eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
-          eivecs.col(l).normalize();
-        }
-        else
-        {
-          tmp = scaledMat;
-          tmp.diagonal().array () -= eivals(l);
-
-          VectorType dummy;
-          extract_kernel(tmp, eivecs.col(l), dummy);
-        }
-
-        // Compute last eigenvector from the other two
-        eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
-      }
-    }
-
-    // Rescale back to the original size.
-    eivals *= scale;
-    eivals.array() += shift;
-    
-    solver.m_info = Success;
-    solver.m_isInitialized = true;
-    solver.m_eigenvectorsOk = computeEigenvectors;
-  }
-};
-
-// 2x2 direct eigenvalues decomposition, code from Hauke Heibel
-template<typename SolverType> 
-struct direct_selfadjoint_eigenvalues<SolverType,2,false>
-{
-  typedef typename SolverType::MatrixType MatrixType;
-  typedef typename SolverType::RealVectorType VectorType;
-  typedef typename SolverType::Scalar Scalar;
-  typedef typename SolverType::EigenvectorsType EigenvectorsType;
-  
-  EIGEN_DEVICE_FUNC
-  static inline void computeRoots(const MatrixType& m, VectorType& roots)
-  {
-    using std::sqrt;
-    const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0)));
-    const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
-    roots(0) = t1 - t0;
-    roots(1) = t1 + t0;
-  }
-  
-  EIGEN_DEVICE_FUNC
-  static inline void run(SolverType& solver, const MatrixType& mat, int options)
-  {
-    EIGEN_USING_STD_MATH(sqrt);
-    EIGEN_USING_STD_MATH(abs);
-    
-    eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
-    eigen_assert((options&~(EigVecMask|GenEigMask))==0
-            && (options&EigVecMask)!=EigVecMask
-            && "invalid option parameter");
-    bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
-    
-    EigenvectorsType& eivecs = solver.m_eivec;
-    VectorType& eivals = solver.m_eivalues;
-  
-    // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
-    Scalar shift = mat.trace() / Scalar(2);
-    MatrixType scaledMat = mat;
-    scaledMat.coeffRef(0,1) = mat.coeff(1,0);
-    scaledMat.diagonal().array() -= shift;
-    Scalar scale = scaledMat.cwiseAbs().maxCoeff();
-    if(scale > Scalar(0))
-      scaledMat /= scale;
-
-    // Compute the eigenvalues
-    computeRoots(scaledMat,eivals);
-
-    // compute the eigen vectors
-    if(computeEigenvectors)
-    {
-      if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon())
-      {
-        eivecs.setIdentity();
-      }
-      else
-      {
-        scaledMat.diagonal().array () -= eivals(1);
-        Scalar a2 = numext::abs2(scaledMat(0,0));
-        Scalar c2 = numext::abs2(scaledMat(1,1));
-        Scalar b2 = numext::abs2(scaledMat(1,0));
-        if(a2>c2)
-        {
-          eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
-          eivecs.col(1) /= sqrt(a2+b2);
-        }
-        else
-        {
-          eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
-          eivecs.col(1) /= sqrt(c2+b2);
-        }
-
-        eivecs.col(0) << eivecs.col(1).unitOrthogonal();
-      }
-    }
-
-    // Rescale back to the original size.
-    eivals *= scale;
-    eivals.array() += shift;
-
-    solver.m_info = Success;
-    solver.m_isInitialized = true;
-    solver.m_eigenvectorsOk = computeEigenvectors;
-  }
-};
-
-}
-
-template<typename MatrixType>
-EIGEN_DEVICE_FUNC
-SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
-::computeDirect(const MatrixType& matrix, int options)
-{
-  internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options);
-  return *this;
-}
-
-namespace internal {
-template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
-EIGEN_DEVICE_FUNC
-static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
-{
-  using std::abs;
-  RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
-  RealScalar e = subdiag[end-1];
-  // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
-  // underflow thus leading to inf/NaN values when using the following commented code:
-//   RealScalar e2 = numext::abs2(subdiag[end-1]);
-//   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
-  // This explain the following, somewhat more complicated, version:
-  RealScalar mu = diag[end];
-  if(td==RealScalar(0))
-    mu -= abs(e);
-  else
-  {
-    RealScalar e2 = numext::abs2(subdiag[end-1]);
-    RealScalar h = numext::hypot(td,e);
-    if(e2==RealScalar(0)) mu -= (e / (td + (td>RealScalar(0) ? RealScalar(1) : RealScalar(-1)))) * (e / h);
-    else                  mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
-  }
-  
-  RealScalar x = diag[start] - mu;
-  RealScalar z = subdiag[start];
-  for (Index k = start; k < end; ++k)
-  {
-    JacobiRotation<RealScalar> rot;
-    rot.makeGivens(x, z);
-
-    // do T = G' T G
-    RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
-    RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
-
-    diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
-    diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
-    subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
-    
-
-    if (k > start)
-      subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
-
-    x = subdiag[k];
-
-    if (k < end - 1)
-    {
-      z = -rot.s() * subdiag[k+1];
-      subdiag[k + 1] = rot.c() * subdiag[k+1];
-    }
-    
-    // apply the givens rotation to the unit matrix Q = Q * G
-    if (matrixQ)
-    {
-      // FIXME if StorageOrder == RowMajor this operation is not very efficient
-      Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
-      q.applyOnTheRight(k,k+1,rot);
-    }
-  }
-}
-
-} // end namespace internal
-
-} // end namespace Eigen
-
-#endif // EIGEN_SELFADJOINTEIGENSOLVER_H
diff --git a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h
deleted file mode 100644
index b0c947d..0000000
--- a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h
+++ /dev/null
@@ -1,87 +0,0 @@
-/*
- Copyright (c) 2011, Intel Corporation. All rights reserved.
-
- Redistribution and use in source and binary forms, with or without modification,
- are permitted provided that the following conditions are met:
-
- * Redistributions of source code must retain the above copyright notice, this
-   list of conditions and the following disclaimer.
- * Redistributions in binary form must reproduce the above copyright notice,
-   this list of conditions and the following disclaimer in the documentation
-   and/or other materials provided with the distribution.
- * Neither the name of Intel Corporation nor the names of its contributors may
-   be used to endorse or promote products derived from this software without
-   specific prior written permission.
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
- WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
- ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
- ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
- ********************************************************************************
- *   Content : Eigen bindings to LAPACKe
- *    Self-adjoint eigenvalues/eigenvectors.
- ********************************************************************************
-*/
-
-#ifndef EIGEN_SAEIGENSOLVER_LAPACKE_H
-#define EIGEN_SAEIGENSOLVER_LAPACKE_H
-
-namespace Eigen { 
-
-/** \internal Specialization for the data types supported by LAPACKe */
-
-#define EIGEN_LAPACKE_EIG_SELFADJ_2(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, EIGCOLROW ) \
-template<> template<typename InputType> inline \
-SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \
-SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, int options) \
-{ \
-  eigen_assert(matrix.cols() == matrix.rows()); \
-  eigen_assert((options&~(EigVecMask|GenEigMask))==0 \
-          && (options&EigVecMask)!=EigVecMask \
-          && "invalid option parameter"); \
-  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; \
-  lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), lda, info; \
-  m_eivalues.resize(n,1); \
-  m_subdiag.resize(n-1); \
-  m_eivec = matrix; \
-\
-  if(n==1) \
-  { \
-    m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0)); \
-    if(computeEigenvectors) m_eivec.setOnes(n,n); \
-    m_info = Success; \
-    m_isInitialized = true; \
-    m_eigenvectorsOk = computeEigenvectors; \
-    return *this; \
-  } \
-\
-  lda = internal::convert_index<lapack_int>(m_eivec.outerStride()); \
-  char jobz, uplo='L'/*, range='A'*/; \
-  jobz = computeEigenvectors ? 'V' : 'N'; \
-\
-  info = LAPACKE_##LAPACKE_NAME( LAPACK_COL_MAJOR, jobz, uplo, n, (LAPACKE_TYPE*)m_eivec.data(), lda, (LAPACKE_RTYPE*)m_eivalues.data() ); \
-  m_info = (info==0) ? Success : NoConvergence; \
-  m_isInitialized = true; \
-  m_eigenvectorsOk = computeEigenvectors; \
-  return *this; \
-}
-
-#define EIGEN_LAPACKE_EIG_SELFADJ(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME )              \
-        EIGEN_LAPACKE_EIG_SELFADJ_2(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, ColMajor )  \
-        EIGEN_LAPACKE_EIG_SELFADJ_2(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, RowMajor ) 
-
-EIGEN_LAPACKE_EIG_SELFADJ(double,   double,                double, dsyev)
-EIGEN_LAPACKE_EIG_SELFADJ(float,    float,                 float,  ssyev)
-EIGEN_LAPACKE_EIG_SELFADJ(dcomplex, lapack_complex_double, double, zheev)
-EIGEN_LAPACKE_EIG_SELFADJ(scomplex, lapack_complex_float,  float,  cheev)
-
-} // end namespace Eigen
-
-#endif // EIGEN_SAEIGENSOLVER_H
diff --git a/eigen/Eigen/src/Eigenvalues/Tridiagonalization.h b/eigen/Eigen/src/Eigenvalues/Tridiagonalization.h
deleted file mode 100644
index 1d102c1..0000000
--- a/eigen/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ /dev/null
@@ -1,556 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#ifndef EIGEN_TRIDIAGONALIZATION_H
-#define EIGEN_TRIDIAGONALIZATION_H
-
-namespace Eigen { 
-
-namespace internal {
-  
-template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
-template<typename MatrixType>
-struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
-  : public traits<typename MatrixType::PlainObject>
-{
-  typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix?
-  enum { Flags = 0 };
-};
-
-template<typename MatrixType, typename CoeffVectorType>
-void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
-}
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  *
-  * \class Tridiagonalization
-  *
-  * \brief Tridiagonal decomposition of a selfadjoint matrix
-  *
-  * \tparam _MatrixType the type of the matrix of which we are computing the
-  * tridiagonal decomposition; this is expected to be an instantiation of the
-  * Matrix class template.
-  *
-  * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that:
-  * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix.
-  *
-  * A tridiagonal matrix is a matrix which has nonzero elements only on the
-  * main diagonal and the first diagonal below and above it. The Hessenberg
-  * decomposition of a selfadjoint matrix is in fact a tridiagonal
-  * decomposition. This class is used in SelfAdjointEigenSolver to compute the
-  * eigenvalues and eigenvectors of a selfadjoint matrix.
-  *
-  * Call the function compute() to compute the tridiagonal decomposition of a
-  * given matrix. Alternatively, you can use the Tridiagonalization(const MatrixType&)
-  * constructor which computes the tridiagonal Schur decomposition at
-  * construction time. Once the decomposition is computed, you can use the
-  * matrixQ() and matrixT() functions to retrieve the matrices Q and T in the
-  * decomposition.
-  *
-  * The documentation of Tridiagonalization(const MatrixType&) contains an
-  * example of the typical use of this class.
-  *
-  * \sa class HessenbergDecomposition, class SelfAdjointEigenSolver
-  */
-template<typename _MatrixType> class Tridiagonalization
-{
-  public:
-
-    /** \brief Synonym for the template parameter \p _MatrixType. */
-    typedef _MatrixType MatrixType;
-
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
-
-    enum {
-      Size = MatrixType::RowsAtCompileTime,
-      SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
-      Options = MatrixType::Options,
-      MaxSize = MatrixType::MaxRowsAtCompileTime,
-      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
-    };
-
-    typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
-    typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
-    typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
-    typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
-    typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
-
-    typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
-              typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
-              const Diagonal<const MatrixType>
-            >::type DiagonalReturnType;
-
-    typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
-              typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,
-              const Diagonal<const MatrixType, -1>
-            >::type SubDiagonalReturnType;
-
-    /** \brief Return type of matrixQ() */
-    typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
-
-    /** \brief Default constructor.
-      *
-      * \param [in]  size  Positive integer, size of the matrix whose tridiagonal
-      * 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.
-      */
-    explicit Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
-      : m_matrix(size,size),
-        m_hCoeffs(size > 1 ? size-1 : 1),
-        m_isInitialized(false)
-    {}
-
-    /** \brief Constructor; computes tridiagonal decomposition of given matrix.
-      *
-      * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
-      * is to be computed.
-      *
-      * This constructor calls compute() to compute the tridiagonal decomposition.
-      *
-      * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
-      * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
-      */
-    template<typename InputType>
-    explicit Tridiagonalization(const EigenBase<InputType>& matrix)
-      : m_matrix(matrix.derived()),
-        m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
-        m_isInitialized(false)
-    {
-      internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
-      m_isInitialized = true;
-    }
-
-    /** \brief Computes tridiagonal decomposition of given matrix.
-      *
-      * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
-      * is to be computed.
-      * \returns    Reference to \c *this
-      *
-      * The tridiagonal decomposition is computed by bringing the columns of
-      * the matrix successively in the required form using Householder
-      * reflections. The cost is \f$ 4n^3/3 \f$ flops, where \f$ n \f$ denotes
-      * the size of the given matrix.
-      *
-      * This method reuses of the allocated data in the Tridiagonalization
-      * object, if the size of the matrix does not change.
-      *
-      * Example: \include Tridiagonalization_compute.cpp
-      * Output: \verbinclude Tridiagonalization_compute.out
-      */
-    template<typename InputType>
-    Tridiagonalization& compute(const EigenBase<InputType>& matrix)
-    {
-      m_matrix = matrix.derived();
-      m_hCoeffs.resize(matrix.rows()-1, 1);
-      internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
-      m_isInitialized = true;
-      return *this;
-    }
-
-    /** \brief Returns the Householder coefficients.
-      *
-      * \returns a const reference to the vector of Householder coefficients
-      *
-      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
-      * the member function compute(const MatrixType&) has been called before
-      * to compute the tridiagonal decomposition of a matrix.
-      *
-      * The Householder coefficients allow the reconstruction of the matrix
-      * \f$ Q \f$ in the tridiagonal decomposition from the packed data.
-      *
-      * Example: \include Tridiagonalization_householderCoefficients.cpp
-      * Output: \verbinclude Tridiagonalization_householderCoefficients.out
-      *
-      * \sa packedMatrix(), \ref Householder_Module "Householder module"
-      */
-    inline CoeffVectorType householderCoefficients() const
-    {
-      eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
-      return m_hCoeffs;
-    }
-
-    /** \brief Returns the internal representation of the decomposition
-      *
-      *	\returns a const reference to a matrix with the internal representation
-      *	         of the decomposition.
-      *
-      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
-      * the member function compute(const MatrixType&) has been called before
-      * to compute the tridiagonal decomposition of a matrix.
-      *
-      * The returned matrix contains the following information:
-      *  - the strict upper triangular part is equal to the input matrix A.
-      *  - the diagonal and lower sub-diagonal represent the real tridiagonal
-      *    symmetric matrix T.
-      *  - the rest of the lower part contains the Householder vectors that,
-      *    combined with Householder coefficients returned by
-      *    householderCoefficients(), allows to reconstruct the matrix Q as
-      *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
-      *    Here, the matrices \f$ H_i \f$ are the Householder transformations
-      *       \f$ H_i = (I - h_i v_i v_i^T) \f$
-      *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
-      *    \f$ v_i \f$ is the Householder vector defined by
-      *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
-      *    with M the matrix returned by this function.
-      *
-      * See LAPACK for further details on this packed storage.
-      *
-      * Example: \include Tridiagonalization_packedMatrix.cpp
-      * Output: \verbinclude Tridiagonalization_packedMatrix.out
-      *
-      * \sa householderCoefficients()
-      */
-    inline const MatrixType& packedMatrix() const
-    {
-      eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
-      return m_matrix;
-    }
-
-    /** \brief Returns the unitary matrix Q in the decomposition
-      *
-      * \returns object representing the matrix Q
-      *
-      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
-      * the member function compute(const MatrixType&) has been called before
-      * to compute the tridiagonal decomposition of a matrix.
-      *
-      * This function returns a light-weight object of template class
-      * HouseholderSequence. You can either apply it directly to a matrix or
-      * you can convert it to a matrix of type #MatrixType.
-      *
-      * \sa Tridiagonalization(const MatrixType&) for an example,
-      *     matrixT(), class HouseholderSequence
-      */
-    HouseholderSequenceType matrixQ() const
-    {
-      eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
-      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
-             .setLength(m_matrix.rows() - 1)
-             .setShift(1);
-    }
-
-    /** \brief Returns an expression of the tridiagonal matrix T in the decomposition
-      *
-      * \returns expression object representing the matrix T
-      *
-      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
-      * the member function compute(const MatrixType&) has been called before
-      * to compute the tridiagonal decomposition of a matrix.
-      *
-      * Currently, this function can be used to extract the matrix T from internal
-      * data and copy it to a dense matrix object. In most cases, it may be
-      * sufficient to directly use the packed matrix or the vector expressions
-      * returned by diagonal() and subDiagonal() instead of creating a new
-      * dense copy matrix with this function.
-      *
-      * \sa Tridiagonalization(const MatrixType&) for an example,
-      * matrixQ(), packedMatrix(), diagonal(), subDiagonal()
-      */
-    MatrixTReturnType matrixT() const
-    {
-      eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
-      return MatrixTReturnType(m_matrix.real());
-    }
-
-    /** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition.
-      *
-      * \returns expression representing the diagonal of T
-      *
-      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
-      * the member function compute(const MatrixType&) has been called before
-      * to compute the tridiagonal decomposition of a matrix.
-      *
-      * Example: \include Tridiagonalization_diagonal.cpp
-      * Output: \verbinclude Tridiagonalization_diagonal.out
-      *
-      * \sa matrixT(), subDiagonal()
-      */
-    DiagonalReturnType diagonal() const;
-
-    /** \brief Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
-      *
-      * \returns expression representing the subdiagonal of T
-      *
-      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
-      * the member function compute(const MatrixType&) has been called before
-      * to compute the tridiagonal decomposition of a matrix.
-      *
-      * \sa diagonal() for an example, matrixT()
-      */
-    SubDiagonalReturnType subDiagonal() const;
-
-  protected:
-
-    MatrixType m_matrix;
-    CoeffVectorType m_hCoeffs;
-    bool m_isInitialized;
-};
-
-template<typename MatrixType>
-typename Tridiagonalization<MatrixType>::DiagonalReturnType
-Tridiagonalization<MatrixType>::diagonal() const
-{
-  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
-  return m_matrix.diagonal().real();
-}
-
-template<typename MatrixType>
-typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
-Tridiagonalization<MatrixType>::subDiagonal() const
-{
-  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
-  return m_matrix.template diagonal<-1>().real();
-}
-
-namespace internal {
-
-/** \internal
-  * Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place.
-  *
-  * \param[in,out] matA On input the selfadjoint matrix. Only the \b lower triangular part is referenced.
-  *                     On output, the strict upper part is left unchanged, and the lower triangular part
-  *                     represents the T and Q matrices in packed format has detailed below.
-  * \param[out]    hCoeffs returned Householder coefficients (see below)
-  *
-  * On output, the tridiagonal selfadjoint matrix T is stored in the diagonal
-  * and lower sub-diagonal of the matrix \a matA.
-  * The unitary matrix Q is represented in a compact way as a product of
-  * Householder reflectors \f$ H_i \f$ such that:
-  *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
-  * The Householder reflectors are defined as
-  *       \f$ H_i = (I - h_i v_i v_i^T) \f$
-  * where \f$ h_i = hCoeffs[i]\f$ is the \f$ i \f$th Householder coefficient and
-  * \f$ v_i \f$ is the Householder vector defined by
-  *       \f$ v_i = [ 0, \ldots, 0, 1, matA(i+2,i), \ldots, matA(N-1,i) ]^T \f$.
-  *
-  * Implemented from Golub's "Matrix Computations", algorithm 8.3.1.
-  *
-  * \sa Tridiagonalization::packedMatrix()
-  */
-template<typename MatrixType, typename CoeffVectorType>
-void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
-{
-  using numext::conj;
-  typedef typename MatrixType::Scalar Scalar;
-  typedef typename MatrixType::RealScalar RealScalar;
-  Index n = matA.rows();
-  eigen_assert(n==matA.cols());
-  eigen_assert(n==hCoeffs.size()+1 || n==1);
-  
-  for (Index i = 0; i<n-1; ++i)
-  {
-    Index remainingSize = n-i-1;
-    RealScalar beta;
-    Scalar h;
-    matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
-
-    // Apply similarity transformation to remaining columns,
-    // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
-    matA.col(i).coeffRef(i+1) = 1;
-
-    hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
-                                  * (conj(h) * matA.col(i).tail(remainingSize)));
-
-    hCoeffs.tail(n-i-1) += (conj(h)*RealScalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1);
-
-    matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
-      .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));
-
-    matA.col(i).coeffRef(i+1) = beta;
-    hCoeffs.coeffRef(i) = h;
-  }
-}
-
-// forward declaration, implementation at the end of this file
-template<typename MatrixType,
-         int Size=MatrixType::ColsAtCompileTime,
-         bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
-struct tridiagonalization_inplace_selector;
-
-/** \brief Performs a full tridiagonalization in place
-  *
-  * \param[in,out]  mat  On input, the selfadjoint matrix whose tridiagonal
-  *    decomposition is to be computed. Only the lower triangular part referenced.
-  *    The rest is left unchanged. On output, the orthogonal matrix Q
-  *    in the decomposition if \p extractQ is true.
-  * \param[out]  diag  The diagonal of the tridiagonal matrix T in the
-  *    decomposition.
-  * \param[out]  subdiag  The subdiagonal of the tridiagonal matrix T in
-  *    the decomposition.
-  * \param[in]  extractQ  If true, the orthogonal matrix Q in the
-  *    decomposition is computed and stored in \p mat.
-  *
-  * Computes the tridiagonal decomposition of the selfadjoint matrix \p mat in place
-  * such that \f$ mat = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real
-  * symmetric tridiagonal matrix.
-  *
-  * The tridiagonal matrix T is passed to the output parameters \p diag and \p subdiag. If
-  * \p extractQ is true, then the orthogonal matrix Q is passed to \p mat. Otherwise the lower
-  * part of the matrix \p mat is destroyed.
-  *
-  * The vectors \p diag and \p subdiag are not resized. The function
-  * assumes that they are already of the correct size. The length of the
-  * vector \p diag should equal the number of rows in \p mat, and the
-  * length of the vector \p subdiag should be one left.
-  *
-  * This implementation contains an optimized path for 3-by-3 matrices
-  * which is especially useful for plane fitting.
-  *
-  * \note Currently, it requires two temporary vectors to hold the intermediate
-  * Householder coefficients, and to reconstruct the matrix Q from the Householder
-  * reflectors.
-  *
-  * Example (this uses the same matrix as the example in
-  *    Tridiagonalization::Tridiagonalization(const MatrixType&)):
-  *    \include Tridiagonalization_decomposeInPlace.cpp
-  * Output: \verbinclude Tridiagonalization_decomposeInPlace.out
-  *
-  * \sa class Tridiagonalization
-  */
-template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
-void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
-{
-  eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
-  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
-}
-
-/** \internal
-  * General full tridiagonalization
-  */
-template<typename MatrixType, int Size, bool IsComplex>
-struct tridiagonalization_inplace_selector
-{
-  typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
-  typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
-  template<typename DiagonalType, typename SubDiagonalType>
-  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
-  {
-    CoeffVectorType hCoeffs(mat.cols()-1);
-    tridiagonalization_inplace(mat,hCoeffs);
-    diag = mat.diagonal().real();
-    subdiag = mat.template diagonal<-1>().real();
-    if(extractQ)
-      mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
-            .setLength(mat.rows() - 1)
-            .setShift(1);
-  }
-};
-
-/** \internal
-  * Specialization for 3x3 real matrices.
-  * Especially useful for plane fitting.
-  */
-template<typename MatrixType>
-struct tridiagonalization_inplace_selector<MatrixType,3,false>
-{
-  typedef typename MatrixType::Scalar Scalar;
-  typedef typename MatrixType::RealScalar RealScalar;
-
-  template<typename DiagonalType, typename SubDiagonalType>
-  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
-  {
-    using std::sqrt;
-    const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
-    diag[0] = mat(0,0);
-    RealScalar v1norm2 = numext::abs2(mat(2,0));
-    if(v1norm2 <= tol)
-    {
-      diag[1] = mat(1,1);
-      diag[2] = mat(2,2);
-      subdiag[0] = mat(1,0);
-      subdiag[1] = mat(2,1);
-      if (extractQ)
-        mat.setIdentity();
-    }
-    else
-    {
-      RealScalar beta = sqrt(numext::abs2(mat(1,0)) + v1norm2);
-      RealScalar invBeta = RealScalar(1)/beta;
-      Scalar m01 = mat(1,0) * invBeta;
-      Scalar m02 = mat(2,0) * invBeta;
-      Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
-      diag[1] = mat(1,1) + m02*q;
-      diag[2] = mat(2,2) - m02*q;
-      subdiag[0] = beta;
-      subdiag[1] = mat(2,1) - m01 * q;
-      if (extractQ)
-      {
-        mat << 1,   0,    0,
-               0, m01,  m02,
-               0, m02, -m01;
-      }
-    }
-  }
-};
-
-/** \internal
-  * Trivial specialization for 1x1 matrices
-  */
-template<typename MatrixType, bool IsComplex>
-struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
-{
-  typedef typename MatrixType::Scalar Scalar;
-
-  template<typename DiagonalType, typename SubDiagonalType>
-  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
-  {
-    diag(0,0) = numext::real(mat(0,0));
-    if(extractQ)
-      mat(0,0) = Scalar(1);
-  }
-};
-
-/** \internal
-  * \eigenvalues_module \ingroup Eigenvalues_Module
-  *
-  * \brief Expression type for return value of Tridiagonalization::matrixT()
-  *
-  * \tparam MatrixType type of underlying dense matrix
-  */
-template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
-: public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
-{
-  public:
-    /** \brief Constructor.
-      *
-      * \param[in] mat The underlying dense matrix
-      */
-    TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
-
-    template <typename ResultType>
-    inline void evalTo(ResultType& result) const
-    {
-      result.setZero();
-      result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
-      result.diagonal() = m_matrix.diagonal();
-      result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
-    }
-
-    Index rows() const { return m_matrix.rows(); }
-    Index cols() const { return m_matrix.cols(); }
-
-  protected:
-    typename MatrixType::Nested m_matrix;
-};
-
-} // end namespace internal
-
-} // end namespace Eigen
-
-#endif // EIGEN_TRIDIAGONALIZATION_H