rename the EigenSolver module to Eigenvalues

Eigen/src/Eigenvalues/CMakeLists.txt

@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_EIGENVALUES_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_EIGENVALUES_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigenvalues COMPONENT Devel
+  )

Eigen/src/Eigenvalues/ComplexEigenSolver.h

@@ -0,0 +1,148 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Claire Maurice
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
+#define EIGEN_COMPLEX_EIGEN_SOLVER_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+  * \nonstableyet
+  *
+  * \class ComplexEigenSolver
+  *
+  * \brief Eigen values/vectors solver for general complex matrices
+  *
+  * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+  *
+  * \sa class EigenSolver, class SelfAdjointEigenSolver
+  */
+template<typename _MatrixType> class ComplexEigenSolver
+{
+  public:
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+    typedef std::complex<RealScalar> Complex;
+    typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType;
+    typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType;
+
+    /**
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via ComplexEigenSolver::compute(const MatrixType&).
+    */
+    ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
+    {}
+
+    ComplexEigenSolver(const MatrixType& matrix)
+            : m_eivec(matrix.rows(),matrix.cols()),
+              m_eivalues(matrix.cols()),
+              m_isInitialized(false)
+    {
+      compute(matrix);
+    }
+
+    EigenvectorType eigenvectors(void) const
+    {
+      ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+      return m_eivec;
+    }
+
+    EigenvalueType eigenvalues() const
+    {
+      ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+      return m_eivalues;
+    }
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+    MatrixType m_eivec;
+    EigenvalueType m_eivalues;
+    bool m_isInitialized;
+};
+
+
+template<typename MatrixType>
+void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
+{
+  // this code is inspired from Jampack
+  assert(matrix.cols() == matrix.rows());
+  int n = matrix.cols();
+  m_eivalues.resize(n,1);
+
+  RealScalar eps = epsilon<RealScalar>();
+
+  // Reduce to complex Schur form
+  ComplexSchur<MatrixType> schur(matrix);
+
+  m_eivalues = schur.matrixT().diagonal();
+
+  m_eivec.setZero();
+
+  Scalar d2, z;
+  RealScalar norm = matrix.norm();
+
+  // compute the (normalized) eigenvectors
+  for(int k=n-1 ; k>=0 ; k--)
+  {
+    d2 = schur.matrixT().coeff(k,k);
+    m_eivec.coeffRef(k,k) = Scalar(1.0,0.0);
+    for(int i=k-1 ; i>=0 ; i--)
+    {
+      m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
+      if(k-i-1>0)
+        m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value();
+      z = schur.matrixT().coeff(i,i) - d2;
+      if(z==Scalar(0))
+        ei_real_ref(z) = eps * norm;
+      m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z;
+
+    }
+    m_eivec.col(k).normalize();
+  }
+
+  m_eivec = schur.matrixU() * m_eivec;
+  m_isInitialized = true;
+
+  // sort the eigenvalues
+  {
+    for (int i=0; i<n; i++)
+    {
+      int k;
+      m_eivalues.cwise().abs().end(n-i).minCoeff(&k);
+      if (k != 0)
+      {
+        k += i;
+        std::swap(m_eivalues[k],m_eivalues[i]);
+        m_eivec.col(i).swap(m_eivec.col(k));
+      }
+    }
+  }
+}
+
+
+
+#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H

Eigen/src/Eigenvalues/ComplexSchur.h

@@ -0,0 +1,237 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Claire Maurice
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_COMPLEX_SCHUR_H
+#define EIGEN_COMPLEX_SCHUR_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+  * \nonstableyet
+  *
+  * \class ComplexShur
+  *
+  * \brief Performs a complex Shur decomposition of a real or complex square matrix
+  *
+  */
+template<typename _MatrixType> class ComplexSchur
+{
+  public:
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+    typedef std::complex<RealScalar> Complex;
+    typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType;
+
+    /**
+      * \brief Default Constructor.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via ComplexSchur::compute(const MatrixType&).
+      */
+    ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false)
+    {}
+
+    ComplexSchur(const MatrixType& matrix)
+            : m_matT(matrix.rows(),matrix.cols()),
+              m_matU(matrix.rows(),matrix.cols()),
+              m_isInitialized(false)
+    {
+      compute(matrix);
+    }
+
+    ComplexMatrixType matrixU() const
+    {
+      ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
+      return m_matU;
+    }
+
+    ComplexMatrixType matrixT() const
+    {
+      ei_assert(m_isInitialized && "ComplexShur is not initialized.");
+      return m_matT;
+    }
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+    ComplexMatrixType m_matT, m_matU;
+    bool m_isInitialized;
+};
+
+/** Computes the principal value of the square root of the complex \a z. */
+template<typename RealScalar>
+std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
+{
+  RealScalar t, tre, tim;
+
+  t = ei_abs(z);
+
+  if (ei_abs(ei_real(z)) <= ei_abs(ei_imag(z)))
+  {
+    // No cancellation in these formulas
+    tre = ei_sqrt(0.5*(t + ei_real(z)));
+    tim = ei_sqrt(0.5*(t - ei_real(z)));
+  }
+  else
+  {
+    // Stable computation of the above formulas
+    if (z.real() > 0)
+    {
+      tre = t + z.real();
+      tim = ei_abs(ei_imag(z))*ei_sqrt(0.5/tre);
+      tre = ei_sqrt(0.5*tre);
+    }
+    else
+    {
+      tim = t - z.real();
+      tre = ei_abs(ei_imag(z))*ei_sqrt(0.5/tim);
+      tim = ei_sqrt(0.5*tim);
+    }
+  }
+  if(z.imag() < 0)
+    tim = -tim;
+
+  return (std::complex<RealScalar>(tre,tim));
+
+}
+
+template<typename MatrixType>
+void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
+{
+  // this code is inspired from Jampack
+  assert(matrix.cols() == matrix.rows());
+  int n = matrix.cols();
+
+  // Reduce to Hessenberg form
+  HessenbergDecomposition<MatrixType> hess(matrix);
+
+  m_matT = hess.matrixH();
+  m_matU = hess.matrixQ();
+
+  int iu = m_matT.cols() - 1;
+  int il;
+  RealScalar d,sd,sf;
+  Complex c,b,disc,r1,r2,kappa;
+
+  RealScalar eps = epsilon<RealScalar>();
+
+  int iter = 0;
+  while(true)
+  {
+    //locate the range in which to iterate
+    while(iu > 0)
+    {
+      d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1));
+      sd = ei_norm1(m_matT.coeffRef(iu,iu-1));
+
+      if(sd >= eps * d) break; // FIXME : precision criterion ??
+
+      m_matT.coeffRef(iu,iu-1) = Complex(0);
+      iter = 0;
+      --iu;
+    }
+    if(iu==0) break;
+    iter++;
+
+    if(iter >= 30)
+    {
+      // FIXME : what to do when iter==MAXITER ??
+      std::cerr << "MAXITER" << std::endl;
+      return;
+    }
+
+    il = iu-1;
+    while( il > 0 )
+    {
+      // check if the current 2x2 block on the diagonal is upper triangular
+      d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1));
+      sd = ei_norm1(m_matT.coeffRef(il,il-1));
+
+      if(sd < eps * d) break; // FIXME : precision criterion ??
+
+      --il;
+    }
+
+    if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0);
+
+    // compute the shift (the normalization by sf is to avoid under/overflow)
+    Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
+    sf = t.cwise().abs().sum();
+    t /= sf;
+
+    c = t.determinant();
+    b = t.diagonal().sum();
+
+    disc = ei_sqrt(b*b - RealScalar(4)*c);
+
+    r1 = (b+disc)/RealScalar(2);
+    r2 = (b-disc)/RealScalar(2);
+
+    if(ei_norm1(r1) > ei_norm1(r2))
+      r2 = c/r1;
+    else
+      r1 = c/r2;
+
+    if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1)))
+      kappa = sf * r1;
+    else
+      kappa = sf * r2;
+
+    // perform the QR step using Givens rotations
+    PlanarRotation<Complex> rot;
+    rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il));
+
+    for(int i=il ; i<iu ; i++)
+    {
+      m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint());
+      m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot);
+      m_matU.applyOnTheRight(i, i+1, rot);
+
+      if(i != iu-1)
+      {
+        int i1 = i+1;
+        int i2 = i+2;
+
+        rot.makeGivens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), &m_matT.coeffRef(i1,i));
+        m_matT.coeffRef(i2,i) = Complex(0);
+      }
+    }
+  }
+
+  // FIXME : is it necessary ?
+  /*
+  for(int i=0 ; i<n ; i++)
+    for(int j=0 ; j<n ; j++)
+    {
+      if(ei_abs(ei_real(m_matT.coeff(i,j))) < eps)
+        ei_real_ref(m_matT.coeffRef(i,j)) = 0;
+      if(ei_imag(ei_abs(m_matT.coeff(i,j))) < eps)
+        ei_imag_ref(m_matT.coeffRef(i,j)) = 0;
+    }
+  */
+
+  m_isInitialized = true;
+}
+
+#endif // EIGEN_COMPLEX_SCHUR_H

