merge

Eigen/Core

@@ -93,7 +93,14 @@
     #endif
 
     // include files
-
+    #if (defined __GNUC__) && (defined __MINGW32__)
+      #include <intrin.h>
+      //including intrin.h works around a MINGW bug http://sourceforge.net/tracker/?func=detail&atid=102435&aid=2962480&group_id=2435
+      //in essence, intrin.h is included by windows.h and also declares intrinsics (just as emmintrin.h etc. below do).  However,
+      //intrin.h uses an extern "C" declaration, and g++ thus complains of duplicate declarations with conflicting linkage.  The linkage for intrinsics
+      //doesn't matter, but at that stage the compiler doesn't know; so, to avoid compile errors when windows.h is included after Eigen/Core,
+      //include intrin here.
+    #endif
     #include <emmintrin.h>
     #include <xmmintrin.h>
     #ifdef  EIGEN_VECTORIZE_SSE3

Eigen/src/CMakeLists.txt

@@ -1,12 +1,6 @@
-ADD_SUBDIRECTORY(Core)
-ADD_SUBDIRECTORY(LU)
-ADD_SUBDIRECTORY(QR)
-ADD_SUBDIRECTORY(SVD)
-ADD_SUBDIRECTORY(Cholesky)
-ADD_SUBDIRECTORY(Geometry)
-ADD_SUBDIRECTORY(Sparse)
-ADD_SUBDIRECTORY(Jacobi)
-ADD_SUBDIRECTORY(Householder)
-ADD_SUBDIRECTORY(Eigenvalues)
-ADD_SUBDIRECTORY(misc)
-ADD_SUBDIRECTORY(plugins)
+file(GLOB Eigen_src_subdirectories "*")
+foreach(f ${Eigen_src_subdirectories})
+  if(NOT f MATCHES ".txt")
+    add_subdirectory(${f})
+  endif()
+endforeach()

Eigen/src/Core/CwiseBinaryOp.h

@@ -216,7 +216,7 @@ EIGEN_STRONG_INLINE Derived &
 MatrixBase<Derived>::operator-=(const MatrixBase<OtherDerived> &other)
 {
   SelfCwiseBinaryOp<ei_scalar_difference_op<Scalar>, Derived, OtherDerived> tmp(derived());
-  tmp = other;
+  tmp = other.derived();
   return derived();
 }
 

Eigen/src/Core/NoAlias.h

@@ -43,6 +43,7 @@
 template<typename ExpressionType, template <typename> class StorageBase>
 class NoAlias
 {
+    typedef typename ExpressionType::Scalar Scalar;
   public:
     NoAlias(ExpressionType& expression) : m_expression(expression) {}
 
@@ -50,17 +51,31 @@ class NoAlias
       * \sa MatrixBase::lazyAssign() */
     template<typename OtherDerived>
     EIGEN_STRONG_INLINE ExpressionType& operator=(const StorageBase<OtherDerived>& other)
-    { return m_expression.lazyAssign(other.derived()); }
+    { return ei_assign_selector<ExpressionType,OtherDerived,false>::run(m_expression,other.derived()); }
 
     /** \sa MatrixBase::operator+= */
     template<typename OtherDerived>
     EIGEN_STRONG_INLINE ExpressionType& operator+=(const StorageBase<OtherDerived>& other)
-    { return m_expression.lazyAssign(m_expression + other.derived()); }
+    {
+      typedef SelfCwiseBinaryOp<ei_scalar_sum_op<Scalar>, ExpressionType, OtherDerived> SelfAdder;
+      SelfAdder tmp(m_expression);
+      typedef typename ei_nested<OtherDerived>::type OtherDerivedNested;
+      typedef typename ei_cleantype<OtherDerivedNested>::type _OtherDerivedNested;
+      ei_assign_selector<SelfAdder,_OtherDerivedNested,false>::run(tmp,OtherDerivedNested(other.derived()));
+      return m_expression;
+    }
 
     /** \sa MatrixBase::operator-= */
     template<typename OtherDerived>
     EIGEN_STRONG_INLINE ExpressionType& operator-=(const StorageBase<OtherDerived>& other)
-    { return m_expression.lazyAssign(m_expression - other.derived()); }
+    {
+      typedef SelfCwiseBinaryOp<ei_scalar_difference_op<Scalar>, ExpressionType, OtherDerived> SelfAdder;
+      SelfAdder tmp(m_expression);
+      typedef typename ei_nested<OtherDerived>::type OtherDerivedNested;
+      typedef typename ei_cleantype<OtherDerivedNested>::type _OtherDerivedNested;
+      ei_assign_selector<SelfAdder,_OtherDerivedNested,false>::run(tmp,OtherDerivedNested(other.derived()));
+      return m_expression;
+    }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename ProductDerived, typename Lhs, typename Rhs>

Eigen/src/Core/SelfCwiseBinaryOp.h

@@ -62,8 +62,6 @@ template<typename BinaryOp, typename Lhs, typename Rhs> class SelfCwiseBinaryOp
 
     typedef typename ei_packet_traits<Scalar>::type Packet;
 
-    using Base::operator=;
-
     inline SelfCwiseBinaryOp(Lhs& xpr, const BinaryOp& func = BinaryOp()) : m_matrix(xpr), m_functor(func) {}
 
     inline Index rows() const { return m_matrix.rows(); }
@@ -142,6 +140,15 @@ template<typename BinaryOp, typename Lhs, typename Rhs> class SelfCwiseBinaryOp
     #endif
       return *this;
     }
