Merged eigen/eigen into default

6 files changed, 33 insertions, 124 deletions

diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h
index 0769a212e..53f934999 100644
--- a/Eigen/src/Core/GeneralProduct.h
+++ b/Eigen/src/Core/GeneralProduct.h
@@ -213,15 +213,18 @@ template<> struct gemv_dense_selector<OnTheRight,ColMajor,true>
     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(lhs)
                                   * RhsBlasTraits::extractScalarFactor(rhs);
 
+    // make sure Dest is a compile-time vector type (bug 1166)
+    typedef typename conditional<Dest::IsVectorAtCompileTime, Dest, typename Dest::ColXpr>::type ActualDest;
+
     enum {
       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
       // on, the other hand it is good for the cache to pack the vector anyways...
-      EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1,
+      EvalToDestAtCompileTime = (ActualDest::InnerStrideAtCompileTime==1),
       ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
-      MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal
+      MightCannotUseDest = (ActualDest::InnerStrideAtCompileTime!=1) || ComplexByReal
     };
 
-    gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
+    gemv_static_vector_if<ResScalar,ActualDest::SizeAtCompileTime,ActualDest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
 
     const bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
     const bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
@@ -314,7 +317,7 @@ template<> struct gemv_dense_selector<OnTheRight,RowMajor,true>
         actualLhs.rows(), actualLhs.cols(),
         LhsMapper(actualLhs.data(), actualLhs.outerStride()),
         RhsMapper(actualRhsPtr, 1),