Eigen/src/Eigenvalues/EigenSolver.h

@@ -0,0 +1,723 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_EIGENSOLVER_H
+#define EIGEN_EIGENSOLVER_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+  * \nonstableyet
+  *
+  * \class EigenSolver
+  *
+  * \brief Eigen values/vectors solver for non selfadjoint matrices
+  *
+  * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+  *
+  * Currently it only support real matrices.
+  *
+  * \note this code was adapted from JAMA (public domain)
+  *
+  * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver
+  */
+template<typename _MatrixType> class EigenSolver
+{
+  public:
+
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+    typedef std::complex<RealScalar> Complex;
+    typedef Matrix<Complex, MatrixType::ColsAtCompileTime, 1> EigenvalueType;
+    typedef Matrix<Complex, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> EigenvectorType;
+    typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
+    typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
+
+    /**
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via EigenSolver::compute(const MatrixType&).
+    */
+    EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
+
+    EigenSolver(const MatrixType& matrix)
+      : m_eivec(matrix.rows(), matrix.cols()),
+        m_eivalues(matrix.cols()),
+        m_isInitialized(false)
+    {
+      compute(matrix);
+    }
+
+
+    EigenvectorType eigenvectors(void) const;
+
+    /** \returns a real matrix V of pseudo eigenvectors.
+      *
+      * Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks,
+      * and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D
+      * and V satisfy A*V = V*D.
+      *
+      * More precisely, if the diagonal matrix of the eigen values is:\n
+      * \f$
+      * \left[ \begin{array}{cccccc}
+      * u+iv &      &      &      &   &   \\
+      *      & u-iv &      &      &   &   \\
+      *      &      & a+ib &      &   &   \\
+      *      &      &      & a-ib &   &   \\
+      *      &      &      &      & x &   \\
+      *      &      &      &      &   & y \\
+      * \end{array} \right]
+      * \f$ \n
+      * then, we have:\n
+      * \f$
+      * D =\left[ \begin{array}{cccccc}
+      *  u & v &    &   &   &   \\
+      * -v & u &    &   &   &   \\
+      *    &   &  a & b &   &   \\
+      *    &   & -b & a &   &   \\
+      *    &   &    &   & x &   \\
+      *    &   &    &   &   & y \\
+      * \end{array} \right]
+      * \f$
+      *
+      * \sa pseudoEigenvalueMatrix()
+      */
+    const MatrixType& pseudoEigenvectors() const
+    {
+      ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+      return m_eivec;
+    }
+
+    MatrixType pseudoEigenvalueMatrix() const;
+
+    /** \returns the eigenvalues as a column vector */
+    EigenvalueType eigenvalues() const
+    {
+      ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+      return m_eivalues;
+    }
+
+    EigenSolver& compute(const MatrixType& matrix);
+
+  private:
+
+    void orthes(MatrixType& matH, RealVectorType& ort);
+    void hqr2(MatrixType& matH);
+
+  protected:
+    MatrixType m_eivec;
+    EigenvalueType m_eivalues;
+    bool m_isInitialized;
+};
+
+/** \returns the real block diagonal matrix D of the eigenvalues.
+  *
+  * See pseudoEigenvectors() for the details.
+  */
+template<typename MatrixType>
+MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
+{
+  ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+  int n = m_eivec.cols();
+  MatrixType matD = MatrixType::Zero(n,n);
+  for (int i=0; i<n; ++i)
+  {
+    if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i))))
+      matD.coeffRef(i,i) = ei_real(m_eivalues.coeff(i));
+    else
+    {
+      matD.template block<2,2>(i,i) <<  ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)),
+                                       -ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i));
+      ++i;
+    }
+  }
+  return matD;
+}
+
+/** \returns the normalized complex eigenvectors as a matrix of column vectors.
+  *
+  * \sa eigenvalues(), pseudoEigenvectors()
+  */
+template<typename MatrixType>
+typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigenvectors(void) const
+{
+  ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+  int n = m_eivec.cols();
+  EigenvectorType matV(n,n);
+  for (int j=0; j<n; ++j)
+  {
+    if (ei_isMuchSmallerThan(ei_abs(ei_imag(m_eivalues.coeff(j))), ei_abs(ei_real(m_eivalues.coeff(j)))))
+    {
+      // we have a real eigen value
+      matV.col(j) = m_eivec.col(j).template cast<Complex>();
+    }
+    else
+    {
+      // we have a pair of complex eigen values
+      for (int i=0; i<n; ++i)
+      {
+        matV.coeffRef(i,j)   = Complex(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
+        matV.coeffRef(i,j+1) = Complex(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>
+EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix)
+{
+  assert(matrix.cols() == matrix.rows());
+  int n = matrix.cols();
+  m_eivalues.resize(n,1);
+
+  MatrixType matH = matrix;
+  RealVectorType ort(n);
+
+  // Reduce to Hessenberg form.
+  orthes(matH, ort);
+
+  // Reduce Hessenberg to real Schur form.
+  hqr2(matH);
+
+  m_isInitialized = true;
+  return *this;
+}
+
+// Nonsymmetric reduction to Hessenberg form.
+template<typename MatrixType>
+void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
+{
+  //  This is derived from the Algol procedures orthes and ortran,
+  //  by Martin and Wilkinson, Handbook for Auto. Comp.,
+  //  Vol.ii-Linear Algebra, and the corresponding
+  //  Fortran subroutines in EISPACK.
+
+  int n = m_eivec.cols();
+  int low = 0;
+  int high = n-1;
+
+  for (int m = low+1; m <= high-1; ++m)
+  {
+    // Scale column.
+    RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum();
+    if (scale != 0.0)
+    {
+      // Compute Householder transformation.
+      RealScalar h = 0.0;
+      // FIXME could be rewritten, but this one looks better wrt cache
+      for (int i = high; i >= m; i--)
+      {
+        ort.coeffRef(i) = matH.coeff(i,m-1)/scale;
+        h += ort.coeff(i) * ort.coeff(i);
+      }
+      RealScalar g = ei_sqrt(h);
+      if (ort.coeff(m) > 0)
+        g = -g;
+      h = h - ort.coeff(m) * g;
+      ort.coeffRef(m) = ort.coeff(m) - g;
+
+      // Apply Householder similarity transformation
+      // H = (I-u*u'/h)*H*(I-u*u')/h)
+      int bSize = high-m+1;
+      matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h)
+        * (ort.segment(m, bSize).transpose() *  matH.block(m, m, bSize, n-m)));
+
+      matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize))
+        * (ort.segment(m, bSize)/h).transpose());
+
+      ort.coeffRef(m) = scale*ort.coeff(m);
+      matH.coeffRef(m,m-1) = scale*g;
+    }
+  }
+
+  // Accumulate transformations (Algol's ortran).
+  m_eivec.setIdentity();
+
+  for (int m = high-1; m >= low+1; m--)
+  {
+    if (matH.coeff(m,m-1) != 0.0)
+    {
+      ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m);
+
+      int bSize = high-m+1;
+      m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) /  (matH.coeff(m,m-1) * ort.coeff(m)))
+        * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) );
+    }
+  }
+}
+
+// Complex scalar division.
+template<typename Scalar>
+std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
+{
+  Scalar r,d;
+  if (ei_abs(yr) > ei_abs(yi))
+  {
+      r = yi/yr;
+      d = yr + r*yi;
+      return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
+  }
+  else
+  {
+      r = yr/yi;
+      d = yi + r*yr;
+      return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
+  }
+}
+
+
+// Nonsymmetric reduction from Hessenberg to real Schur form.
+template<typename MatrixType>
+void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
+{
+  //  This is derived from the Algol procedure hqr2,
+  //  by Martin and Wilkinson, Handbook for Auto. Comp.,
+  //  Vol.ii-Linear Algebra, and the corresponding
+  //  Fortran subroutine in EISPACK.
+
+  // Initialize
+  int nn = m_eivec.cols();
+  int n = nn-1;
+  int low = 0;
+  int high = nn-1;
+  Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
+  Scalar exshift = 0.0;
+  Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y;
+
+  // Store roots isolated by balanc and compute matrix norm
+  // FIXME to be efficient the following would requires a triangular reduxion code
+  // Scalar norm = matH.upper().cwise().abs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwise().abs().sum();
+  Scalar norm = 0.0;
+  for (int j = 0; j < nn; ++j)
+  {
+    // FIXME what's the purpose of the following since the condition is always false
+    if ((j < low) || (j > high))
+    {
+      m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0);
+    }
+    norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwise().abs().sum();
+  }
+
+  // Outer loop over eigenvalue index
+  int iter = 0;
+  while (n >= low)
+  {
+    // Look for single small sub-diagonal element
+    int l = n;
+    while (l > low)
+    {
+      s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l));
+      if (s == 0.0)
+          s = norm;
+      if (ei_abs(matH.coeff(l,l-1)) < eps * s)
+        break;
+      l--;
+    }
+
+    // Check for convergence
+    // One root found
+    if (l == n)
+    {
+      matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
+      m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0);
+      n--;
+      iter = 0;
+    }
+    else if (l == n-1) // Two roots found
+    {
+      w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
+      p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5);
+      q = p * p + w;
+      z = ei_sqrt(ei_abs(q));
+      matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
+      matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift;
+      x = matH.coeff(n,n);
+
+      // Scalar pair
+      if (q >= 0)
+      {
+        if (p >= 0)
+          z = p + z;
+        else
+          z = p - z;
+
+        m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0);
+        m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0);
+
+        x = matH.coeff(n,n-1);
+        s = ei_abs(x) + ei_abs(z);
+        p = x / s;
+        q = z / s;
+        r = ei_sqrt(p * p+q * q);
+        p = p / r;
+        q = q / r;
+
+        // Row modification
+        for (int j = n-1; j < nn; ++j)
+        {
+          z = matH.coeff(n-1,j);
+          matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j);
+          matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z;
+        }
+
+        // Column modification
+        for (int i = 0; i <= n; ++i)
+        {
+          z = matH.coeff(i,n-1);
+          matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n);
+          matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z;
+        }
+
+        // Accumulate transformations
+        for (int i = low; i <= high; ++i)
+        {
+          z = m_eivec.coeff(i,n-1);
+          m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n);
+          m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z;
+        }
+      }
+      else // Complex pair
+      {
+        m_eivalues.coeffRef(n-1) = Complex(x + p, z);
+        m_eivalues.coeffRef(n)   = Complex(x + p, -z);
+      }
+      n = n - 2;
+      iter = 0;
+    }
+    else // No convergence yet
+    {
+      // Form shift
+      x = matH.coeff(n,n);
+      y = 0.0;
+      w = 0.0;
+      if (l < n)
+      {
+          y = matH.coeff(n-1,n-1);
+          w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
+      }
+
+      // Wilkinson's original ad hoc shift
+      if (iter == 10)
+      {
+        exshift += x;
+        for (int i = low; i <= n; ++i)
+          matH.coeffRef(i,i) -= x;
+        s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2));
+        x = y = Scalar(0.75) * s;
+        w = Scalar(-0.4375) * s * s;
+      }
+
+      // MATLAB's new ad hoc shift
+      if (iter == 30)
+      {
+        s = Scalar((y - x) / 2.0);
+        s = s * s + w;
+        if (s > 0)
+        {
+          s = ei_sqrt(s);
+          if (y < x)
+            s = -s;
+          s = Scalar(x - w / ((y - x) / 2.0 + s));
+          for (int i = low; i <= n; ++i)
+            matH.coeffRef(i,i) -= s;
+          exshift += s;
+          x = y = w = Scalar(0.964);
+        }
+      }
+
+      iter = iter + 1;   // (Could check iteration count here.)
+
+      // Look for two consecutive small sub-diagonal elements
+      int m = n-2;
+      while (m >= l)
+      {
+        z = matH.coeff(m,m);
+        r = x - z;
+        s = y - z;
+        p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1);
+        q = matH.coeff(m+1,m+1) - z - r - s;
+        r = matH.coeff(m+2,m+1);
+        s = ei_abs(p) + ei_abs(q) + ei_abs(r);
+        p = p / s;
+        q = q / s;
+        r = r / s;
+        if (m == l) {
+          break;
+        }
+        if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) <
+          eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) +
+          ei_abs(matH.coeff(m+1,m+1)))))
+        {
+          break;
+        }
+        m--;
+      }
+
+      for (int i = m+2; i <= n; ++i)
+      {
+        matH.coeffRef(i,i-2) = 0.0;
+        if (i > m+2)
+          matH.coeffRef(i,i-3) = 0.0;
+      }
+
+      // Double QR step involving rows l:n and columns m:n
+      for (int k = m; k <= n-1; ++k)
+      {
+        int notlast = (k != n-1);
+        if (k != m) {
+          p = matH.coeff(k,k-1);
+          q = matH.coeff(k+1,k-1);
+          r = notlast ? matH.coeff(k+2,k-1) : Scalar(0);
+          x = ei_abs(p) + ei_abs(q) + ei_abs(r);
+          if (x != 0.0)
+          {
+            p = p / x;
+            q = q / x;
+            r = r / x;
+          }
+        }
+
+        if (x == 0.0)
+          break;
+
+        s = ei_sqrt(p * p + q * q + r * r);
+
+        if (p < 0)
+          s = -s;
+
+        if (s != 0)
+        {
+          if (k != m)
+            matH.coeffRef(k,k-1) = -s * x;
+          else if (l != m)
+            matH.coeffRef(k,k-1) = -matH.coeff(k,k-1);
+
+          p = p + s;
+          x = p / s;
+          y = q / s;
+          z = r / s;
+          q = q / p;
+          r = r / p;
+
+          // Row modification
+          for (int j = k; j < nn; ++j)
+          {
+            p = matH.coeff(k,j) + q * matH.coeff(k+1,j);
+            if (notlast)
+            {
+              p = p + r * matH.coeff(k+2,j);
+              matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z;
+            }
+            matH.coeffRef(k,j) = matH.coeff(k,j) - p * x;
+            matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y;
+          }
+
+          // Column modification
+          for (int i = 0; i <= std::min(n,k+3); ++i)
+          {
+            p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1);
+            if (notlast)
+            {
+              p = p + z * matH.coeff(i,k+2);
+              matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r;
+            }
+            matH.coeffRef(i,k) = matH.coeff(i,k) - p;
+            matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q;
+          }
+
+          // Accumulate transformations
+          for (int i = low; i <= high; ++i)
+          {
+            p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1);
+            if (notlast)
+            {
+              p = p + z * m_eivec.coeff(i,k+2);
+              m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r;
+            }
+            m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p;
+            m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q;
+          }
+        }  // (s != 0)
+      }  // k loop
+    }  // check convergence
+  }  // while (n >= low)
+
+  // Backsubstitute to find vectors of upper triangular form
+  if (norm == 0.0)
+  {
+      return;
+  }
+
+  for (n = nn-1; n >= 0; n--)
+  {
+    p = m_eivalues.coeff(n).real();
+    q = m_eivalues.coeff(n).imag();
+
+    // Scalar vector
+    if (q == 0)
+    {
+      int l = n;
+      matH.coeffRef(n,n) = 1.0;
+      for (int i = n-1; i >= 0; i--)
+      {
+        w = matH.coeff(i,i) - p;
+        r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1));
+
+        if (m_eivalues.coeff(i).imag() < 0.0)
+        {
+          z = w;
+          s = r;
+        }
+        else
+        {
+          l = i;
+          if (m_eivalues.coeff(i).imag() == 0.0)
+          {
+            if (w != 0.0)
+              matH.coeffRef(i,n) = -r / w;
+            else
+              matH.coeffRef(i,n) = -r / (eps * norm);
+          }
+          else // Solve real equations
+          {
+            x = matH.coeff(i,i+1);
+            y = matH.coeff(i+1,i);
+            q = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
+            t = (x * s - z * r) / q;
+            matH.coeffRef(i,n) = t;
+            if (ei_abs(x) > ei_abs(z))
+              matH.coeffRef(i+1,n) = (-r - w * t) / x;
+            else
+              matH.coeffRef(i+1,n) = (-s - y * t) / z;
+          }
+
+          // Overflow control
+          t = ei_abs(matH.coeff(i,n));
+          if ((eps * t) * t > 1)
+            matH.col(n).end(nn-i) /= t;
+        }
+      }
+    }
+    else if (q < 0) // Complex vector
+    {
+      std::complex<Scalar> cc;
+      int l = n-1;
+
+      // Last vector component imaginary so matrix is triangular
+      if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n)))
+      {
+        matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1);
+        matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1);
+      }
+      else
+      {
+        cc = cdiv<Scalar>(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q);
+        matH.coeffRef(n-1,n-1) = ei_real(cc);
+        matH.coeffRef(n-1,n) = ei_imag(cc);
+      }
+      matH.coeffRef(n,n-1) = 0.0;
+      matH.coeffRef(n,n) = 1.0;
+      for (int i = n-2; i >= 0; i--)
+      {
+        Scalar ra,sa,vr,vi;
+        ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1));
+        sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1));
+        w = matH.coeff(i,i) - p;
+
+        if (m_eivalues.coeff(i).imag() < 0.0)
+        {
+          z = w;
+          r = ra;
+          s = sa;
+        }
+        else
+        {
+          l = i;
+          if (m_eivalues.coeff(i).imag() == 0)
+          {
+            cc = cdiv(-ra,-sa,w,q);
+            matH.coeffRef(i,n-1) = ei_real(cc);
+            matH.coeffRef(i,n) = ei_imag(cc);
+          }
+          else
+          {
+            // Solve complex equations
+            x = matH.coeff(i,i+1);
+            y = matH.coeff(i+1,i);
+            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;
+            vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
+            if ((vr == 0.0) && (vi == 0.0))
+              vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(z));
+
+            cc= cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
+            matH.coeffRef(i,n-1) = ei_real(cc);
+            matH.coeffRef(i,n) = ei_imag(cc);
+            if (ei_abs(x) > (ei_abs(z) + ei_abs(q)))
+            {
+              matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x;
+              matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x;
+            }
+            else
+            {
+              cc = cdiv(-r-y*matH.coeff(i,n-1),-s-y*matH.coeff(i,n),z,q);
+              matH.coeffRef(i+1,n-1) = ei_real(cc);
+              matH.coeffRef(i+1,n) = ei_imag(cc);
+            }
+          }
+
+          // Overflow control
+          t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n)));
+          if ((eps * t) * t > 1)
+            matH.block(i, n-1, nn-i, 2) /= t;
+
+        }
+      }
+    }
+  }
+
+  // Vectors of isolated roots
+  for (int i = 0; i < nn; ++i)
+  {
+    // FIXME again what's the purpose of this test ?
+    // in this algo low==0 and high==nn-1 !!
+    if (i < low || i > high)
+    {
+      m_eivec.row(i).end(nn-i) = matH.row(i).end(nn-i);
+    }
+  }
+
+  // Back transformation to get eigenvectors of original matrix
+  int bRows = high-low+1;
+  for (int j = nn-1; j >= low; j--)
+  {
+    int bSize = std::min(j,high)-low+1;
+    m_eivec.col(j).segment(low, bRows) = (m_eivec.block(low, low, bRows, bSize) * matH.col(j).segment(low, bSize));
+  }
+}
+
+#endif // EIGEN_EIGENSOLVER_H