+    
+    // overloaded to honor evaluation of special matrices
+    // maybe another solution would be to not use SelfCwiseBinaryOp
+    // at first...
+    SelfCwiseBinaryOp& operator=(const Rhs& _rhs)
+    {
+      typename ei_nested<Rhs>::type rhs(_rhs);
+      return Base::operator=(rhs);
+    }
 
   protected:
     Lhs& m_matrix;

Eigen/src/Eigen2Support/CMakeLists.txt

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

Eigen/src/Geometry/OrthoMethods.h

@@ -54,7 +54,7 @@ MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
 
 template< int Arch,typename VectorLhs,typename VectorRhs,
           typename Scalar = typename VectorLhs::Scalar,
-          int Vectorizable = (VectorLhs::Flags&VectorRhs::Flags)&PacketAccessBit>
+          bool Vectorizable = (VectorLhs::Flags&VectorRhs::Flags)&PacketAccessBit>
 struct ei_cross3_impl {
   inline static typename ei_plain_matrix_type<VectorLhs>::type
   run(const VectorLhs& lhs, const VectorRhs& rhs)

Eigen/src/Geometry/Quaternion.h

@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
 // Copyright (C) 2009 Mathieu Gautier <mathieu.gautier@cea.fr>
 //
 // Eigen is free software; you can redistribute it and/or
@@ -180,6 +180,10 @@ public:
     return typename ei_cast_return_type<Derived,Quaternion<NewScalarType> >::type(
       coeffs().template cast<NewScalarType>());
   }