-        dest.data(), dest.innerStride(),
+        dest.data(), dest.col(0).innerStride(), //NOTE  if dest is not a vector at compile-time, then dest.innerStride() might be wrong. (bug 1166)
         actualAlpha);
   }
 };
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index e87b60f8f..447f1b834 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -496,7 +496,7 @@ template<typename Scalar, bool IsInteger>
 struct pow_default_impl
 {
   typedef Scalar retval;
-  static inline Scalar run(const Scalar& x, const Scalar& y)
+  static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y)
   {
     EIGEN_USING_STD_MATH(pow);
     return pow(x, y);
@@ -506,7 +506,7 @@ struct pow_default_impl
 template<typename Scalar>
 struct pow_default_impl<Scalar, true>
 {
-  static inline Scalar run(Scalar x, Scalar y)
+  static EIGEN_DEVICE_FUNC inline Scalar run(Scalar x, Scalar y)
   {
     Scalar res(1);
     eigen_assert(!NumTraits<Scalar>::IsSigned || y >= 0);
diff --git a/Eigen/src/Core/SelfCwiseBinaryOp.h b/Eigen/src/Core/SelfCwiseBinaryOp.h
index 38185d9d7..78fff1549 100644
--- a/Eigen/src/Core/SelfCwiseBinaryOp.h
+++ b/Eigen/src/Core/SelfCwiseBinaryOp.h
@@ -13,7 +13,7 @@
 namespace Eigen { 
 
 template<typename Derived>
-inline Derived& DenseBase<Derived>::operator*=(const Scalar& other)
+EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::operator*=(const Scalar& other)
 {
   typedef typename Derived::PlainObject PlainObject;
   internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::mul_assign_op<Scalar>());
@@ -21,7 +21,7 @@ inline Derived& DenseBase<Derived>::operator*=(const Scalar& other)
 }
 
 template<typename Derived>
-inline Derived& ArrayBase<Derived>::operator+=(const Scalar& other)
+EIGEN_STRONG_INLINE Derived& ArrayBase<Derived>::operator+=(const Scalar& other)
 {
   typedef typename Derived::PlainObject PlainObject;
   internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::add_assign_op<Scalar>());
@@ -29,7 +29,7 @@ inline Derived& ArrayBase<Derived>::operator+=(const Scalar& other)
 }
 
 template<typename Derived>
-inline Derived& ArrayBase<Derived>::operator-=(const Scalar& other)
+EIGEN_STRONG_INLINE Derived& ArrayBase<Derived>::operator-=(const Scalar& other)
 {
   typedef typename Derived::PlainObject PlainObject;
   internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::sub_assign_op<Scalar>());
@@ -37,7 +37,7 @@ inline Derived& ArrayBase<Derived>::operator-=(const Scalar& other)
 }
 
 template<typename Derived>
-inline Derived& DenseBase<Derived>::operator/=(const Scalar& other)
+EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::operator/=(const Scalar& other)
 {
   typedef typename Derived::PlainObject PlainObject;
   internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::div_assign_op<Scalar>());
diff --git a/Eigen/src/Geometry/Umeyama.h b/Eigen/src/Geometry/Umeyama.h
index 6943f719e..7e933fca1 100644
--- a/Eigen/src/Geometry/Umeyama.h
+++ b/Eigen/src/Geometry/Umeyama.h
@@ -139,6 +139,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
   if  ( svd.matrixU().determinant() * svd.matrixV().determinant() < 0 )
     S(m-1) = -1;
 
+  // Eq. (40) and (43)
   Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose();
 
   if (with_scaling)
diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h
index 402733cce..b92bb17e2 100644
--- a/Eigen/src/SparseCore/SparseSelfAdjointView.h
+++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h
@@ -387,7 +387,10 @@ void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename Matri
   typedef typename MatrixType::Scalar Scalar;
   typedef SparseMatrix<Scalar,DestOrder,StorageIndex> Dest;
   typedef Matrix<StorageIndex,Dynamic,1> VectorI;
+  typedef evaluator<MatrixType> MatEval;
+  typedef typename evaluator<MatrixType>::InnerIterator MatIterator;
   
+  MatEval matEval(mat);
   Dest& dest(_dest.derived());
   enum {
     StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor)
@@ -401,7 +404,7 @@ void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename Matri
   for(Index j = 0; j<size; ++j)
   {
     Index jp = perm ? perm[j] : j;
-    for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
+    for(MatIterator it(matEval,j); it; ++it)
     {
       Index i = it.index();
       Index r = it.row();
@@ -431,7 +434,7 @@ void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix<typename Matri
   // copy data
   for(StorageIndex j = 0; j<size; ++j)
   {
-    for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
+    for(MatIterator it(matEval,j); it; ++it)
     {
       StorageIndex i = internal::convert_index<StorageIndex>(it.index());
       Index r = it.row();
@@ -474,12 +477,17 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixTyp
   typedef typename MatrixType::Scalar Scalar;
   SparseMatrix<Scalar,DstOrder,StorageIndex>& dest(_dest.derived());
   typedef Matrix<StorageIndex,Dynamic,1> VectorI;
+  typedef evaluator<MatrixType> MatEval;
+  typedef typename evaluator<MatrixType>::InnerIterator MatIterator;
+
   enum {
     SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor,
     StorageOrderMatch = int(SrcOrder) == int(DstOrder),
     DstMode = DstOrder==RowMajor ? (_DstMode==Upper ? Lower : Upper) : _DstMode,
     SrcMode = SrcOrder==RowMajor ? (_SrcMode==Upper ? Lower : Upper) : _SrcMode
   };
+
+  MatEval matEval(mat);
   
   Index size = mat.rows();
   VectorI count(size);
@@ -488,7 +496,7 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixTyp
   for(StorageIndex j = 0; j<size; ++j)
   {
     StorageIndex jp = perm ? perm[j] : j;
-    for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
+    for(MatIterator it(matEval,j); it; ++it)
     {
       StorageIndex i = it.index();
       if((int(SrcMode)==int(Lower) && i<j) || (int(SrcMode)==int(Upper) && i>j))
@@ -508,7 +516,7 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixTyp
   for(StorageIndex j = 0; j<size; ++j)
   {
     
-    for(typename MatrixType::InnerIterator it(mat,j); it; ++it)
+    for(MatIterator it(matEval,j); it; ++it)
     {
       StorageIndex i = it.index();
       if((int(SrcMode)==int(Lower) && i<j) || (int(SrcMode)==int(Upper) && i>j))
diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h
index 7c718e4e1..0c27855d5 100644
--- a/Eigen/src/SparseCore/SparseTriangularView.h
+++ b/Eigen/src/SparseCore/SparseTriangularView.h
@@ -43,9 +43,6 @@ template<typename MatrixType, unsigned int Mode> class TriangularViewImpl<Matrix
     
     EIGEN_SPARSE_PUBLIC_INTERFACE(TriangularViewType)
     
-    class InnerIterator;
-    class ReverseInnerIterator;
-
     typedef typename MatrixType::Nested MatrixTypeNested;
     typedef typename internal::remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
     typedef typename internal::remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
@@ -63,108 +60,6 @@ template<typename MatrixType, unsigned int Mode> class TriangularViewImpl<Matrix
   
 };
 
-template<typename MatrixType, unsigned int Mode>
-class TriangularViewImpl<MatrixType,Mode,Sparse>::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator
-{
-    typedef typename MatrixTypeNestedCleaned::InnerIterator Base;
-  public:
-
-    EIGEN_STRONG_INLINE InnerIterator(const TriangularViewImpl& view, Index outer)
-      : Base(view.derived().nestedExpression(), outer), m_returnOne(false)
-    {
-      if(SkipFirst)
-      {
-        while((*this) && ((HasUnitDiag||SkipDiag)  ? this->index()<=outer : this->index()<outer))
-          Base::operator++();
-        if(HasUnitDiag)
-          m_returnOne = true;
-      }
-      else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
-      {
-        if((!SkipFirst) && Base::operator bool())
-          Base::operator++();
-        m_returnOne = true;
-      }
-    }
-
-    EIGEN_STRONG_INLINE InnerIterator& operator++()
-    {
-      if(HasUnitDiag && m_returnOne)
-        m_returnOne = false;
-      else
-      {
-        Base::operator++();
-        if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer()))
-        {
-          if((!SkipFirst) && Base::operator bool())
-            Base::operator++();
-          m_returnOne = true;
-        }
-      }
-      return *this;
-    }
-
-    inline Index row() const { return (MatrixType::Flags&RowMajorBit ? Base::outer() : this->index()); }
-    inline Index col() const { return (MatrixType::Flags&RowMajorBit ? this->index() : Base::outer()); }
-    inline StorageIndex index() const
-    {
-      if(HasUnitDiag && m_returnOne)  return Base::outer();
-      else                            return Base::index();
-    }
-    inline Scalar value() const
-    {
-      if(HasUnitDiag && m_returnOne)  return Scalar(1);
-      else                            return Base::value();
-    }
-
-    EIGEN_STRONG_INLINE operator bool() const
-    {
-      if(HasUnitDiag && m_returnOne)
-        return true;
-      if(SkipFirst) return  Base::operator bool();
-      else
-      {
-        if (SkipDiag) return (Base::operator bool() && this->index() < this->outer());
-        else return (Base::operator bool() && this->index() <= this->outer());
-      }
-    }
-  protected:
-    bool m_returnOne;
-};
-
-template<typename MatrixType, unsigned int Mode>
-class TriangularViewImpl<MatrixType,Mode,Sparse>::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator
-{
-    typedef typename MatrixTypeNestedCleaned::ReverseInnerIterator Base;
-  public:
-
-    EIGEN_STRONG_INLINE ReverseInnerIterator(const TriangularViewType& view, Index outer)
-      : Base(view.derived().nestedExpression(), outer)
-    {
-      eigen_assert((!HasUnitDiag) && "ReverseInnerIterator does not support yet triangular views with a unit diagonal");
-      if(SkipLast) {
-        while((*this) && (SkipDiag ? this->index()>=outer : this->index()>outer))
-          --(*this);
-      }
-    }
-
-    EIGEN_STRONG_INLINE ReverseInnerIterator& operator--()
-    { Base::operator--(); return *this; }
-
-    inline Index row() const { return Base::row(); }
-    inline Index col() const { return Base::col(); }
-
-    EIGEN_STRONG_INLINE operator bool() const
-    {
-      if (SkipLast) return Base::operator bool() ;
-      else
-      {
-        if(SkipDiag) return (Base::operator bool() && this->index() > this->outer());
-        else return (Base::operator bool() && this->index() >= this->outer());
-      }
-    }
-};
-
 namespace internal {
 
 template<typename ArgType, unsigned int Mode>
@@ -193,7 +88,7 @@ public:
     Flags = XprType::Flags
   };
     
-  explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()) {}
+  explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()), m_arg(xpr.nestedExpression()) {}
   
   inline Index nonZerosEstimate() const {
     return m_argImpl.nonZerosEstimate();
@@ -205,20 +100,20 @@ public:
     public:
 
       EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& xprEval, Index outer)
-        : Base(xprEval.m_argImpl,outer), m_returnOne(false)
+        : Base(xprEval.m_argImpl,outer), m_returnOne(false), m_containsDiag(Base::outer()<xprEval.m_arg.innerSize())
       {
         if(SkipFirst)
         {
           while((*this) && ((HasUnitDiag||SkipDiag)  ? this->index()<=outer : this->index()<outer))
             Base::operator++();
           if(HasUnitDiag)
-            m_returnOne = true;
+            m_returnOne = m_containsDiag;
         }
         else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
         {
           if((!SkipFirst) && Base::operator bool())
             Base::operator++();
-          m_returnOne = true; // FIXME check innerSize()>outer();
+          m_returnOne = m_containsDiag;
         }
       }
 
@@ -233,7 +128,7 @@ public:
           {
             if((!SkipFirst) && Base::operator bool())
               Base::operator++();
-            m_returnOne = true; // FIXME check innerSize()>outer();
+            m_returnOne = m_containsDiag;
           }
         }
         return *this;
@@ -266,12 +161,14 @@ public:
 
     protected:
       bool m_returnOne;
+      bool m_containsDiag;
     private:
       Scalar& valueRef();
   };
   
 protected:
   evaluator<ArgType> m_argImpl;
+  const ArgType& m_arg;
 };
 
 } // end namespace internal