Eigen/src/Eigenvalues/HessenbergDecomposition.h

@@ -0,0 +1,206 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_HESSENBERGDECOMPOSITION_H
+#define EIGEN_HESSENBERGDECOMPOSITION_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+  * \nonstableyet
+  *
+  * \class HessenbergDecomposition
+  *
+  * \brief Reduces a squared matrix to an Hessemberg form
+  *
+  * \param 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$ such that:
+  * \f$ A = Q H Q^* \f$ where \f$ Q \f$ is unitary and \f$ H \f$ a Hessenberg matrix.
+  *
+  * \sa class Tridiagonalization, class Qr
+  */
+template<typename _MatrixType> class HessenbergDecomposition
+{
+  public:
+
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+
+    enum {
+      Size = MatrixType::RowsAtCompileTime,
+      SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic
+                   ? Dynamic
+                   : MatrixType::RowsAtCompileTime-1
+    };
+
+    typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType;
+    typedef Matrix<RealScalar, Size, 1> DiagonalType;
+    typedef Matrix<RealScalar, SizeMinusOne, 1> SubDiagonalType;
+
+    typedef typename NestByValue<Diagonal<MatrixType,0> >::RealReturnType DiagonalReturnType;
+
+    typedef typename NestByValue<Diagonal<
+        NestByValue<Block<MatrixType,SizeMinusOne,SizeMinusOne> >,0 > >::RealReturnType SubDiagonalReturnType;
+
+    /** This constructor initializes a HessenbergDecomposition object for
+      * further use with HessenbergDecomposition::compute()
+      */
+    HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size)
+      : m_matrix(size,size), m_hCoeffs(size-1)
+    {}
+
+    HessenbergDecomposition(const MatrixType& matrix)
+      : m_matrix(matrix),
+        m_hCoeffs(matrix.cols()-1)
+    {
+      _compute(m_matrix, m_hCoeffs);
+    }
+
+    /** Computes or re-compute the Hessenberg decomposition for the matrix \a matrix.
+      *
+      * This method allows to re-use the allocated data.
+      */
+    void compute(const MatrixType& matrix)
+    {
+      m_matrix = matrix;
+      m_hCoeffs.resize(matrix.rows()-1,1);
+      _compute(m_matrix, m_hCoeffs);
+    }
+
+    /** \returns the householder coefficients allowing to
+      * reconstruct the matrix Q from the packed data.
+      *
+      * \sa packedMatrix()
+      */
+    CoeffVectorType householderCoefficients() const { return m_hCoeffs; }
+
+    /** \returns the internal result of the decomposition.
+      *
+      * 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 follow:
+      *       Q = H_{N-1} ... H_1 H_0
+      *    where the matrices H are the Householder transformation:
+      *       H_i = (I - h_i * v_i * v_i')
+      *    where h_i == householderCoefficients()[i] and v_i is a Householder vector:
+      *       v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ]
+      *
+      * See LAPACK for further details on this packed storage.
+      */
+    const MatrixType& packedMatrix(void) const { return m_matrix; }
+
+    MatrixType matrixQ() const;
+    MatrixType matrixH() const;
+
+  private:
+
+    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
+
+  protected:
+    MatrixType m_matrix;
+    CoeffVectorType m_hCoeffs;
+};
+
+#ifndef EIGEN_HIDE_HEAVY_CODE
+
+/** \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)
+{
+  assert(matA.rows()==matA.cols());
+  int n = matA.rows();
+  Matrix<Scalar,1,Dynamic> temp(n);
+  for (int i = 0; i<n-1; ++i)
+  {
+    // let's consider the vector v = i-th column starting at position i+1
+    int remainingSize = n-i-1;
+    RealScalar beta;
+    Scalar h;
+    matA.col(i).end(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.corner(BottomRight, remainingSize, remainingSize)
+        .applyHouseholderOnTheLeft(matA.col(i).end(remainingSize-1), h, &temp.coeffRef(0));
+
+    // A = A H'
+    matA.corner(BottomRight, n, remainingSize)
+        .applyHouseholderOnTheRight(matA.col(i).end(remainingSize-1).conjugate(), ei_conj(h), &temp.coeffRef(0));
+  }
+}
+
+/** reconstructs and returns the matrix Q */
+template<typename MatrixType>
+typename HessenbergDecomposition<MatrixType>::MatrixType
+HessenbergDecomposition<MatrixType>::matrixQ() const
+{
+  int n = m_matrix.rows();
+  MatrixType matQ = MatrixType::Identity(n,n);
+  Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(n);
+  for (int i = n-2; i>=0; i--)
+  {
+    matQ.corner(BottomRight,n-i-1,n-i-1)
+        .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0));
+  }
+  return matQ;
+}
+
+#endif // EIGEN_HIDE_HEAVY_CODE
+
+/** constructs and returns the matrix H.
+  * Note that the matrix H is equivalent to the upper part of the packed matrix
+  * (including the lower sub-diagonal). Therefore, it might be often sufficient
+  * to directly use the packed matrix instead of creating a new one.
+  */
+template<typename MatrixType>
+typename HessenbergDecomposition<MatrixType>::MatrixType
+HessenbergDecomposition<MatrixType>::matrixH() const
+{
+  // FIXME should this function (and other similar) rather take a matrix as argument
+  // and fill it (to avoid temporaries)
+  int n = m_matrix.rows();
+  MatrixType matH = m_matrix;
+  if (n>2)
+    matH.corner(BottomLeft,n-2, n-2).template triangularView<LowerTriangular>().setZero();
+  return matH;
+}
+
+#endif // EIGEN_HESSENBERGDECOMPOSITION_H

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -0,0 +1,366 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
+#define EIGEN_SELFADJOINTEIGENSOLVER_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+  * \nonstableyet
+  *
+  * \class SelfAdjointEigenSolver
+  *
+  * \brief Eigen values/vectors solver for selfadjoint matrix
+  *
+  * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+  *
+  * \note MatrixType must be an actual Matrix type, it can't be an expression type.
+  *
+  * \sa MatrixBase::eigenvalues(), class EigenSolver
+  */
+template<typename _MatrixType> class SelfAdjointEigenSolver
+{
+  public:
+
+    enum {Size = _MatrixType::RowsAtCompileTime };
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+    typedef std::complex<RealScalar> Complex;
+    typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
+    typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
+    typedef Tridiagonalization<MatrixType> TridiagonalizationType;
+//     typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType;
+
+    SelfAdjointEigenSolver()
+        : m_eivec(int(Size), int(Size)),
+          m_eivalues(int(Size))
+    {
+      ei_assert(Size!=Dynamic);
+    }
+
+    SelfAdjointEigenSolver(int size)
+        : m_eivec(size, size),
+          m_eivalues(size)
+    {}
+
+    /** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix,
+      * as well as the eigenvectors if \a computeEigenvectors is true.
+      *
+      * \sa compute(MatrixType,bool), SelfAdjointEigenSolver(MatrixType,MatrixType,bool)
+      */
+    SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
+      : m_eivec(matrix.rows(), matrix.cols()),
+        m_eivalues(matrix.cols())
+    {
+      compute(matrix, computeEigenvectors);
+    }
+
+    /** Constructors computing the eigenvalues of the generalized eigen problem
+      * \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$ . The eigenvectors
+      * are computed if \a computeEigenvectors is true.
+      *
+      * \sa compute(MatrixType,MatrixType,bool), SelfAdjointEigenSolver(MatrixType,bool)
+      */
+    SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
+      : m_eivec(matA.rows(), matA.cols()),
+        m_eivalues(matA.cols())
+    {
+      compute(matA, matB, computeEigenvectors);
+    }
+
+    SelfAdjointEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+
+    SelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true);
+
+    /** \returns the computed eigen vectors as a matrix of column vectors */
+    MatrixType eigenvectors(void) const
+    {
+      #ifndef NDEBUG
+      ei_assert(m_eigenvectorsOk);
+      #endif
+      return m_eivec;
+    }
+
+    /** \returns the computed eigen values */
+    RealVectorType eigenvalues(void) const { return m_eivalues; }
+
+    /** \returns the positive square root of the matrix
+      *
+      * \note the matrix itself must be positive in order for this to make sense.
+      */
+    MatrixType operatorSqrt() const
+    {
+      return m_eivec * m_eivalues.cwise().sqrt().asDiagonal() * m_eivec.adjoint();
+    }
+
+    /** \returns the positive inverse square root of the matrix
+      *
+      * \note the matrix itself must be positive definite in order for this to make sense.
+      */
+    MatrixType operatorInverseSqrt() const
+    {
+      return m_eivec * m_eivalues.cwise().inverse().cwise().sqrt().asDiagonal() * m_eivec.adjoint();
+    }
+
+
+  protected:
+    MatrixType m_eivec;
+    RealVectorType m_eivalues;
+    #ifndef NDEBUG
+    bool m_eigenvectorsOk;
+    #endif
+};
+
+#ifndef EIGEN_HIDE_HEAVY_CODE
+
+/** \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 matA the input selfadjoint matrix
+  * \param hCoeffs returned Householder coefficients
+  *
+  * 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<typename RealScalar, typename Scalar>
+static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n);
+
+/** Computes the eigenvalues of the selfadjoint matrix \a matrix,
+  * as well as the eigenvectors if \a computeEigenvectors is true.
+  *
+  * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool)
+  */
+template<typename MatrixType>
+SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+{
+  #ifndef NDEBUG
+  m_eigenvectorsOk = computeEigenvectors;
+  #endif
+  assert(matrix.cols() == matrix.rows());
+  int n = matrix.cols();
+  m_eivalues.resize(n,1);
+
+  if(n==1)
+  {
+    m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
+    m_eivec.setOnes();
+    return *this;
+  }
+
+  m_eivec = matrix;
+
+  // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ?
+  // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times...
+  // (same for diag and subdiag)
+  RealVectorType& diag = m_eivalues;
+  typename TridiagonalizationType::SubDiagonalType subdiag(n-1);
+  TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, computeEigenvectors);
+
+  int end = n-1;
+  int start = 0;
+  while (end>0)
+  {
+    for (int i = start; i<end; ++i)
+      if (ei_isMuchSmallerThan(ei_abs(subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1]))))
+        subdiag[i] = 0;
+
+    // find the largest unreduced block
+    while (end>0 && subdiag[end-1]==0)
+      end--;
+    if (end<=0)
+      break;
+    start = end - 1;
+    while (start>0 && subdiag[start-1]!=0)
+      start--;
+
+    ei_tridiagonal_qr_step(diag.data(), subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
+  }
+
+  // Sort eigenvalues and corresponding vectors.
+  // TODO make the sort optional ?
+  // TODO use a better sort algorithm !!
+  for (int i = 0; i < n-1; ++i)
+  {
+    int k;
+    m_eivalues.segment(i,n-i).minCoeff(&k);
+    if (k > 0)
+    {
+      std::swap(m_eivalues[i], m_eivalues[k+i]);
+      m_eivec.col(i).swap(m_eivec.col(k+i));
+    }
+  }
+  return *this;
+}
+
+/** Computes the eigenvalues of the generalized eigen problem
+  * \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$ . The eigenvectors
+  * are computed if \a computeEigenvectors is true.
+  *
+  * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool)
+  */
+template<typename MatrixType>
+SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::
+compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors)
+{
+  ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
+
+  // Compute the cholesky decomposition of matB = L L'
+  LLT<MatrixType> cholB(matB);
+
+  // compute C = inv(L) A inv(L')
+  MatrixType matC = matA;
+  cholB.matrixL().solveInPlace(matC);
+  // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' :
+  matC.adjointInPlace();
+  cholB.matrixL().solveInPlace(matC);
+  matC.adjointInPlace();
+  // this version works too:
+//   matC = matC.transpose();
+//   cholB.matrixL().conjugate().template marked<LowerTriangular>().solveTriangularInPlace(matC);
+//   matC = matC.transpose();
+  // FIXME: this should work: (currently it only does for small matrices)
+//   Transpose<MatrixType> trMatC(matC);
+//   cholB.matrixL().conjugate().eval().template marked<LowerTriangular>().solveTriangularInPlace(trMatC);
+
+  compute(matC, computeEigenvectors);
+
+  if (computeEigenvectors)
+  {
+    // transform back the eigen vectors: evecs = inv(U) * evecs
+    cholB.matrixU().solveInPlace(m_eivec);
+    for (int i=0; i<m_eivec.cols(); ++i)
+      m_eivec.col(i) = m_eivec.col(i).normalized();
+  }
+  return *this;
+}
+
+#endif // EIGEN_HIDE_HEAVY_CODE
+
+/** \eigenvalues_module
+  *
+  * \returns a vector listing the eigenvalues of this matrix.
+  */
+template<typename Derived>
+inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1>
+MatrixBase<Derived>::eigenvalues() const
+{
+  ei_assert(Flags&SelfAdjointBit);
+  return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(eval(),false).eigenvalues();
+}
+
+template<typename Derived, bool IsSelfAdjoint>
+struct ei_operatorNorm_selector
+{
+  static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+  operatorNorm(const MatrixBase<Derived>& m)
+  {
+    // 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 m.eigenvalues().cwise().abs().maxCoeff();
+  }
+};
+
+template<typename Derived> struct ei_operatorNorm_selector<Derived, false>
+{
+  static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+  operatorNorm(const MatrixBase<Derived>& m)
+  {
+    typename Derived::PlainMatrixType m_eval(m);
+    // 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 ei_sqrt(
+             (m_eval*m_eval.adjoint())
+             .template marked<SelfAdjoint>()
+             .eigenvalues()
+             .maxCoeff()
+           );
+  }
+};
+
+/** \eigenvalues_module
+  *
+  * \returns the matrix norm of this matrix.
+  */
+template<typename Derived>
+inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+MatrixBase<Derived>::operatorNorm() const
+{
+  return ei_operatorNorm_selector<Derived, Flags&SelfAdjointBit>
+       ::operatorNorm(derived());
+}
+
+#ifndef EIGEN_EXTERN_INSTANTIATIONS
+template<typename RealScalar, typename Scalar>
+static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n)
+{
+  RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
+  RealScalar e2 = ei_abs2(subdiag[end-1]);
+  RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * ei_sqrt(td*td + e2));
+  RealScalar x = diag[start] - mu;
+  RealScalar z = subdiag[start];
+
+  for (int k = start; k < end; ++k)
+  {
+    PlanarRotation<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)
+    {
+      Map<Matrix<Scalar,Dynamic,Dynamic> > q(matrixQ,n,n);
+      q.applyOnTheRight(k,k+1,rot);
+    }
+  }
+}
+#endif
+
+#endif // EIGEN_SELFADJOINTEIGENSOLVER_H