+  
+#ifdef EIGEN_QUATERNIONBASE_PLUGIN
+# include EIGEN_QUATERNIONBASE_PLUGIN
+#endif
 };
 
 /***************************************************************************
@@ -395,7 +399,8 @@ QuaternionBase<Derived>::_transformVector(Vector3 v) const
     // It appears to be much faster than the common algorithm found
     // in the litterature (30 versus 39 flops). It also requires two
     // Vector3 as temporaries.
-    Vector3 uv = Scalar(2) * this->vec().cross(v);
+    Vector3 uv = this->vec().cross(v);
+    uv += uv;
     return v + this->w() * uv + this->vec().cross(uv);
 }
 

Eigen/src/Geometry/arch/Geometry_SSE.h

@@ -54,8 +54,8 @@ struct ei_cross3_impl<Architecture::SSE,VectorLhs,VectorRhs,float,true>
   inline static typename ei_plain_matrix_type<VectorLhs>::type
   run(const VectorLhs& lhs, const VectorRhs& rhs)
   {
-    __m128 a = lhs.coeffs().packet<VectorLhs::Flags&AlignedBit ? Aligned : Unaligned>(0);
-    __m128 b = rhs.coeffs().packet<VectorRhs::Flags&AlignedBit ? Aligned : Unaligned>(0);
+    __m128 a = lhs.template packet<VectorLhs::Flags&AlignedBit ? Aligned : Unaligned>(0);
+    __m128 b = rhs.template packet<VectorRhs::Flags&AlignedBit ? Aligned : Unaligned>(0);
     __m128 mul1=_mm_mul_ps(ei_vec4f_swizzle1(a,1,2,0,3),ei_vec4f_swizzle1(b,2,0,1,3));
     __m128 mul2=_mm_mul_ps(ei_vec4f_swizzle1(a,2,0,1,3),ei_vec4f_swizzle1(b,1,2,0,3));
     typename ei_plain_matrix_type<VectorLhs>::type res;

Eigen/src/StlSupport/CMakeLists.txt

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

bench/geometry.cpp

@@ -0,0 +1,126 @@
+
+#include <iostream>
+#include <Eigen/Geometry>
+#include <bench/BenchTimer.h>
+
+using namespace std;
+using namespace Eigen;
+
+#ifndef SCALAR
+#define SCALAR float
+#endif
+
+#ifndef SIZE
+#define SIZE 8
+#endif
+
+typedef SCALAR Scalar;
+typedef NumTraits<Scalar>::Real RealScalar;
+typedef Matrix<RealScalar,Dynamic,Dynamic> A;
+typedef Matrix</*Real*/Scalar,Dynamic,Dynamic> B;
+typedef Matrix<Scalar,Dynamic,Dynamic> C;
+typedef Matrix<RealScalar,Dynamic,Dynamic> M;
+
+template<typename Transformation, typename Data>
+EIGEN_DONT_INLINE void transform(const Transformation& t, Data& data)
+{
+  EIGEN_ASM_COMMENT("begin");
+  data = t * data;
+  EIGEN_ASM_COMMENT("end");
+}
+
+template<typename Scalar, typename Data>
+EIGEN_DONT_INLINE void transform(const Quaternion<Scalar>& t, Data& data)
+{
+  EIGEN_ASM_COMMENT("begin quat");
+  for(int i=0;i<data.cols();++i)
+    data.col(i) = t * data.col(i);
+  EIGEN_ASM_COMMENT("end quat");
+}
+
+template<typename T> struct ToRotationMatrixWrapper
+{
+  enum {Dim = T::Dim};
+  typedef typename T::Scalar Scalar;
+  ToRotationMatrixWrapper(const T& o) : object(o) {}
+  T object;
+};
+
+template<typename QType, typename Data>
+EIGEN_DONT_INLINE void transform(const ToRotationMatrixWrapper<QType>& t, Data& data)
+{
+  EIGEN_ASM_COMMENT("begin quat via mat");
+  data = t.object.toRotationMatrix() * data;
+  EIGEN_ASM_COMMENT("end quat via mat");
+}
+
+template<typename Scalar, int Dim, typename Data>
+EIGEN_DONT_INLINE void transform(const Transform<Scalar,Dim,Projective>& t, Data& data)
+{
+  data = (t * data.colwise().homogeneous()).template block<Dim,Data::ColsAtCompileTime>(0,0);
+}
+
+template<typename T> struct get_dim { enum { Dim = T::Dim }; };
+template<typename S, int R, int C, int O, int MR, int MC>
+struct get_dim<Matrix<S,R,C,O,MR,MC> > { enum { Dim = R }; };
+
+template<typename Transformation, int N>
+struct bench_impl
+{
+  static EIGEN_DONT_INLINE void run(const Transformation& t)
+  {
+    Matrix<typename Transformation::Scalar,get_dim<Transformation>::Dim,N> data;
+    data.setRandom();
+    bench_impl<Transformation,N-1>::run(t);
+    BenchTimer timer;
+    BENCH(timer,10,100000,transform(t,data));
+    cout.width(9);
+    cout << timer.best() << " ";
+  }
+};
+
+
+template<typename Transformation>
+struct bench_impl<Transformation,0>
+{
+  static EIGEN_DONT_INLINE void run(const Transformation&) {}
+};
+
+template<typename Transformation>
+EIGEN_DONT_INLINE void bench(const std::string& msg, const Transformation& t)
+{
+  cout << msg << " ";
+  bench_impl<Transformation,SIZE>::run(t);
+  std::cout << "\n";
+}
+
+int main(int argc, char ** argv)
+{
+  Matrix<Scalar,3,4> mat34; mat34.setRandom();
+  Transform<Scalar,3,Isometry> iso3(mat34);
+  Transform<Scalar,3,Affine> aff3(mat34);
+  Transform<Scalar,3,AffineCompact> caff3(mat34);
+  Transform<Scalar,3,Projective> proj3(mat34);
+  Quaternion<Scalar> quat;quat.setIdentity();
+  ToRotationMatrixWrapper<Quaternion<Scalar> > quatmat(quat);
+  Matrix<Scalar,3,3> mat33; mat33.setRandom();
+  
+  cout.precision(4);
+  std::cout
+     << "N          ";
+  for(int i=0;i<SIZE;++i)
+  {
+    cout.width(9);
+    cout << i+1 << " ";
+  }
+  cout << "\n";
+  
+  bench("matrix 3x3", mat33);
+  bench("quaternion", quat);
+  bench("quat-mat  ", quatmat);
+  bench("isometry3 ", iso3);
+  bench("affine3   ", aff3);
+  bench("c affine3 ", caff3);
+  bench("proj3     ", proj3);
+}
+

test/basicstuff.cpp

@@ -31,6 +31,7 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   typedef typename MatrixType::Index Index;
   typedef typename MatrixType::Scalar Scalar;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
 
   Index rows = m.rows();
   Index cols = m.cols();
@@ -47,6 +48,7 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   VectorType v1 = VectorType::Random(rows),
              v2 = VectorType::Random(rows),
              vzero = VectorType::Zero(rows);
+  SquareMatrixType sm1 = SquareMatrixType::Random(rows,rows), sm2(rows,rows);
 
   Scalar x = ei_random<Scalar>();
 
