* added a Tridiagonalization class for selfadjoint matrices

* added MatrixBase::real()
* added the ability to extract a selfadjoint matrix from the
  lower or upper part of a matrix, e.g.:
    m.extract<Upper|SelfAdjoint>()
  will ignore the strict lower part and return a selfadjoint.
  This is compatible with ZeroDiag and UnitDiag.

7 files changed, 268 insertions, 4 deletions

diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h
index 893c8abf6..8715f931e 100644
--- a/Eigen/src/Core/CwiseUnaryOp.h
+++ b/Eigen/src/Core/CwiseUnaryOp.h
@@ -139,7 +139,7 @@ MatrixBase<Derived>::cwiseAbs2() const
   return derived();
 }
 
-/** \returns an expression of the complex conjugate of *this.
+/** \returns an expression of the complex conjugate of \c *this.
   *
   * \sa adjoint() */
 template<typename Derived>
@@ -149,6 +149,16 @@ MatrixBase<Derived>::conjugate() const
   return ConjugateReturnType(derived());
 }
 
+/** \returns an expression of the real part of \c *this.
+  *
+  * \sa adjoint() */
+template<typename Derived>
+inline const typename MatrixBase<Derived>::RealReturnType
+MatrixBase<Derived>::real() const
+{
+  return derived();
+}
+
 /** \returns an expression of *this with the \a Scalar type casted to
   * \a NewScalar.
   *
diff --git a/Eigen/src/Core/Extract.h b/Eigen/src/Core/Extract.h
index ea7240ada..6671e45b2 100755
--- a/Eigen/src/Core/Extract.h
+++ b/Eigen/src/Core/Extract.h
@@ -77,7 +77,7 @@ template<typename MatrixType, unsigned int Mode> class Extract
     inline Scalar _coeff(int row, int col) const
     {
       if(Flags & LowerTriangularBit ? col>row : row>col)
-        return (Scalar)0;
+        return (Flags & SelfAdjointBit) ? ei_conj(m_matrix.coeff(col, row)) : (Scalar)0;
       if(Flags & UnitDiagBit)
         return col==row ? (Scalar)1 : m_matrix.coeff(row, col);
       else if(Flags & ZeroDiagBit)
diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h
index 6a4a1213d..f3ffe10aa 100644
--- a/Eigen/src/Core/Functors.h
+++ b/Eigen/src/Core/Functors.h
@@ -204,6 +204,19 @@ template<typename Scalar, typename NewType>
 struct ei_functor_traits<ei_scalar_cast_op<Scalar,NewType> >
 { enum { Cost = ei_is_same_type<Scalar, NewType>::ret ? 0 : NumTraits<NewType>::AddCost, IsVectorizable = false }; };
 
+/** \internal
+  * \brief Template functor to extract the real part of a complex
+  *
+  * \sa class CwiseUnaryOp, MatrixBase::real()
+  */
+template<typename Scalar>
+struct ei_scalar_real_op EIGEN_EMPTY_STRUCT {
+  typedef typename NumTraits<Scalar>::Real result_type;
+  inline result_type operator() (const Scalar& a) const { return ei_real(a); }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_real_op<Scalar> >
+{ enum { Cost =  0, IsVectorizable = false }; };
 
 /** \internal
   * \brief Template functor to multiply a scalar by a fixed other one
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index ab35b94a6..a29504b33 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -197,6 +197,8 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived>
                         CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Derived>,
                         Derived&
                      >::ret ConjugateReturnType;
+    /** the return type of MatrixBase::real() */
+    typedef CwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived> RealReturnType;
     /** the return type of MatrixBase::adjoint() */
     typedef Transpose<NestByValue<typename ei_unref<ConjugateReturnType>::type> >
             AdjointReturnType;
@@ -477,6 +479,7 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived>
     /// \name Coefficient-wise operations
     //@{
     const ConjugateReturnType conjugate() const;
+    const RealReturnType real() const;
 
     template<typename OtherDerived>
     const CwiseBinaryOp<ei_scalar_product_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index e17563c9a..fd3451922 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -59,8 +59,6 @@ const unsigned int Upper = UpperTriangularBit;
 const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit;
 const unsigned int Lower = LowerTriangularBit;
 const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit;
-
-// additional possible values for the Mode parameter of part()
 const unsigned int SelfAdjoint = SelfAdjointBit;
 
 // additional possible values for the Mode parameter of extract()
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index d48d4c325..a9974c38b 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -59,6 +59,7 @@ template<typename Scalar> struct ei_scalar_product_op;
 template<typename Scalar> struct ei_scalar_quotient_op;
 template<typename Scalar> struct ei_scalar_opposite_op;
 template<typename Scalar> struct ei_scalar_conjugate_op;
+template<typename Scalar> struct ei_scalar_real_op;
 template<typename Scalar> struct ei_scalar_abs_op;
 template<typename Scalar> struct ei_scalar_abs2_op;
 template<typename Scalar> struct ei_scalar_sqrt_op;
diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h
new file mode 100755
index 000000000..a85068e38
--- /dev/null
+++ b/Eigen/src/QR/Tridiagonalization.h
@@ -0,0 +1,239 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// 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
+
+/** \class Tridiagonalization
+  *
+  * \brief Trigiagonal decomposition of a selfadjoint matrix
+  *
+  * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+  *
+  * 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 unitatry 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;
+
+    enum {SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic
+                        ? Dynamic
+                        : MatrixType::RowsAtCompileTime-1};
+
+    typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType;
+
+    typedef typename NestByValue<DiagonalCoeffs<MatrixType> >::RealReturnType DiagonalType;
+
+    typedef typename NestByValue<DiagonalCoeffs<
+        NestByValue<Block<
+          MatrixType,SizeMinusOne,SizeMinusOne> > > >::RealReturnType SubDiagonalType;
+
+    Tridiagonalization()
+    {}
+
+    Tridiagonalization(int rows, int cols)
+      : m_matrix(rows,cols), m_hCoeffs(rows-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);
+      _compute(m_matrix, m_hCoeffs);
+    }
+
+    /** \returns the householder coefficients allowing to
+      * reconstruct the matrix Q from the packed data.
+      *
+      * \sa packedMatrix()
+      */
+    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 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(void) const;
+    const DiagonalType diagonal(void) const;
+    const SubDiagonalType subDiagonal(void) const;
+
+  private:
+
+    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
+
+  protected:
+    MatrixType m_matrix;
+    CoeffVectorType m_hCoeffs;
+};
+
+
+/** \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:
+  *
+  * \sa packedMatrix()
+  */
+template<typename MatrixType>
+void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs)
+{
+  assert(matA.rows()==matA.cols());
+  int n = matA.rows();
+  for (int i = 0; i<n-2; ++i)
+  {
+    // let's consider the vector v = i-th column starting at position i+1
+
+    // start of the householder transformation
+    // squared norm of the vector v skipping the first element
+    RealScalar v1norm2 = matA.col(i).end(n-(i+2)).norm2();
+
+    if (ei_isMuchSmallerThan(v1norm2,static_cast<Scalar>(1)))
+    {
+      hCoeffs.coeffRef(i) = 0.;
+    }
+    else
+    {
+      Scalar v0 = matA.col(i).coeff(i+1);
+      RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2);
+      if (ei_real(v0)>=0.)
+        beta = -beta;
+      matA.col(i).end(n-(i+2)) *= (1./(v0-beta));
+      matA.col(i).coeffRef(i+1) = beta;
+      Scalar h = (beta - v0) / beta;
+      // end of the householder transformation
+
+      // 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;
+      // let's use the end of hCoeffs to store temporary values
+      hCoeffs.end(n-i-1) = h * (matA.corner(BottomRight,n-i-1,n-i-1).template extract<Lower|SelfAdjoint>()
+                                * matA.col(i).end(n-i-1));
+
+
+      hCoeffs.end(n-i-1) += (h * (-0.5) * matA.col(i).end(n-i-1).dot(hCoeffs.end(n-i-1)))
+                            * matA.col(i).end(n-i-1);
+
+      matA.corner(BottomRight,n-i-1,n-i-1).template part<Lower>() =
+        matA.corner(BottomRight,n-i-1,n-i-1) - (
+            (matA.col(i).end(n-i-1) * hCoeffs.end(n-i-1).adjoint()).lazy()
+          + (hCoeffs.end(n-i-1) * matA.col(i).end(n-i-1).adjoint()).lazy() );
+      // FIXME check that the above expression does follow the lazy path (no temporary and
+      // only lower products are evaluated)
+      // FIXME can we avoid to evaluate twice the diagonal products ?
+      // (in a simple way otherwise it's overkill)
+
+      matA.col(i).coeffRef(i+1) = beta;
+
+      hCoeffs.coeffRef(i) = h;
+    }
+  }
+  if (NumTraits<Scalar>::IsComplex)
+  {
+    // householder transformation on the remaining single scalar
+    int i = n-2;
+    Scalar v0 = matA.col(i).coeff(i+1);
+    RealScalar beta = ei_abs(v0);
+    if (ei_real(v0)>=0.)
+      beta = -beta;
+    matA.col(i).coeffRef(i+1) = beta;
+    hCoeffs.coeffRef(i) = (beta - v0) / beta;
+  }
+}
+
+/** reconstructs and returns the matrix Q */
+template<typename MatrixType>
+typename Tridiagonalization<MatrixType>::MatrixType
+Tridiagonalization<MatrixType>::matrixQ(void) const
+{
+  int n = m_matrix.rows();
+  MatrixType matQ = MatrixType::identity(n,n);
+  for (int i = n-2; i>=0; i--)
+  {
+    Scalar tmp = m_matrix.coeff(i+1,i);
+    m_matrix.const_cast_derived().coeffRef(i+1,i) = 1;
+
+    matQ.corner(BottomRight,n-i-1,n-i-1) -=
+      ((m_hCoeffs[i] * m_matrix.col(i).end(n-i-1)) *
+      (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy();
+
+    m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp;
+  }
+  return matQ;
+}
+
+/** \returns an expression of the diagonal vector */
+template<typename MatrixType>
+const typename Tridiagonalization<MatrixType>::DiagonalType
+Tridiagonalization<MatrixType>::diagonal(void) const
+{
+  return m_matrix.diagonal().nestByValue().real();
+}
+
+/** \returns an expression of the sub-diagonal vector */
+template<typename MatrixType>
+const typename Tridiagonalization<MatrixType>::SubDiagonalType
+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().real();
+}
+
+#endif // EIGEN_TRIDIAGONALIZATION_H