Eigen/src/Eigenvalues/Tridiagonalization.h

@@ -0,0 +1,317 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TRIDIAGONALIZATION_H
+#define EIGEN_TRIDIAGONALIZATION_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+  * \nonstableyet
+  *
+  * \class Tridiagonalization
+  *
+  * \brief Trigiagonal decomposition of a selfadjoint matrix
+  *
+  * \param MatrixType the type of the matrix of which we are performing the tridiagonalization
+  *
+  * 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.
+  *
+  * \sa MatrixBase::tridiagonalize()
+  */
+template<typename _MatrixType> class Tridiagonalization
+{
+  public:
+
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+    typedef typename ei_packet_traits<Scalar>::type Packet;
+
+    enum {
+      Size = MatrixType::RowsAtCompileTime,
+      SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic
+                   ? Dynamic
+                   : MatrixType::RowsAtCompileTime-1,
+      PacketSize = ei_packet_traits<Scalar>::size
+    };
+
+    typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType;
+    typedef Matrix<RealScalar, Size, 1> DiagonalType;
+    typedef Matrix<RealScalar, SizeMinusOne, 1> SubDiagonalType;
+
+    typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+              typename NestByValue<Diagonal<MatrixType,0> >::RealReturnType,
+              Diagonal<MatrixType,0>
+            >::ret DiagonalReturnType;
+
+    typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+              typename NestByValue<Diagonal<
+                NestByValue<Block<MatrixType,SizeMinusOne,SizeMinusOne> >,0 > >::RealReturnType,
+              Diagonal<
+                NestByValue<Block<MatrixType,SizeMinusOne,SizeMinusOne> >,0 >
+            >::ret SubDiagonalReturnType;
+
+    /** This constructor initializes a Tridiagonalization object for
+      * further use with Tridiagonalization::compute()
+      */
+    Tridiagonalization(int size = Size==Dynamic ? 2 : Size)
+      : m_matrix(size,size), m_hCoeffs(size-1)
+    {}
+
+    Tridiagonalization(const MatrixType& matrix)
+      : m_matrix(matrix), m_hCoeffs(matrix.cols()-1)
+    {
+      _compute(m_matrix, m_hCoeffs);
+    }
+
+    /** Computes or re-compute the tridiagonalization for the matrix \a matrix.
+      *
+      * This method allows to re-use the allocated data.
+      */
+    void compute(const MatrixType& matrix)
+    {
+      m_matrix = matrix;
+      m_hCoeffs.resize(matrix.rows()-1, 1);
+      _compute(m_matrix, m_hCoeffs);
+    }
+
+    /** \returns the householder coefficients allowing to
+      * reconstruct the matrix Q from the packed data.
+      *
+      * \sa packedMatrix()
+      */
+    inline CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; }
+
+    /** \returns the internal result of the decomposition.
+      *
+      * The returned matrix contains the following information:
+      *  - the strict upper part is equal to the input matrix A
+      *  - the diagonal and lower sub-diagonal represent the tridiagonal symmetric matrix (real).
+      *  - 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 follow:
+      *       Q = H_{N-1} ... H_1 H_0
+      *    where the matrices H are the Householder transformations:
+      *       H_i = (I - h_i * v_i * v_i')
+      *    where h_i == householderCoefficients()[i] and v_i is a Householder vector:
+      *       v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ]
+      *
+      * See LAPACK for further details on this packed storage.
+      */
+    inline const MatrixType& packedMatrix(void) const { return m_matrix; }
+
+    MatrixType matrixQ() const;
+    template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const;
+    MatrixType matrixT() const;
+    const DiagonalReturnType diagonal(void) const;
+    const SubDiagonalReturnType subDiagonal(void) const;
+
+    static void decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true);
+
+    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
+
+  protected:
+
+    static void _decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true);
+
+    MatrixType m_matrix;
+    CoeffVectorType m_hCoeffs;
+};
+
+/** \returns an expression of the diagonal vector */
+template<typename MatrixType>
+const typename Tridiagonalization<MatrixType>::DiagonalReturnType
+Tridiagonalization<MatrixType>::diagonal(void) const
+{
+  return m_matrix.diagonal().nestByValue();
+}
+
+/** \returns an expression of the sub-diagonal vector */
+template<typename MatrixType>
+const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
+Tridiagonalization<MatrixType>::subDiagonal(void) const
+{
+  int n = m_matrix.rows();
+  return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1)
+    .nestByValue().diagonal().nestByValue();
+}
+
+/** constructs and returns the tridiagonal matrix T.
+  * Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix.
+  * Therefore, it might be often sufficient to directly use the packed matrix, or the vector
+  * expressions returned by diagonal() and subDiagonal() instead of creating a new matrix.
+  */
+template<typename MatrixType>
+typename Tridiagonalization<MatrixType>::MatrixType
+Tridiagonalization<MatrixType>::matrixT(void) const
+{
+  // FIXME should this function (and other similar ones) rather take a matrix as argument
+  // and fill it ? (to avoid temporaries)
+  int n = m_matrix.rows();
+  MatrixType matT = m_matrix;
+  matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate();
+  if (n>2)
+  {
+    matT.corner(TopRight,n-2, n-2).template triangularView<UpperTriangular>().setZero();
+    matT.corner(BottomLeft,n-2, n-2).template triangularView<LowerTriangular>().setZero();
+  }
+  return matT;
+}
+
+#ifndef EIGEN_HIDE_HEAVY_CODE
+
+/** \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 Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs)
+{
+  assert(matA.rows()==matA.cols());
+  int n = matA.rows();
+  Matrix<Scalar,1,Dynamic> aux(n);
+  for (int i = 0; i<n-1; ++i)
+  {
+    int remainingSize = n-i-1;
+    RealScalar beta;
+    Scalar h;
+    matA.col(i).end(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).end(n-i-1)
+    matA.col(i).coeffRef(i+1) = 1;
+
+    hCoeffs.end(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<LowerTriangular>()
+                        * (ei_conj(h) * matA.col(i).end(remainingSize)));
+
+    hCoeffs.end(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.end(remainingSize).dot(matA.col(i).end(remainingSize)))) * matA.col(i).end(n-i-1);
+
+    matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView<LowerTriangular>()
+      .rankUpdate(matA.col(i).end(remainingSize), hCoeffs.end(remainingSize), -1);
+
+    matA.col(i).coeffRef(i+1) = beta;
+    hCoeffs.coeffRef(i) = h;
+  }
+}
+
+/** reconstructs and returns the matrix Q */
+template<typename MatrixType>
+typename Tridiagonalization<MatrixType>::MatrixType
+Tridiagonalization<MatrixType>::matrixQ(void) const
+{
+  MatrixType matQ;
+  matrixQInPlace(&matQ);
+  return matQ;
+}
+
+template<typename MatrixType>
+template<typename QDerived>
+void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) const
+{
+  QDerived& matQ = q->derived();
+  int n = m_matrix.rows();
+  matQ = MatrixType::Identity(n,n);
+  Matrix<Scalar,1,Dynamic> aux(n);
+  for (int i = n-2; i>=0; i--)
+  {
+    matQ.corner(BottomRight,n-i-1,n-i-1)
+        .applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0));
+  }
+}
+
+/** Performs a full decomposition in place */
+template<typename MatrixType>
+void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
+{
+  int n = mat.rows();
+  ei_assert(mat.cols()==n && diag.size()==n && subdiag.size()==n-1);
+  if (n==3 && (!NumTraits<Scalar>::IsComplex) )
+  {
+    _decomposeInPlace3x3(mat, diag, subdiag, extractQ);
+  }
+  else
+  {
+    Tridiagonalization tridiag(mat);
+    diag = tridiag.diagonal();
+    subdiag = tridiag.subDiagonal();
+    if (extractQ)
+      tridiag.matrixQInPlace(&mat);
+  }
+}
+
+/** \internal
+  * Optimized path for 3x3 matrices.
+  * Especially useful for plane fitting.
+  */
+template<typename MatrixType>
+void Tridiagonalization<MatrixType>::_decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
+{
+  diag[0] = ei_real(mat(0,0));
+  RealScalar v1norm2 = ei_abs2(mat(0,2));
+  if (ei_isMuchSmallerThan(v1norm2, RealScalar(1)))
+  {
+    diag[1] = ei_real(mat(1,1));
+    diag[2] = ei_real(mat(2,2));
+    subdiag[0] = ei_real(mat(0,1));
+    subdiag[1] = ei_real(mat(1,2));
+    if (extractQ)
+      mat.setIdentity();
+  }
+  else
+  {
+    RealScalar beta = ei_sqrt(ei_abs2(mat(0,1))+v1norm2);
+    RealScalar invBeta = RealScalar(1)/beta;
+    Scalar m01 = mat(0,1) * invBeta;
+    Scalar m02 = mat(0,2) * invBeta;
+    Scalar q = RealScalar(2)*m01*mat(1,2) + m02*(mat(2,2) - mat(1,1));
+    diag[1] = ei_real(mat(1,1) + m02*q);
+    diag[2] = ei_real(mat(2,2) - m02*q);
+    subdiag[0] = beta;
+    subdiag[1] = ei_real(mat(1,2) - m01 * q);
+    if (extractQ)
+    {
+      mat(0,0) = 1;
+      mat(0,1) = 0;
+      mat(0,2) = 0;
+      mat(1,0) = 0;
+      mat(1,1) = m01;
+      mat(1,2) = m02;
+      mat(2,0) = 0;
+      mat(2,1) = m02;
+      mat(2,2) = -m01;
+    }
+  }
+}
+
+#endif // EIGEN_HIDE_HEAVY_CODE
+
+#endif // EIGEN_TRIDIAGONALIZATION_H