@@ -121,6 +123,27 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   m1 = m2;
   VERIFY(m1==m2);
   VERIFY(!(m1!=m2));
+  
+  // check automatic transposition
+  sm2.setZero();
+  for(typename MatrixType::Index i=0;i<rows;++i)
+    sm2.col(i) = sm1.row(i);
+  VERIFY_IS_APPROX(sm2,sm1.transpose());
+  
+  sm2.setZero();
+  for(typename MatrixType::Index i=0;i<rows;++i)
+    sm2.col(i).noalias() = sm1.row(i);
+  VERIFY_IS_APPROX(sm2,sm1.transpose());
+  
+  sm2.setZero();
+  for(typename MatrixType::Index i=0;i<rows;++i)
+    sm2.col(i).noalias() += sm1.row(i);
+  VERIFY_IS_APPROX(sm2,sm1.transpose());
+  
+  sm2.setZero();
+  for(typename MatrixType::Index i=0;i<rows;++i)
+    sm2.col(i).noalias() -= sm1.row(i);
+  VERIFY_IS_APPROX(sm2,-sm1.transpose());
 }
 
 template<typename MatrixType> void basicStuffComplex(const MatrixType& m)
@@ -177,14 +200,14 @@ void test_basicstuff()
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( basicStuff(Matrix<float, 1, 1>()) );
     CALL_SUBTEST_2( basicStuff(Matrix4d()) );
-    CALL_SUBTEST_3( basicStuff(MatrixXcf(3, 3)) );
-    CALL_SUBTEST_4( basicStuff(MatrixXi(8, 12)) );
-    CALL_SUBTEST_5( basicStuff(MatrixXcd(20, 20)) );
+    CALL_SUBTEST_3( basicStuff(MatrixXcf(ei_random<int>(1,100), ei_random<int>(1,100))) );
+    CALL_SUBTEST_4( basicStuff(MatrixXi(ei_random<int>(1,100), ei_random<int>(1,100))) );
+    CALL_SUBTEST_5( basicStuff(MatrixXcd(ei_random<int>(1,100), ei_random<int>(1,100))) );
     CALL_SUBTEST_6( basicStuff(Matrix<float, 100, 100>()) );
-    CALL_SUBTEST_7( basicStuff(Matrix<long double,Dynamic,Dynamic>(10,10)) );
+    CALL_SUBTEST_7( basicStuff(Matrix<long double,Dynamic,Dynamic>(ei_random<int>(1,100),ei_random<int>(1,100))) );
 
-    CALL_SUBTEST_3( basicStuffComplex(MatrixXcf(21, 17)) );
-    CALL_SUBTEST_5( basicStuffComplex(MatrixXcd(2, 3)) );
+    CALL_SUBTEST_3( basicStuffComplex(MatrixXcf(ei_random<int>(1,100), ei_random<int>(1,100))) );
+    CALL_SUBTEST_5( basicStuffComplex(MatrixXcd(ei_random<int>(1,100), ei_random<int>(1,100))) );
   }
 
   CALL_SUBTEST_2(casting());

test/cholesky.cpp

@@ -168,6 +168,21 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     }
   }
 
+  // test some special use cases of SelfCwiseBinaryOp:
+  MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols);
+  m2 = m1;
+  m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB);
+  VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
+  m2 = m1;
+  m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
+  VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
+  m2 = m1;
+  m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB);
+  VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
+  m2 = m1;
+  m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
+  VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
+  
 }
 
 template<typename MatrixType> void cholesky_cplx(const MatrixType& m)

test/linearstructure.cpp

@@ -27,7 +27,7 @@
 template<typename MatrixType> void linearStructure(const MatrixType& m)
 {
   /* this test covers the following files:
-     Sum.h Difference.h Opposite.h ScalarMultiple.h
+     CwiseUnaryOp.h, CwiseBinaryOp.h, SelfCwiseBinaryOp.h 
   */
   typedef typename MatrixType::Index Index;
   typedef typename MatrixType::Scalar Scalar;

unsupported/Eigen/CholmodSupport

@@ -5,11 +5,9 @@
 
 #include "src/Core/util/DisableMSVCWarnings.h"
 
-#ifdef EIGEN_CHOLMOD_SUPPORT
-  extern "C" {
-    #include <cholmod.h>
-  }
-#endif
+extern "C" {
+  #include <cholmod.h>
+}
 
 namespace Eigen {
 

unsupported/Eigen/TaucsSupport

@@ -39,9 +39,7 @@ namespace Eigen {
   * \endcode
   */
 
-#ifdef EIGEN_TAUCS_SUPPORT
 # include "src/SparseExtra/TaucsSupport.h"
-#endif
 
 } // namespace Eigen