merge

Eigen/src/Array/VectorwiseOp.h

@@ -355,10 +355,22 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
 
 /////////// Artithmetic operators ///////////
 
+    /** Copies the vector \a other to each subvector of \c *this */
+    template<typename OtherDerived>
+    ExpressionType& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
+      //ei_assert((m_matrix.isNull()) == (other.isNull())); FIXME
+      for(int j=0; j<subVectors(); ++j)
+        subVector(j) = other;
+      return const_cast<ExpressionType&>(m_matrix);
+    }
+
     /** Adds the vector \a other to each subvector of \c *this */
     template<typename OtherDerived>
     ExpressionType& operator+=(const MatrixBase<OtherDerived>& other)
     {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
       for(int j=0; j<subVectors(); ++j)
         subVector(j) += other;
       return const_cast<ExpressionType&>(m_matrix);
@@ -368,6 +380,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
     template<typename OtherDerived>
     ExpressionType& operator-=(const MatrixBase<OtherDerived>& other)
     {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
       for(int j=0; j<subVectors(); ++j)
         subVector(j) -= other;
       return const_cast<ExpressionType&>(m_matrix);

Eigen/src/Core/Assign.h

@@ -412,8 +412,10 @@ template<typename Derived, typename OtherDerived,
          bool EvalBeforeAssigning = (int(OtherDerived::Flags) & EvalBeforeAssigningBit) != 0,
          bool NeedToTranspose = Derived::IsVectorAtCompileTime
                 && OtherDerived::IsVectorAtCompileTime
-                && int(Derived::RowsAtCompileTime) == int(OtherDerived::ColsAtCompileTime)
-                && int(Derived::ColsAtCompileTime) == int(OtherDerived::RowsAtCompileTime)
+                && ((int(Derived::RowsAtCompileTime) == 1 && int(OtherDerived::ColsAtCompileTime) == 1)
+                      |  // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&".
+                         // revert to || as soon as not needed anymore.
+                    (int(Derived::ColsAtCompileTime) == 1 && int(OtherDerived::RowsAtCompileTime) == 1))
                 && int(Derived::SizeAtCompileTime) != 1>
 struct ei_assign_selector;
 

Eigen/src/Core/Fuzzy.h

@@ -49,7 +49,7 @@ template<typename Derived>
 template<typename OtherDerived>
 bool MatrixBase<Derived>::isApprox(
   const MatrixBase<OtherDerived>& other,
-  typename NumTraits<Scalar>::Real prec
+  RealScalar prec
 ) const
 {
   const typename ei_nested<Derived,2>::type nested(derived());
@@ -73,7 +73,7 @@ bool MatrixBase<Derived>::isApprox(
 template<typename Derived>
 bool MatrixBase<Derived>::isMuchSmallerThan(
   const typename NumTraits<Scalar>::Real& other,
-  typename NumTraits<Scalar>::Real prec
+  RealScalar prec
 ) const
 {
   return cwise().abs2().sum() <= prec * prec * other * other;
@@ -93,7 +93,7 @@ template<typename Derived>
 template<typename OtherDerived>
 bool MatrixBase<Derived>::isMuchSmallerThan(
   const MatrixBase<OtherDerived>& other,
-  typename NumTraits<Scalar>::Real prec
+  RealScalar prec
 ) const
 {
   return this->cwise().abs2().sum() <= prec * prec * other.cwise().abs2().sum();
@@ -124,7 +124,7 @@ template<typename Derived>
 template<typename OtherDerived>
 bool MatrixBase<Derived>::isApprox(
   const MatrixBase<OtherDerived>& other,
-  typename NumTraits<Scalar>::Real prec
+  RealScalar prec
 ) const
 {
   return ei_fuzzy_selector<Derived,OtherDerived>::isApprox(derived(), other.derived(), prec);
@@ -143,7 +143,7 @@ bool MatrixBase<Derived>::isApprox(
 template<typename Derived>
 bool MatrixBase<Derived>::isMuchSmallerThan(
   const typename NumTraits<Scalar>::Real& other,
-  typename NumTraits<Scalar>::Real prec
+  RealScalar prec
 ) const
 {
   return ei_fuzzy_selector<Derived>::isMuchSmallerThan(derived(), other, prec);
@@ -163,7 +163,7 @@ template<typename Derived>
 template<typename OtherDerived>
 bool MatrixBase<Derived>::isMuchSmallerThan(
   const MatrixBase<OtherDerived>& other,
-  typename NumTraits<Scalar>::Real prec
+  RealScalar prec
 ) const
 {
   return ei_fuzzy_selector<Derived,OtherDerived>::isMuchSmallerThan(derived(), other.derived(), prec);

Eigen/src/Core/MatrixBase.h

@@ -45,7 +45,7 @@ template<typename Derived> struct AnyMatrixBase
   inline int rows() const { return derived().rows(); }
   /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
   inline int cols() const { return derived().cols(); }
-
+  
   template<typename Dest> inline void evalTo(Dest& dst) const
   { derived().evalTo(dst); }
 
@@ -765,6 +765,19 @@ template<typename Derived> class MatrixBase
     template<typename Derived1, typename Derived2>
     Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product);
 
+////////// Householder module ///////////
+
+    template<typename EssentialPart>
+    void makeHouseholder(EssentialPart *essential,
+                         RealScalar *beta) const;
+    template<typename EssentialPart>
+    void applyHouseholderOnTheLeft(const EssentialPart& essential,
+                                   const RealScalar& beta);
+    template<typename EssentialPart>
+    void applyHouseholderOnTheRight(const EssentialPart& essential,
+                                    const RealScalar& beta);
+    
+
     #ifdef EIGEN_MATRIXBASE_PLUGIN
     #include EIGEN_MATRIXBASE_PLUGIN
     #endif

Eigen/src/Householder/CMakeLists.txt

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

Eigen/src/Householder/Householder.h

@@ -0,0 +1,88 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// 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_HOUSEHOLDER_H
+#define EIGEN_HOUSEHOLDER_H
+
+template<int n> struct ei_decrement_size
+{
+  enum {
+    ret = (n==1 || n==Dynamic) ? n : n-1
+  };
+};
+
+template<typename Derived>
+template<typename EssentialPart>
+void MatrixBase<Derived>::makeHouseholder(
+  EssentialPart *essential,
+  RealScalar *beta) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
+  RealScalar _squaredNorm = squaredNorm();
+  Scalar c0;
+  if(ei_abs2(coeff(0)) <= ei_abs2(precision<Scalar>()) * _squaredNorm)
+  {
+    c0 = ei_sqrt(_squaredNorm);
+  }
+  else
+  {
+    Scalar sign = coeff(0) / ei_abs(coeff(0));
+    c0 = coeff(0) + sign * ei_sqrt(_squaredNorm);
+  }
+  *essential = end(size()-1) / c0; // FIXME take advantage of fixed size
+  const RealScalar c0abs2 = ei_abs2(c0);
+  *beta = RealScalar(2) * c0abs2 / (c0abs2 + _squaredNorm - ei_abs2(coeff(0)));
+}
+
+template<typename Derived>
+template<typename EssentialPart>
+void MatrixBase<Derived>::applyHouseholderOnTheLeft(
+  const EssentialPart& essential,
+  const RealScalar& beta)
+{
+  Matrix<Scalar, 1, ColsAtCompileTime, PlainMatrixType::Options, 1, MaxColsAtCompileTime> tmp(cols());
+  tmp = row(0) + essential.adjoint() * block(1,0,rows()-1,cols());
+  // FIXME take advantage of fixed size
+  // FIXME play with lazy()
+  // FIXME maybe not a good idea to use matrix product
+  row(0) -= beta * tmp;
+  block(1,0,rows()-1,cols()) -= beta * essential * tmp;
+}
+
+template<typename Derived>
+template<typename EssentialPart>
+void MatrixBase<Derived>::applyHouseholderOnTheRight(
+  const EssentialPart& essential,
+  const RealScalar& beta)
+{
+  Matrix<Scalar, RowsAtCompileTime, 1, PlainMatrixType::Options, MaxRowsAtCompileTime, 1> tmp(rows());
+  tmp = col(0) + block(0,1,rows(),cols()-1) * essential.conjugate();
+  // FIXME take advantage of fixed size
+  // FIXME play with lazy()
+  // FIXME maybe not a good idea to use matrix product
+  col(0) -= beta * tmp;
+  block(0,1,rows(),cols()-1) -= beta * tmp * essential.transpose();
+}
+
+#endif // EIGEN_HOUSEHOLDER_H