* add .imag() function

* fix a very old bug in EigenSolver that I had completely forgotten
  (thanks to Timothy to refresh my mind)
* fix doc of Matrix::Map

Eigen/src/Core/CwiseUnaryOp.h

@@ -174,13 +174,17 @@ MatrixBase<Derived>::conjugate() const
 
 /** \returns an expression of the real part of \c *this.
   *
-  * \sa adjoint() */
+  * \sa imag() */
 template<typename Derived>
 inline const typename MatrixBase<Derived>::RealReturnType
-MatrixBase<Derived>::real() const
+MatrixBase<Derived>::real() const { return derived(); }
-{
+
-  return derived();
+/** \returns an expression of the imaginary part of \c *this.
-}
+  *
+  * \sa real() */
+template<typename Derived>
+inline const typename MatrixBase<Derived>::ImagReturnType
+MatrixBase<Derived>::imag() const { return derived(); }
 
 /** \returns an expression of *this with the \a Scalar type casted to
   * \a NewScalar.

Eigen/src/Core/Functors.h

@@ -240,7 +240,21 @@ struct ei_scalar_real_op EIGEN_EMPTY_STRUCT {
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_real_op<Scalar> >
-{ enum { Cost =  0, PacketAccess = false }; };
+{ enum { Cost = 0, PacketAccess = false }; };
+
+/** \internal
+  * \brief Template functor to extract the imaginary part of a complex
+  *
+  * \sa class CwiseUnaryOp, MatrixBase::imag()
+  */
+template<typename Scalar>
+struct ei_scalar_imag_op EIGEN_EMPTY_STRUCT {
+  typedef typename NumTraits<Scalar>::Real result_type;
+  inline result_type operator() (const Scalar& a) const { return ei_imag(a); }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_imag_op<Scalar> >
+{ enum { Cost = 0, PacketAccess = false }; };
 
 /** \internal
   * \brief Template functor to multiply a scalar by a fixed other one

Eigen/src/Core/Matrix.h

@@ -437,14 +437,14 @@ class Matrix
         this->Base::swap(other);
     }
 
-    /**
+    /** \name Map
       * These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
       * while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned
       * \a data pointers.
-      * 
+      *
       * \see class Map
       */
-    //@}
+    //@{
     inline static const UnalignedMapType Map(const Scalar* data)
     { return UnalignedMapType(data); }
     inline static UnalignedMapType Map(Scalar* data)

Eigen/src/Core/MatrixBase.h

@@ -189,6 +189,8 @@ template<typename Derived> class MatrixBase
                      >::ret ConjugateReturnType;
     /** \internal the return type of MatrixBase::real() */
     typedef CwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived> RealReturnType;
+    /** \internal the return type of MatrixBase::imag() */
+    typedef CwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
     /** \internal the return type of MatrixBase::adjoint() */
     typedef Eigen::Transpose<NestByValue<typename ei_cleantype<ConjugateReturnType>::type> >
             AdjointReturnType;
@@ -496,6 +498,7 @@ template<typename Derived> class MatrixBase
 
     ConjugateReturnType conjugate() const;
     const RealReturnType real() const;
+    const ImagReturnType imag() const;
 
     template<typename CustomUnaryOp>
     const CwiseUnaryOp<CustomUnaryOp, Derived> unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const;

Eigen/src/Core/util/ForwardDeclarations.h

@@ -64,6 +64,7 @@ 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_imag_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;

Eigen/src/QR/EigenSolver.h

@@ -59,7 +59,7 @@ template<typename _MatrixType> class EigenSolver
       compute(matrix);
     }
 
-    
+
     EigenvectorType eigenvectors(void) const;
 
     /** \returns a real matrix V of pseudo eigenvectors.
@@ -200,18 +200,18 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
   for (int m = low+1; m <= high-1; m++)
   {
     // Scale column.
-    Scalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum();
+    RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum();
     if (scale != 0.0)
     {
       // Compute Householder transformation.
-      Scalar h = 0.0;
+      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);
       }
-      Scalar g = ei_sqrt(h);
+      RealScalar g = ei_sqrt(h);
       if (ort.coeff(m) > 0)
         g = -g;
       h = h - ort.coeff(m) * g;
@@ -247,7 +247,6 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
   }
 }
 
-
 // Complex scalar division.
 template<typename Scalar>
 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
@@ -298,7 +297,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
       m_eivalues.coeffRef(j).real() = matH.coeff(j,j);
       m_eivalues.coeffRef(j).imag() = 0.0;
     }
-    norm += matH.col(j).start(std::min(j+1,nn)).cwise().abs().sum();
+    norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwise().abs().sum();
   }
 
   // Outer loop over eigenvalue index
@@ -571,7 +570,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
       for (int i = n-1; i >= 0; i--)
       {
         w = matH.coeff(i,i) - p;
-        r = (matH.row(i).end(nn-l) * matH.col(n).end(nn-l))(0,0);
+        r = (matH.row(i).segment(l,n-l+1) * matH.col(n).segment(l, n-l+1))(0,0);
 
         if (m_eivalues.coeff(i).imag() < 0.0)
         {
@@ -630,8 +629,8 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
       for (int i = n-2; i >= 0; i--)
       {
         Scalar ra,sa,vr,vi;
-        ra = (matH.row(i).end(nn-l) * matH.col(n-1).end(nn-l)).lazy()(0,0);
+        ra = (matH.block(i,l, 1, n-l+1) * matH.block(l,n-1, n-l+1, 1)).lazy()(0,0);
-        sa = (matH.row(i).end(nn-l) * matH.col(n).end(nn-l)).lazy()(0,0);
+        sa = (matH.block(i,l, 1, n-l+1) * matH.block(l,n, n-l+1, 1)).lazy()(0,0);
         w = matH.coeff(i,i) - p;
 
         if (m_eivalues.coeff(i).imag() < 0.0)

test/eigensolver.cpp

@@ -130,27 +130,13 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   MatrixType a1 = MatrixType::Random(rows,cols);
   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
 
-//   MatrixType b = MatrixType::Random(rows,cols);
-//   MatrixType b1 = MatrixType::Random(rows,cols);
-//   MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
-
   EigenSolver<MatrixType> ei0(symmA);
+  VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix());
   VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
     (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
 
-//   a = a /*+ symmA*/;
   EigenSolver<MatrixType> ei1(a);
-
+  VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
-  IOFormat OctaveFmt(4, AlignCols, ", ", ";\n", "", "", "[", "]");
-//   std::cerr << "==============\n" << a.format(OctaveFmt) << "\n\n" << ei1.eigenvalues().transpose() << "\n\n";
-//   std::cerr << a * ei1.pseudoEigenvectors() << "\n\n" << ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix() << "\n\n\n";
-
-//   VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
-
-  
-//   std::cerr << a.format(OctaveFmt) << "\n\n";
-//   std::cerr << ei1.eigenvectors().format(OctaveFmt) << "\n\n";
-//   std::cerr << a.template cast<Complex>() * ei1.eigenvectors() << "\n\n" << ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval() << "\n\n";
   VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
                    ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval());
 
@@ -167,8 +153,7 @@ void test_eigensolver()
     CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) );
 
     CALL_SUBTEST( eigensolver(Matrix4f()) );
-    // FIXME the test fails for larger matrices
+    CALL_SUBTEST( eigensolver(MatrixXd(17,17)) );
-//     CALL_SUBTEST( eigensolver(MatrixXd(7,7)) );
   }
 }