bump to 3.0.1

(transplanted from 9464745385
)

CMakeLists.txt

@@ -101,6 +101,8 @@ if(EIGEN_DEFAULT_TO_ROW_MAJOR)
   add_definitions("-DEIGEN_DEFAULT_TO_ROW_MAJOR")
 endif()
 
+add_definitions("-DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS")
+
 if(CMAKE_COMPILER_IS_GNUCXX)
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wnon-virtual-dtor -Wno-long-long -ansi -Wundef -Wcast-align -Wchar-subscripts -Wall -W -Wpointer-arith -Wwrite-strings -Wformat-security -fexceptions -fno-check-new -fno-common -fstrict-aliasing")
   set(CMAKE_CXX_FLAGS_DEBUG "-g3")
@@ -205,6 +207,7 @@ endif(MSVC)
 
 option(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION "Disable explicit vectorization in tests/examples" OFF)
 option(EIGEN_TEST_X87 "Force using X87 instructions. Implies no vectorization." OFF)
+option(EIGEN_TEST_32BIT "Force generating 32bit code." OFF)
 
 if(EIGEN_TEST_X87)
   set(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION ON)
@@ -216,6 +219,15 @@ if(EIGEN_TEST_X87)
   endif()
 endif()
 
+if(EIGEN_TEST_32BIT)
+  if(CMAKE_COMPILER_IS_GNUCXX)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -m32")
+    message(STATUS "Forcing generation of 32-bit code in tests/examples")
+  else()
+    message(STATUS "EIGEN_TEST_32BIT ignored on your compiler")
+  endif()
+endif()
+
 if(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
   add_definitions(-DEIGEN_DONT_VECTORIZE=1)
   message(STATUS "Disabling vectorization in tests/examples")
@@ -277,21 +289,52 @@ add_subdirectory(Eigen)
 
 add_subdirectory(doc EXCLUDE_FROM_ALL)
 
+add_custom_target(buildtests)
+add_custom_target(check COMMAND "ctest")
+add_dependencies(check buildtests)
+
+# CMake/Ctest does not allow us to change the build command,
+# so we have to workaround by directly editing the generated DartConfiguration.tcl file
+# save CMAKE_MAKE_PROGRAM
+set(CMAKE_MAKE_PROGRAM_SAVE ${CMAKE_MAKE_PROGRAM})
+# and set a fake one
+set(CMAKE_MAKE_PROGRAM "@EIGEN_MAKECOMMAND_PLACEHOLDER@")
+
 include(CTest)
 enable_testing() # must be called from the root CMakeLists, see man page
 include(EigenTesting)
 ei_init_testing()
 
+# overwrite default DartConfiguration.tcl
+# The worarounds are different for each version of the MSVC IDE
+if(MSVC_IDE)
+  if(MSVC_VERSION EQUAL 1600) # MSVC 2010
+    set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} buildtests.vcxproj /p:Configuration=\${CTEST_CONFIGURATION_TYPE} \n # ")
+  else() # MSVC 2008 (TODO check MSVC 2005)
+    set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} /project buildtests")
+  endif()
+else()
+  # for make and nmake
+  set(EIGEN_MAKECOMMAND_PLACEHOLDER "${CMAKE_MAKE_PROGRAM_SAVE} buildtests")
+endif()
+
+configure_file(${CMAKE_BINARY_DIR}/DartConfiguration.tcl ${CMAKE_BINARY_DIR}/DartConfiguration.tcl)
+# restore default CMAKE_MAKE_PROGRAM
+set(CMAKE_MAKE_PROGRAM ${CMAKE_MAKE_PROGRAM_SAVE})
+# un-set temporary variables so that it is like they never existed. 
+# CMake 2.6.3 introduces the more logical unset() syntax for this.
+set(CMAKE_MAKE_PROGRAM_SAVE) 
+set(EIGEN_MAKECOMMAND_PLACEHOLDER)
+
+configure_file(${CMAKE_SOURCE_DIR}/CTestCustom.cmake.in ${CMAKE_BINARY_DIR}/CTestCustom.cmake)
+
+
 if(EIGEN_LEAVE_TEST_IN_ALL_TARGET)
   add_subdirectory(test) # can't do EXCLUDE_FROM_ALL here, breaks CTest
 else()
   add_subdirectory(test EXCLUDE_FROM_ALL)
 endif()
 
-add_subdirectory(unsupported)
-
-add_subdirectory(demos EXCLUDE_FROM_ALL)
-
 if(NOT MSVC)
   if(EIGEN_LEAVE_TEST_IN_ALL_TARGET)
     add_subdirectory(blas)
@@ -302,6 +345,10 @@ if(NOT MSVC)
   endif()
 endif(NOT MSVC)
 
+add_subdirectory(unsupported)
+
+add_subdirectory(demos EXCLUDE_FROM_ALL)
+
 # must be after test and unsupported, for configuring buildtests.in
 add_subdirectory(scripts EXCLUDE_FROM_ALL)
 

CTestConfig.cmake

@@ -11,7 +11,3 @@ set(CTEST_DROP_METHOD "http")
 set(CTEST_DROP_SITE "eigen.tuxfamily.org")
 set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen")
 set(CTEST_DROP_SITE_CDASH TRUE)
-
-## A tribute to Dynamic!
-set(CTEST_CUSTOM_MAXIMUM_NUMBER_OF_WARNINGS "33331")
-set(CTEST_CUSTOM_MAXIMUM_NUMBER_OF_ERRORS "33331")

CTestCustom.cmake.in

@@ -0,0 +1,4 @@
+
+## A tribute to Dynamic!
+set(CTEST_CUSTOM_MAXIMUM_NUMBER_OF_WARNINGS "33331")
+set(CTEST_CUSTOM_MAXIMUM_NUMBER_OF_ERRORS "33331")

Eigen/Cholesky

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 namespace Eigen {
 
@@ -27,7 +27,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_CHOLESKY_MODULE_H
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

Eigen/Core

@@ -26,8 +26,8 @@
 #ifndef EIGEN_CORE_H
 #define EIGEN_CORE_H
 
-// first thing Eigen does: prevent MSVC from committing suicide
-#include "src/Core/util/DisableMSVCWarnings.h"
+// first thing Eigen does: stop the compiler from committing suicide
+#include "src/Core/util/DisableStupidWarnings.h"
 
 // then include this file where all our macros are defined. It's really important to do it first because
 // it's where we do all the alignment settings (platform detection and honoring the user's will if he
@@ -51,16 +51,16 @@
       #define EIGEN_SSE2_ON_MSVC_2008_OR_LATER
     #endif
   #endif
-#endif
-
-// Remember that usage of defined() in a #define is undefined by the standard
-#if (defined __SSE2__) && ( (!defined __GNUC__) || EIGEN_GNUC_AT_LEAST(4,2) )
-  #define EIGEN_SSE2_BUT_NOT_OLD_GCC
+#else
+  // Remember that usage of defined() in a #define is undefined by the standard
+  #if (defined __SSE2__) && ( (!defined __GNUC__) || EIGEN_GNUC_AT_LEAST(4,2) )
+    #define EIGEN_SSE2_ON_NON_MSVC_BUT_NOT_OLD_GCC
+  #endif
 #endif
 
 #ifndef EIGEN_DONT_VECTORIZE
 
-  #if defined (EIGEN_SSE2_BUT_NOT_OLD_GCC) || defined(EIGEN_SSE2_ON_MSVC_2008_OR_LATER)
+  #if defined (EIGEN_SSE2_ON_NON_MSVC_BUT_NOT_OLD_GCC) || defined(EIGEN_SSE2_ON_MSVC_2008_OR_LATER)
 
     // Defines symbols for compile-time detection of which instructions are
     // used.
@@ -143,6 +143,7 @@
 #ifdef EIGEN_HAS_ERRNO
 #include <cerrno>
 #endif
+#include <cstddef>
 #include <cstdlib>
 #include <cmath>
 #include <complex>
@@ -158,7 +159,7 @@
 
 // for outputting debug info
 #ifdef EIGEN_DEBUG_ASSIGN
-#include<iostream>
+#include <iostream>
 #endif
 
 // required for __cpuid, needs to be included after cmath
@@ -184,6 +185,7 @@
 // defined in bits/termios.h
 #undef B0
 
+/** \brief Namespace containing all symbols from the %Eigen library. */
 namespace Eigen {
 
 inline static const char *SimdInstructionSetsInUse(void) {
@@ -239,6 +241,8 @@ inline static const char *SimdInstructionSetsInUse(void) {
 // we use size_t frequently and we'll never remember to prepend it with std:: everytime just to
 // ensure QNX/QCC support
 using std::size_t;
+// gcc 4.6.0 wants std:: for ptrdiff_t 
+using std::ptrdiff_t;
 
 /** \defgroup Core_Module Core module
   * This is the main module of Eigen providing dense matrix and vector support
@@ -354,7 +358,7 @@ using std::size_t;
 
 #include "src/Core/GlobalFunctions.h"
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #ifdef EIGEN2_SUPPORT
 #include "Eigen2Support"

Eigen/Eigen2Support

@@ -29,7 +29,7 @@
 #error Eigen2 support must be enabled by defining EIGEN2_SUPPORT before including any Eigen header
 #endif
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 namespace Eigen {
 
@@ -58,7 +58,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 // Eigen2 used to include iostream
 #include<iostream>

Eigen/Eigenvalues

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 #include "Cholesky"
 #include "Jacobi"
@@ -38,7 +38,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_EIGENVALUES_MODULE_H
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

Eigen/Geometry

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 #include "SVD"
 #include "LU"
@@ -60,7 +60,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_GEOMETRY_MODULE_H
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

Eigen/Householder

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 namespace Eigen {
 
@@ -21,7 +21,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_HOUSEHOLDER_MODULE_H
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

Eigen/Jacobi

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 namespace Eigen {
 
@@ -23,7 +23,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_JACOBI_MODULE_H
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

Eigen/LU

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 namespace Eigen {
 
@@ -36,7 +36,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_LU_MODULE_H
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

Eigen/LeastSquares

@@ -5,9 +5,12 @@
 #error LeastSquares is only available in Eigen2 support mode (define EIGEN2_SUPPORT)
 #endif
 
+// exclude from normal eigen3-only documentation
+#ifdef EIGEN2_SUPPORT
+
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 #include "Eigenvalues"
 #include "Geometry"
@@ -26,6 +29,8 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif // EIGEN2_SUPPORT
 
 #endif // EIGEN_REGRESSION_MODULE_H

Eigen/QR

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 #include "Cholesky"
 #include "Jacobi"
@@ -35,7 +35,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #ifdef EIGEN2_SUPPORT
 #include "Eigenvalues"

Eigen/QtAlignedMalloc

@@ -6,7 +6,7 @@
 
 #if (!EIGEN_MALLOC_ALREADY_ALIGNED)
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 void *qMalloc(size_t size)
 {
@@ -26,7 +26,7 @@ void *qRealloc(void *ptr, size_t size)
   return newPtr;
 }
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif
 

Eigen/SVD

@@ -5,7 +5,7 @@
 #include "Householder"
 #include "Jacobi"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 namespace Eigen {
 
@@ -32,7 +32,7 @@ namespace Eigen {
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_SVD_MODULE_H
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

Eigen/Sparse

@@ -3,7 +3,7 @@
 
 #include "Core"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "src/Core/util/DisableStupidWarnings.h"
 
 #include <vector>
 #include <map>
@@ -63,7 +63,7 @@ struct Sparse {};
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "src/Core/util/ReenableStupidWarnings.h"
 
 #endif // EIGEN_SPARSE_MODULE_H
 

Eigen/src/Core/Array.h

@@ -37,6 +37,9 @@
   * API for the %Matrix class provides easy access to linear-algebra
   * operations.
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_ARRAY_PLUGIN.
+  *
   * \sa \ref TutorialArrayClass, \ref TopicClassHierarchy
   */
 namespace internal {

Eigen/src/Core/ArrayBase.h

@@ -42,7 +42,10 @@ template<typename ExpressionType> class MatrixWrapper;
   *
   * This class is the base that is inherited by all array expression types.
   *
-  * \param Derived is the derived type, e.g., an array or an expression type.
+  * \tparam Derived is the derived type, e.g., an array or an expression type.
+  *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_ARRAYBASE_PLUGIN.
   *
   * \sa class MatrixBase, \ref TopicClassHierarchy
   */

Eigen/src/Core/ArrayWrapper.h

@@ -53,6 +53,12 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
     EIGEN_DENSE_PUBLIC_INTERFACE(ArrayWrapper)
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(ArrayWrapper)
 
+    typedef typename internal::conditional<
+                       internal::is_lvalue<ExpressionType>::value,
+                       Scalar,
+                       const Scalar
+                     >::type ScalarWithConstIfNotLvalue;
+
     typedef typename internal::nested<ExpressionType>::type NestedExpressionType;
 
     inline ArrayWrapper(const ExpressionType& matrix) : m_expression(matrix) {}
@@ -62,6 +68,9 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
     inline Index outerStride() const { return m_expression.outerStride(); }
     inline Index innerStride() const { return m_expression.innerStride(); }
 
+    inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
+    inline const Scalar* data() const { return m_expression.data(); }
+
     inline const CoeffReturnType coeff(Index row, Index col) const
     {
       return m_expression.coeff(row, col);
@@ -151,6 +160,12 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
     EIGEN_DENSE_PUBLIC_INTERFACE(MatrixWrapper)
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(MatrixWrapper)
 
+    typedef typename internal::conditional<
+                       internal::is_lvalue<ExpressionType>::value,
+                       Scalar,
+                       const Scalar
+                     >::type ScalarWithConstIfNotLvalue;
+
     typedef typename internal::nested<ExpressionType>::type NestedExpressionType;
 
     inline MatrixWrapper(const ExpressionType& matrix) : m_expression(matrix) {}
@@ -160,6 +175,9 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
     inline Index outerStride() const { return m_expression.outerStride(); }
     inline Index innerStride() const { return m_expression.innerStride(); }
 
+    inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
+    inline const Scalar* data() const { return m_expression.data(); }
+
     inline const CoeffReturnType coeff(Index row, Index col) const
     {
       return m_expression.coeff(row, col);

Eigen/src/Core/Assign.h

@@ -41,7 +41,7 @@ public:
     DstIsAligned = Derived::Flags & AlignedBit,
     DstHasDirectAccess = Derived::Flags & DirectAccessBit,
     SrcIsAligned = OtherDerived::Flags & AlignedBit,
-    JointAlignment = DstIsAligned && SrcIsAligned ? Aligned : Unaligned
+    JointAlignment = bool(DstIsAligned) && bool(SrcIsAligned) ? Aligned : Unaligned
   };
 
 private:
@@ -106,9 +106,9 @@ public:
                                              : int(NoUnrolling)
                   )
               : int(Traversal) == int(LinearVectorizedTraversal)
-                ? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) )
+                ? ( bool(MayUnrollCompletely) && bool(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) )
               : int(Traversal) == int(LinearTraversal)
-                ? ( int(MayUnrollCompletely) ? int(CompleteUnrolling) : int(NoUnrolling) )
+                ? ( bool(MayUnrollCompletely) ? int(CompleteUnrolling) : int(NoUnrolling) )
               : int(NoUnrolling)
   };
 
@@ -474,7 +474,7 @@ struct assign_impl<Derived1, Derived2, SliceVectorizedTraversal, NoUnrolling>
 
       // do the vectorizable part of the assignment
       for(Index inner = alignedStart; inner<alignedEnd; inner+=packetSize)
-        dst.template copyPacketByOuterInner<Derived2, Aligned, Unaligned>(outer, inner, src);
+        dst.template copyPacketByOuterInner<Derived2, dstAlignment, Unaligned>(outer, inner, src);
 
       // do the non-vectorizable part of the assignment
       for(Index inner = alignedEnd; inner<innerSize ; ++inner)

Eigen/src/Core/BandMatrix.h

@@ -180,7 +180,7 @@ class BandMatrixBase : public EigenBase<Derived>
   * \param Cols Number of columns, or \b Dynamic
   * \param Supers Number of super diagonal
   * \param Subs Number of sub diagonal
-  * \param _Options A combination of either \b RowMajor or \b ColMajor, and of \b SelfAdjoint
+  * \param _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
   *                 The former controls \ref TopicStorageOrders "storage order", and defaults to
   *                 column-major. The latter controls whether the matrix represents a selfadjoint 
   *                 matrix in which case either Supers of Subs have to be null.

Eigen/src/Core/DenseBase.h

@@ -34,7 +34,10 @@
   * This class is the base that is inherited by all dense objects (matrix, vector, arrays,
   * and related expression types). The common Eigen API for dense objects is contained in this class.
   *
-  * \param Derived is the derived type, e.g., a matrix type or an expression.
+  * \tparam Derived is the derived type, e.g., a matrix type or an expression.
+  *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_DENSEBASE_PLUGIN.
   *
   * \sa \ref TopicClassHierarchy
   */
@@ -53,7 +56,13 @@ template<typename Derived> class DenseBase
     class InnerIterator;
 
     typedef typename internal::traits<Derived>::StorageKind StorageKind;
-    typedef typename internal::traits<Derived>::Index Index; /**< The type of indices */
+
+    /** \brief The type of indices 
+      * \details To change this, \c \#define the preprocessor symbol \c EIGEN_DEFAULT_DENSE_INDEX_TYPE.
+      * \sa \ref TopicPreprocessorDirectives.
+      */
+    typedef typename internal::traits<Derived>::Index Index; 
+
     typedef typename internal::traits<Derived>::Scalar Scalar;
     typedef typename internal::packet_traits<Scalar>::type PacketScalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;

Eigen/src/Core/DenseCoeffsBase.h

@@ -35,7 +35,7 @@ template<typename T> struct add_const_on_value_type_if_arithmetic
 /** \brief Base class providing read-only coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam ReadOnlyAccessors Constant indicating read-only access
+  * \tparam #ReadOnlyAccessors Constant indicating read-only access
   *
   * This class defines the \c operator() \c const function and friends, which can be used to read specific
   * entries of a matrix or array.
@@ -212,7 +212,7 @@ class DenseCoeffsBase<Derived,ReadOnlyAccessors> : public EigenBase<Derived>
       * to ensure that a packet really starts there. This method is only available on expressions having the
       * PacketAccessBit.
       *
-      * The \a LoadMode parameter may have the value \a Aligned or \a Unaligned. Its effect is to select
+      * The \a LoadMode parameter may have the value \a #Aligned or \a #Unaligned. Its effect is to select
       * the appropriate vectorization instruction. Aligned access is faster, but is only possible for packets
       * starting at an address which is a multiple of the packet size.
       */
@@ -239,7 +239,7 @@ class DenseCoeffsBase<Derived,ReadOnlyAccessors> : public EigenBase<Derived>
       * to ensure that a packet really starts there. This method is only available on expressions having the
       * PacketAccessBit and the LinearAccessBit.
       *
-      * The \a LoadMode parameter may have the value \a Aligned or \a Unaligned. Its effect is to select
+      * The \a LoadMode parameter may have the value \a #Aligned or \a #Unaligned. Its effect is to select
       * the appropriate vectorization instruction. Aligned access is faster, but is only possible for packets
       * starting at an address which is a multiple of the packet size.
       */
@@ -275,7 +275,7 @@ class DenseCoeffsBase<Derived,ReadOnlyAccessors> : public EigenBase<Derived>
 /** \brief Base class providing read/write coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam WriteAccessors Constant indicating read/write access
+  * \tparam #WriteAccessors Constant indicating read/write access
   *
   * This class defines the non-const \c operator() function and friends, which can be used to write specific
   * entries of a matrix or array. This class inherits DenseCoeffsBase<Derived, ReadOnlyAccessors> which
@@ -433,7 +433,7 @@ class DenseCoeffsBase<Derived, WriteAccessors> : public DenseCoeffsBase<Derived,
       * to ensure that a packet really starts there. This method is only available on expressions having the
       * PacketAccessBit.
       *
-      * The \a LoadMode parameter may have the value \a Aligned or \a Unaligned. Its effect is to select
+      * The \a LoadMode parameter may have the value \a #Aligned or \a #Unaligned. Its effect is to select
       * the appropriate vectorization instruction. Aligned access is faster, but is only possible for packets
       * starting at an address which is a multiple of the packet size.
       */
@@ -567,7 +567,7 @@ class DenseCoeffsBase<Derived, WriteAccessors> : public DenseCoeffsBase<Derived,
 /** \brief Base class providing direct read-only coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam DirectAccessors Constant indicating direct access
+  * \tparam #DirectAccessors Constant indicating direct access
   *
   * This class defines functions to work with strides which can be used to access entries directly. This class
   * inherits DenseCoeffsBase<Derived, ReadOnlyAccessors> which defines functions to access entries read-only using
@@ -637,7 +637,7 @@ class DenseCoeffsBase<Derived, DirectAccessors> : public DenseCoeffsBase<Derived
 /** \brief Base class providing direct read/write coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam DirectAccessors Constant indicating direct access
+  * \tparam #DirectWriteAccessors Constant indicating direct access
   *
   * This class defines functions to work with strides which can be used to access entries directly. This class
   * inherits DenseCoeffsBase<Derived, WriteAccessors> which defines functions to access entries read/write using

Eigen/src/Core/DenseStorage.h

@@ -58,14 +58,14 @@ struct plain_array
   #define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask) \
     eigen_assert((reinterpret_cast<size_t>(array) & sizemask) == 0 \
               && "this assertion is explained here: " \
-              "http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html" \
+              "http://eigen.tuxfamily.org/dox-devel/TopicUnalignedArrayAssert.html" \
               " **** READ THIS WEB PAGE !!! ****");
 #endif
 
 template <typename T, int Size, int MatrixOrArrayOptions>
 struct plain_array<T, Size, MatrixOrArrayOptions, 16>
 {
-  EIGEN_ALIGN16 T array[Size];
+  EIGEN_USER_ALIGN16 T array[Size];
   plain_array() { EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(0xf) }
   plain_array(constructor_without_unaligned_array_assert) {}
 };
@@ -73,7 +73,7 @@ struct plain_array<T, Size, MatrixOrArrayOptions, 16>
 template <typename T, int MatrixOrArrayOptions, int Alignment>
 struct plain_array<T, 0, MatrixOrArrayOptions, Alignment>
 {
-  EIGEN_ALIGN16 T array[1];
+  EIGEN_USER_ALIGN16 T array[1];
   plain_array() {}
   plain_array(constructor_without_unaligned_array_assert) {}
 };

Eigen/src/Core/Functors.h

@@ -116,7 +116,7 @@ struct functor_traits<scalar_conj_product_op<LhsScalar,RhsScalar> > {
   */
 template<typename Scalar> struct scalar_min_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_min_op)
-  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return std::min(a, b); }
+  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { using std::min; return min(a, b); }
   template<typename Packet>
   EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
   { return internal::pmin(a,b); }
@@ -139,7 +139,7 @@ struct functor_traits<scalar_min_op<Scalar> > {
   */
 template<typename Scalar> struct scalar_max_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_max_op)
-  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return std::max(a, b); }
+  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { using std::max; return max(a, b); }
   template<typename Packet>
   EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
   { return internal::pmax(a,b); }
@@ -165,8 +165,10 @@ template<typename Scalar> struct scalar_hypot_op {
 //   typedef typename NumTraits<Scalar>::Real result_type;
   EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const
   {
-    Scalar p = std::max(_x, _y);
-    Scalar q = std::min(_x, _y);
+    using std::max;
+    using std::min;
+    Scalar p = max(_x, _y);
+    Scalar q = min(_x, _y);
     Scalar qp = q/p;
     return p * sqrt(Scalar(1) + qp*qp);
   }
@@ -605,7 +607,7 @@ template <typename Scalar, bool RandomAccess> struct linspaced_op
   EIGEN_STRONG_INLINE const Packet packetOp(Index row, Index col) const
   {
     eigen_assert(col==0 || row==0);
-    return impl(col + row);
+    return impl.packetOp(col + row);
   }
 
   // This proxy object handles the actual required temporaries, the different
@@ -669,7 +671,7 @@ struct functor_traits<scalar_sqrt_op<Scalar> >
 
 /** \internal
   * \brief Template functor to compute the cosine of a scalar
-  * \sa class CwiseUnaryOp, Cwise::cos()
+  * \sa class CwiseUnaryOp, ArrayBase::cos()
   */
 template<typename Scalar> struct scalar_cos_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_cos_op)
@@ -688,7 +690,7 @@ struct functor_traits<scalar_cos_op<Scalar> >
 
 /** \internal
   * \brief Template functor to compute the sine of a scalar
-  * \sa class CwiseUnaryOp, Cwise::sin()
+  * \sa class CwiseUnaryOp, ArrayBase::sin()
   */
 template<typename Scalar> struct scalar_sin_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_sin_op)
@@ -708,7 +710,7 @@ struct functor_traits<scalar_sin_op<Scalar> >
 
 /** \internal
   * \brief Template functor to compute the tan of a scalar
-  * \sa class CwiseUnaryOp, Cwise::tan()
+  * \sa class CwiseUnaryOp, ArrayBase::tan()
   */
 template<typename Scalar> struct scalar_tan_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_tan_op)
@@ -725,6 +727,44 @@ struct functor_traits<scalar_tan_op<Scalar> >
   };
 };
 
+/** \internal
+  * \brief Template functor to compute the arc cosine of a scalar
+  * \sa class CwiseUnaryOp, ArrayBase::acos()
+  */
+template<typename Scalar> struct scalar_acos_op {
+  EIGEN_EMPTY_STRUCT_CTOR(scalar_acos_op)
+  inline const Scalar operator() (const Scalar& a) const { return acos(a); }
+  typedef typename packet_traits<Scalar>::type Packet;
+  inline Packet packetOp(const Packet& a) const { return internal::pacos(a); }
+};
+template<typename Scalar>
+struct functor_traits<scalar_acos_op<Scalar> >
+{
+  enum {
+    Cost = 5 * NumTraits<Scalar>::MulCost,
+    PacketAccess = packet_traits<Scalar>::HasACos
+  };
+};
+
+/** \internal
+  * \brief Template functor to compute the arc sine of a scalar
+  * \sa class CwiseUnaryOp, ArrayBase::asin()
+  */
+template<typename Scalar> struct scalar_asin_op {
+  EIGEN_EMPTY_STRUCT_CTOR(scalar_asin_op)
+  inline const Scalar operator() (const Scalar& a) const { return asin(a); }
+  typedef typename packet_traits<Scalar>::type Packet;
+  inline Packet packetOp(const Packet& a) const { return internal::pasin(a); }
+};
+template<typename Scalar>
+struct functor_traits<scalar_asin_op<Scalar> >
+{
+  enum {
+    Cost = 5 * NumTraits<Scalar>::MulCost,
+    PacketAccess = packet_traits<Scalar>::HasASin
+  };
+};
+
 /** \internal
   * \brief Template functor to raise a scalar to a power
   * \sa class CwiseUnaryOp, Cwise::pow

Eigen/src/Core/Fuzzy.h

@@ -34,9 +34,10 @@ struct isApprox_selector
 {
   static bool run(const Derived& x, const OtherDerived& y, typename Derived::RealScalar prec)
   {
+    using std::min;
     const typename internal::nested<Derived,2>::type nested(x);
     const typename internal::nested<OtherDerived,2>::type otherNested(y);
-    return (nested - otherNested).cwiseAbs2().sum() <= prec * prec * std::min(nested.cwiseAbs2().sum(), otherNested.cwiseAbs2().sum());
+    return (nested - otherNested).cwiseAbs2().sum() <= prec * prec * min(nested.cwiseAbs2().sum(), otherNested.cwiseAbs2().sum());
   }
 };
 

Eigen/src/Core/GenericPacketMath.h

@@ -134,12 +134,12 @@ pdiv(const Packet& a,
 /** \internal \returns the min of \a a and \a b  (coeff-wise) */
 template<typename Packet> inline Packet
 pmin(const Packet& a,
-        const Packet& b) { return std::min(a, b); }
+        const Packet& b) { using std::min; return min(a, b); }
 
 /** \internal \returns the max of \a a and \a b  (coeff-wise) */
 template<typename Packet> inline Packet
 pmax(const Packet& a,
-        const Packet& b) { return std::max(a, b); }
+        const Packet& b) { using std::max; return max(a, b); }
 
 /** \internal \returns the absolute value of \a a */
 template<typename Packet> inline Packet
@@ -225,15 +225,20 @@ template<typename Packet> inline typename unpacket_traits<Packet>::type predux_m
 template<typename Packet> inline Packet preverse(const Packet& a)
 { return a; }
 
+
+/** \internal \returns \a a with real and imaginary part flipped (for complex type only) */
+template<typename Packet> inline Packet pcplxflip(const Packet& a)
+{ return Packet(imag(a),real(a)); }
+
 /**************************
 * Special math functions
 ***************************/
 
-/** \internal \returns the sin of \a a (coeff-wise) */
+/** \internal \returns the sine of \a a (coeff-wise) */
 template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 Packet psin(const Packet& a) { return sin(a); }
 
-/** \internal \returns the cos of \a a (coeff-wise) */
+/** \internal \returns the cosine of \a a (coeff-wise) */
 template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 Packet pcos(const Packet& a) { return cos(a); }
 
@@ -241,6 +246,14 @@ Packet pcos(const Packet& a) { return cos(a); }
 template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 Packet ptan(const Packet& a) { return tan(a); }
 
+/** \internal \returns the arc sine of \a a (coeff-wise) */
+template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pasin(const Packet& a) { return asin(a); }
+
+/** \internal \returns the arc cosine of \a a (coeff-wise) */
+template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pacos(const Packet& a) { return acos(a); }
+
 /** \internal \returns the exp of \a a (coeff-wise) */
 template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 Packet pexp(const Packet& a) { return exp(a); }
@@ -257,6 +270,14 @@ Packet psqrt(const Packet& a) { return sqrt(a); }
 * The following functions might not have to be overwritten for vectorized types
 ***************************************************************************/
 
+/** \internal copy a packet with constant coeficient \a a (e.g., [a,a,a,a]) to \a *to. \a to must be 16 bytes aligned */
+// NOTE: this function must really be templated on the packet type (think about different packet types for the same scalar type)
+template<typename Packet>
+inline void pstore1(typename unpacket_traits<Packet>::type* to, const typename unpacket_traits<Packet>::type& a)
+{
+  pstore(to, pset1<Packet>(a));
+}
+
 /** \internal \returns a * b + c (coeff-wise) */
 template<typename Packet> inline Packet
 pmadd(const Packet&  a,
@@ -265,7 +286,7 @@ pmadd(const Packet&  a,
 { return padd(pmul(a, b),c); }
 
 /** \internal \returns a packet version of \a *from.
-  * \If LoadMode equals Aligned, \a from must be 16 bytes aligned */
+  * If LoadMode equals #Aligned, \a from must be 16 bytes aligned */
 template<typename Packet, int LoadMode>
 inline Packet ploadt(const typename unpacket_traits<Packet>::type* from)
 {
@@ -276,7 +297,7 @@ inline Packet ploadt(const typename unpacket_traits<Packet>::type* from)
 }
 
 /** \internal copy the packet \a from to \a *to.
-  * If StoreMode equals Aligned, \a to must be 16 bytes aligned */
+  * If StoreMode equals #Aligned, \a to must be 16 bytes aligned */
 template<typename Scalar, typename Packet, int LoadMode>
 inline void pstoret(Scalar* to, const Packet& from)
 {

Eigen/src/Core/GlobalFunctions.h

@@ -56,6 +56,8 @@ namespace std
   EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(imag,scalar_imag_op)
   EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(sin,scalar_sin_op)
   EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(cos,scalar_cos_op)
+  EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(asin,scalar_asin_op)
+  EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(acos,scalar_acos_op)
   EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(tan,scalar_tan_op)
   EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(exp,scalar_exp_op)
   EIGEN_ARRAY_DECLARE_GLOBAL_STD_UNARY(log,scalar_log_op)
@@ -77,6 +79,8 @@ namespace Eigen
     EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(imag,scalar_imag_op)
     EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(sin,scalar_sin_op)
     EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(cos,scalar_cos_op)
+    EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(asin,scalar_asin_op)
+    EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(acos,scalar_acos_op)
     EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(tan,scalar_tan_op)
     EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(exp,scalar_exp_op)
     EIGEN_ARRAY_DECLARE_GLOBAL_EIGEN_UNARY(log,scalar_log_op)

Eigen/src/Core/IO.h

@@ -141,7 +141,8 @@ struct significant_decimals_default_impl
   typedef typename NumTraits<Scalar>::Real RealScalar;
   static inline int run()
   {
-    return cast<RealScalar,int>(std::ceil(-log(NumTraits<RealScalar>::epsilon())/log(RealScalar(10))));
+    using std::ceil;
+    return cast<RealScalar,int>(ceil(-log(NumTraits<RealScalar>::epsilon())/log(RealScalar(10))));
   }
 };
 

Eigen/src/Core/Map.h

@@ -31,10 +31,10 @@
   *
   * \brief A matrix or vector expression mapping an existing array of data.
   *
-  * \param PlainObjectType the equivalent matrix type of the mapped data
-  * \param MapOptions specifies whether the pointer is \c Aligned, or \c Unaligned.
-  *                The default is \c Unaligned.
-  * \param StrideType optionnally specifies strides. By default, Map assumes the memory layout
+  * \tparam PlainObjectType the equivalent matrix type of the mapped data
+  * \tparam MapOptions specifies whether the pointer is \c #Aligned, or \c #Unaligned.
+  *                The default is \c #Unaligned.
+  * \tparam StrideType optionnally specifies strides. By default, Map assumes the memory layout
   *                   of an ordinary, contiguous array. This can be overridden by specifying strides.
   *                   The type passed here must be a specialization of the Stride template, see examples below.
   *
@@ -95,7 +95,7 @@ struct traits<Map<PlainObjectType, MapOptions, StrideType> >
     HasNoInnerStride = InnerStrideAtCompileTime == 1,
     HasNoOuterStride = StrideType::OuterStrideAtCompileTime == 0,
     HasNoStride = HasNoInnerStride && HasNoOuterStride,
-    IsAligned = int(int(MapOptions)&Aligned)==Aligned,
+    IsAligned = bool(EIGEN_ALIGN) && ((int(MapOptions)&Aligned)==Aligned),
     IsDynamicSize = PlainObjectType::SizeAtCompileTime==Dynamic,
     KeepsPacketAccess = bool(HasNoInnerStride)
                         && ( bool(IsDynamicSize)
@@ -192,14 +192,14 @@ template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int
 inline Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
   ::Array(const Scalar *data)
 {
-  _set_noalias(Eigen::Map<const Array>(data));
+  this->_set_noalias(Eigen::Map<const Array>(data));
 }
 
 template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
 inline Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
   ::Matrix(const Scalar *data)
 {
-  _set_noalias(Eigen::Map<const Matrix>(data));
+  this->_set_noalias(Eigen::Map<const Matrix>(data));
 }
 
 #endif // EIGEN_MAP_H

Eigen/src/Core/MapBase.h

@@ -85,6 +85,8 @@ template<typename Derived> class MapBase<Derived, ReadOnlyAccessors>
     using Base::rowStride;
     using Base::colStride;
 
+    // bug 217 - compile error on ICC 11.1
+    using Base::operator=;
 
     typedef typename Base::CoeffReturnType CoeffReturnType;
 
@@ -236,7 +238,7 @@ template<typename Derived> class MapBase<Derived, WriteAccessors>
                 (this->m_data + index * innerStride(), x);
     }
 
-    inline MapBase(PointerType data) : Base(data) {}
+    explicit inline MapBase(PointerType data) : Base(data) {}
     inline MapBase(PointerType data, Index size) : Base(data, size) {}
     inline MapBase(PointerType data, Index rows, Index cols) : Base(data, rows, cols) {}
 

Eigen/src/Core/MathFunctions.h

@@ -87,7 +87,8 @@ struct real_impl<std::complex<RealScalar> >
 {
   static inline RealScalar run(const std::complex<RealScalar>& x)
   {
-    return std::real(x);
+    using std::real;
+    return real(x);
   }
 };
 
@@ -122,7 +123,8 @@ struct imag_impl<std::complex<RealScalar> >
 {
   static inline RealScalar run(const std::complex<RealScalar>& x)
   {
-    return std::imag(x);
+    using std::imag;
+    return imag(x);
   }
 };
 
@@ -244,7 +246,8 @@ struct conj_impl<std::complex<RealScalar> >
 {
   static inline std::complex<RealScalar> run(const std::complex<RealScalar>& x)
   {
-    return std::conj(x);
+    using std::conj;
+    return conj(x);
   }
 };
 
@@ -270,7 +273,8 @@ struct abs_impl
   typedef typename NumTraits<Scalar>::Real RealScalar;
   static inline RealScalar run(const Scalar& x)
   {
-    return std::abs(x);
+    using std::abs;
+    return abs(x);
   }
 };
 
@@ -305,7 +309,8 @@ struct abs2_impl<std::complex<RealScalar> >
 {
   static inline RealScalar run(const std::complex<RealScalar>& x)
   {
-    return std::norm(x);
+    using std::norm;
+    return norm(x);
   }
 };
 
@@ -369,10 +374,12 @@ struct hypot_impl
   typedef typename NumTraits<Scalar>::Real RealScalar;
   static inline RealScalar run(const Scalar& x, const Scalar& y)
   {
+    using std::max;
+    using std::min;
     RealScalar _x = abs(x);
     RealScalar _y = abs(y);
-    RealScalar p = std::max(_x, _y);
-    RealScalar q = std::min(_x, _y);
+    RealScalar p = max(_x, _y);
+    RealScalar q = min(_x, _y);
     RealScalar qp = q/p;
     return p * sqrt(RealScalar(1) + qp*qp);
   }
@@ -420,7 +427,8 @@ struct sqrt_default_impl
 {
   static inline Scalar run(const Scalar& x)
   {
-    return std::sqrt(x);
+    using std::sqrt;
+    return sqrt(x);
   }
 };
 
@@ -454,194 +462,36 @@ inline EIGEN_MATHFUNC_RETVAL(sqrt, Scalar) sqrt(const Scalar& x)
 }
 
 /****************************************************************************
-* Implementation of exp                                                  *
+* Implementation of standard unary real functions (exp, log, sin, cos, ...  *
 ****************************************************************************/
 
-template<typename Scalar, bool IsInteger>
-struct exp_default_impl
-{
-  static inline Scalar run(const Scalar& x)
-  {
-    return std::exp(x);
+// This macro instanciate all the necessary template mechanism which is common to all unary real functions.
+#define EIGEN_MATHFUNC_STANDARD_REAL_UNARY(NAME) \
+  template<typename Scalar, bool IsInteger> struct NAME##_default_impl {            \
+    static inline Scalar run(const Scalar& x) { using std::NAME; return NAME(x); }  \
+  };                                                                                \
+  template<typename Scalar> struct NAME##_default_impl<Scalar, true> {              \
+    static inline Scalar run(const Scalar&) {                                       \
+      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)                                       \
+      return Scalar(0);                                                             \
+    }                                                                               \
+  };                                                                                \
+  template<typename Scalar> struct NAME##_impl                                      \
+    : NAME##_default_impl<Scalar, NumTraits<Scalar>::IsInteger>                     \
+  {};                                                                               \
+  template<typename Scalar> struct NAME##_retval { typedef Scalar type; };          \
+  template<typename Scalar>                                                         \
+  inline EIGEN_MATHFUNC_RETVAL(NAME, Scalar) NAME(const Scalar& x) {                \
+    return EIGEN_MATHFUNC_IMPL(NAME, Scalar)::run(x);                               \
   }
-};
 
-template<typename Scalar>
-struct exp_default_impl<Scalar, true>
-{
-  static inline Scalar run(const Scalar&)
-  {
-    EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
-    return Scalar(0);
-  }
-};
-
-template<typename Scalar>
-struct exp_impl : exp_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
-
-template<typename Scalar>
-struct exp_retval
-{
-  typedef Scalar type;
-};
-
-template<typename Scalar>
-inline EIGEN_MATHFUNC_RETVAL(exp, Scalar) exp(const Scalar& x)
-{
-  return EIGEN_MATHFUNC_IMPL(exp, Scalar)::run(x);
-}
-
-/****************************************************************************
-* Implementation of cos                                                  *
-****************************************************************************/
-
-template<typename Scalar, bool IsInteger>
-struct cos_default_impl
-{
-  static inline Scalar run(const Scalar& x)
-  {
-    return std::cos(x);
-  }
-};
-
-template<typename Scalar>
-struct cos_default_impl<Scalar, true>
-{
-  static inline Scalar run(const Scalar&)
-  {
-    EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
-    return Scalar(0);
-  }
-};
-
-template<typename Scalar>
-struct cos_impl : cos_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
-
-template<typename Scalar>
-struct cos_retval
-{
-  typedef Scalar type;
-};
-
-template<typename Scalar>
-inline EIGEN_MATHFUNC_RETVAL(cos, Scalar) cos(const Scalar& x)
-{
-  return EIGEN_MATHFUNC_IMPL(cos, Scalar)::run(x);
-}
-
-/****************************************************************************
-* Implementation of sin                                                  *
-****************************************************************************/
-
-template<typename Scalar, bool IsInteger>
-struct sin_default_impl
-{
-  static inline Scalar run(const Scalar& x)
-  {
-    return std::sin(x);
-  }
-};
-
-template<typename Scalar>
-struct sin_default_impl<Scalar, true>
-{
-  static inline Scalar run(const Scalar&)
-  {
-    EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
-    return Scalar(0);
-  }
-};
-
-template<typename Scalar>
-struct sin_impl : sin_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
-
-template<typename Scalar>
-struct sin_retval
-{
-  typedef Scalar type;
-};
-
-template<typename Scalar>
-inline EIGEN_MATHFUNC_RETVAL(sin, Scalar) sin(const Scalar& x)
-{
-  return EIGEN_MATHFUNC_IMPL(sin, Scalar)::run(x);
-}
-
-/****************************************************************************
-* Implementation of tan                                                  *
-****************************************************************************/
-
-template<typename Scalar, bool IsInteger>
-struct tan_default_impl
-{
-  static inline Scalar run(const Scalar& x)
-  {
-    return std::tan(x);
-  }
-};
-
-template<typename Scalar>
-struct tan_default_impl<Scalar, true>
-{
-  static inline Scalar run(const Scalar&)
-  {
-    EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
-    return Scalar(0);
-  }
-};
-
-template<typename Scalar>
-struct tan_impl : tan_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
-
-template<typename Scalar>
-struct tan_retval
-{
-  typedef Scalar type;
-};
-
-template<typename Scalar>
-inline EIGEN_MATHFUNC_RETVAL(tan, Scalar) tan(const Scalar& x)
-{
-  return EIGEN_MATHFUNC_IMPL(tan, Scalar)::run(x);
-}
-
-/****************************************************************************
-* Implementation of log                                                  *
-****************************************************************************/
-
-template<typename Scalar, bool IsInteger>
-struct log_default_impl
-{
-  static inline Scalar run(const Scalar& x)
-  {
-    return std::log(x);
-  }
-};
-
-template<typename Scalar>
-struct log_default_impl<Scalar, true>
-{
-  static inline Scalar run(const Scalar&)
-  {
-    EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
-    return Scalar(0);
-  }
-};
-
-template<typename Scalar>
-struct log_impl : log_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
-
-template<typename Scalar>
-struct log_retval
-{
-  typedef Scalar type;
-};
-
-template<typename Scalar>
-inline EIGEN_MATHFUNC_RETVAL(log, Scalar) log(const Scalar& x)
-{
-  return EIGEN_MATHFUNC_IMPL(log, Scalar)::run(x);
-}
+EIGEN_MATHFUNC_STANDARD_REAL_UNARY(exp)
+EIGEN_MATHFUNC_STANDARD_REAL_UNARY(log)
+EIGEN_MATHFUNC_STANDARD_REAL_UNARY(sin)
+EIGEN_MATHFUNC_STANDARD_REAL_UNARY(cos)
+EIGEN_MATHFUNC_STANDARD_REAL_UNARY(tan)
+EIGEN_MATHFUNC_STANDARD_REAL_UNARY(asin)
+EIGEN_MATHFUNC_STANDARD_REAL_UNARY(acos)
 
 /****************************************************************************
 * Implementation of atan2                                                *
@@ -653,7 +503,8 @@ struct atan2_default_impl
   typedef Scalar retval;
   static inline Scalar run(const Scalar& x, const Scalar& y)
   {
-    return std::atan2(x, y);
+    using std::atan2;
+    return atan2(x, y);
   }
 };
 
@@ -692,7 +543,8 @@ struct pow_default_impl
   typedef Scalar retval;
   static inline Scalar run(const Scalar& x, const Scalar& y)
   {
-    return std::pow(x, y);
+    using std::pow;
+    return pow(x, y);
   }
 };
 
@@ -884,7 +736,8 @@ struct scalar_fuzzy_default_impl<Scalar, false, false>
   }
   static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
   {
-    return abs(x - y) <= std::min(abs(x), abs(y)) * prec;
+    using std::min;
+    return abs(x - y) <= min(abs(x), abs(y)) * prec;
   }
   static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar& prec)
   {
@@ -922,7 +775,8 @@ struct scalar_fuzzy_default_impl<Scalar, true, false>
   }
   static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
   {
-    return abs2(x - y) <= std::min(abs2(x), abs2(y)) * prec * prec;
+    using std::min;
+    return abs2(x - y) <= min(abs2(x), abs2(y)) * prec * prec;
   }
 };
 

Eigen/src/Core/Matrix.h

@@ -43,8 +43,8 @@
   * \tparam _Cols Number of columns, or \b Dynamic
   *
   * The remaining template parameters are optional -- in most cases you don't have to worry about them.
-  * \tparam _Options \anchor matrix_tparam_options A combination of either \b RowMajor or \b ColMajor, and of either
-  *                 \b AutoAlign or \b DontAlign.
+  * \tparam _Options \anchor matrix_tparam_options A combination of either \b #RowMajor or \b #ColMajor, and of either
+  *                 \b #AutoAlign or \b #DontAlign.
   *                 The former controls \ref TopicStorageOrders "storage order", and defaults to column-major. The latter controls alignment, which is required
   *                 for vectorization. It defaults to aligning matrices except for fixed sizes that aren't a multiple of the packet size.
   * \tparam _MaxRows Maximum number of rows. Defaults to \a _Rows (\ref maxrows "note").
@@ -79,6 +79,9 @@
   * m(0, 3) = 3;
   * \endcode
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_MATRIX_PLUGIN.
+  *
   * <i><b>Some notes:</b></i>
   *
   * <dl>

Eigen/src/Core/MatrixBase.h

@@ -38,7 +38,7 @@
   * Note that some methods are defined in other modules such as the \ref LU_Module LU module
   * for all functions related to matrix inversions.
   *
-  * \param Derived is the derived type, e.g. a matrix type, or an expression, etc.
+  * \tparam Derived is the derived type, e.g. a matrix type, or an expression, etc.
   *
   * When writing a function taking Eigen objects as argument, if you want your function
   * to take as argument any matrix, vector, or expression, just let it take a
@@ -53,6 +53,9 @@
     }
   * \endcode
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_MATRIXBASE_PLUGIN.
+  *
   * \sa \ref TopicClassHierarchy
   */
 template<typename Derived> class MatrixBase

Eigen/src/Core/PlainObjectBase.h

@@ -42,6 +42,10 @@ template<typename MatrixTypeA, typename MatrixTypeB, bool SwapPointers> struct m
 
 /**
   * \brief %Dense storage base class for matrices and arrays.
+  *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_PLAINOBJECTBASE_PLUGIN.
+  *
   * \sa \ref TopicClassHierarchy
   */
 template<typename Derived>
@@ -283,33 +287,47 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type
       else resize(other.rows(), other.cols());
     }
 
-    /** Resizes \c *this to a \a rows x \a cols matrix while leaving old values of \c *this untouched.
+    /** Resizes the matrix to \a rows x \a cols while leaving old values untouched.
       *
-      * This method is intended for dynamic-size matrices. If you only want to change the number
-      * of rows and/or of columns, you can use conservativeResize(NoChange_t, Index),
+      * The method is intended for matrices of dynamic size. If you only want to change the number
+      * of rows and/or of columns, you can use conservativeResize(NoChange_t, Index) or
       * conservativeResize(Index, NoChange_t).
       *
-      * The top-left part of the resized matrix will be the same as the overlapping top-left corner
-      * of \c *this. In case values need to be appended to the matrix they will be uninitialized.
+      * Matrices are resized relative to the top-left element. In case values need to be 
+      * appended to the matrix they will be uninitialized.
       */
     EIGEN_STRONG_INLINE void conservativeResize(Index rows, Index cols)
     {
       internal::conservative_resize_like_impl<Derived>::run(*this, rows, cols);
     }
 
+    /** Resizes the matrix to \a rows x \a cols while leaving old values untouched.
+      *
+      * As opposed to conservativeResize(Index rows, Index cols), this version leaves
+      * the number of columns unchanged.
+      *
+      * In case the matrix is growing, new rows will be uninitialized.
+      */
     EIGEN_STRONG_INLINE void conservativeResize(Index rows, NoChange_t)
     {
       // Note: see the comment in conservativeResize(Index,Index)
       conservativeResize(rows, cols());
     }
 
+    /** Resizes the matrix to \a rows x \a cols while leaving old values untouched.
+      *
+      * As opposed to conservativeResize(Index rows, Index cols), this version leaves
+      * the number of rows unchanged.
+      *
+      * In case the matrix is growing, new columns will be uninitialized.
+      */
     EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, Index cols)
     {
       // Note: see the comment in conservativeResize(Index,Index)
       conservativeResize(rows(), cols);
     }
 
-    /** Resizes \c *this to a vector of length \a size while retaining old values of *this.
+    /** Resizes the vector to \a size while retaining old values.
       *
       * \only_for_vectors. This method does not work for
       * partially dynamic matrices when the static dimension is anything other
@@ -322,6 +340,15 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type
       internal::conservative_resize_like_impl<Derived>::run(*this, size);
     }
 
+    /** Resizes the matrix to \a rows x \a cols of \c other, while leaving old values untouched.
+      *
+      * The method is intended for matrices of dynamic size. If you only want to change the number
+      * of rows and/or of columns, you can use conservativeResize(NoChange_t, Index) or
+      * conservativeResize(Index, NoChange_t).
+      *
+      * Matrices are resized relative to the top-left element. In case values need to be 
+      * appended to the matrix they will copied from \c other.
+      */
     template<typename OtherDerived>
     EIGEN_STRONG_INLINE void conservativeResizeLike(const DenseBase<OtherDerived>& other)
     {

Eigen/src/Core/Product.h

@@ -375,8 +375,23 @@ struct gemv_static_vector_if<Scalar,Size,Dynamic,true>
 template<typename Scalar,int Size,int MaxSize>
 struct gemv_static_vector_if<Scalar,Size,MaxSize,true>
 {
+  #if EIGEN_ALIGN_STATICALLY
   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data;
   EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; }
+  #else
+  // Some architectures cannot align on the stack,
+  // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
+  enum {
+    ForceAlignment  = internal::packet_traits<Scalar>::Vectorizable,
+    PacketSize      = internal::packet_traits<Scalar>::size
+  };
+  internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data;
+  EIGEN_STRONG_INLINE Scalar* data() {
+    return ForceAlignment
+            ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16)
+            : m_data.array;
+  }
+  #endif
 };
 
 template<> struct gemv_selector<OnTheRight,ColMajor,true>
@@ -411,28 +426,21 @@ template<> struct gemv_selector<OnTheRight,ColMajor,true>
 
     gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
 
-    bool alphaIsCompatible = (!ComplexByReal) || (imag(actualAlpha)==RealScalar(0));
+    // this is written like this (i.e., with a ?:) to workaround an ICE with ICC 12
+    bool alphaIsCompatible = (!ComplexByReal) ? true : (imag(actualAlpha)==RealScalar(0));
     bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
     
     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
 
-    ResScalar* actualDestPtr;
-    bool freeDestPtr = false;
-    if (evalToDest)
-    {
-      actualDestPtr = &dest.coeffRef(0);
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
+                                                  evalToDest ? dest.data() : static_dest.data());
+    
+    if(!evalToDest)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = dest.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualDestPtr = static_dest.data())==0)
-      {
-        freeDestPtr = true;
-        actualDestPtr = ei_aligned_stack_new(ResScalar,dest.size());
-      }
       if(!alphaIsCompatible)
       {
         MappedDest(actualDestPtr, dest.size()).setZero();
@@ -456,7 +464,6 @@ template<> struct gemv_selector<OnTheRight,ColMajor,true>
         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
       else
         dest = MappedDest(actualDestPtr, dest.size());
-      if(freeDestPtr) ei_aligned_stack_delete(ResScalar, actualDestPtr, dest.size());
     }
   }
 };
@@ -490,23 +497,15 @@ template<> struct gemv_selector<OnTheRight,RowMajor,true>
 
     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
 
-    RhsScalar* actualRhsPtr;
-    bool freeRhsPtr = false;
-    if (DirectlyUseRhs)
-    {
-      actualRhsPtr = const_cast<RhsScalar*>(&actualRhs.coeffRef(0));
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
+        DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
+
+    if(!DirectlyUseRhs)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = actualRhs.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualRhsPtr = static_rhs.data())==0)
-      {
-        freeRhsPtr = true;
-        actualRhsPtr = ei_aligned_stack_new(RhsScalar, actualRhs.size());
-      }
       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
     }
 
@@ -517,8 +516,6 @@ template<> struct gemv_selector<OnTheRight,RowMajor,true>
         actualRhsPtr, 1,
         &dest.coeffRef(0,0), dest.innerStride(),
         actualAlpha);
-
-    if((!DirectlyUseRhs) && freeRhsPtr) ei_aligned_stack_delete(RhsScalar, actualRhsPtr, prod.rhs().size());
   }
 };
 

Eigen/src/Core/Redux.h

@@ -216,7 +216,7 @@ struct redux_impl<Func, Derived, LinearVectorizedTraversal, NoUnrolling>
     const Index packetSize = packet_traits<Scalar>::size;
     const Index alignedStart = first_aligned(mat);
     enum {
-      alignment = (Derived::Flags & DirectAccessBit) || (Derived::Flags & AlignedBit)
+      alignment = bool(Derived::Flags & DirectAccessBit) || bool(Derived::Flags & AlignedBit)
                 ? Aligned : Unaligned
     };
     const Index alignedSize = ((size-alignedStart)/packetSize)*packetSize;

Eigen/src/Core/ReturnByValue.h

@@ -71,9 +71,9 @@ template<typename Derived> class ReturnByValue
 
     template<typename Dest>
     inline void evalTo(Dest& dst) const
-    { static_cast<const Derived* const>(this)->evalTo(dst); }
-    inline Index rows() const { return static_cast<const Derived* const>(this)->rows(); }
-    inline Index cols() const { return static_cast<const Derived* const>(this)->cols(); }
+    { static_cast<const Derived*>(this)->evalTo(dst); }
+    inline Index rows() const { return static_cast<const Derived*>(this)->rows(); }
+    inline Index cols() const { return static_cast<const Derived*>(this)->cols(); }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
 #define Unusable YOU_ARE_TRYING_TO_ACCESS_A_SINGLE_COEFFICIENT_IN_A_SPECIAL_EXPRESSION_WHERE_THAT_IS_NOT_ALLOWED_BECAUSE_THAT_WOULD_BE_INEFFICIENT

Eigen/src/Core/SelfAdjointView.h

@@ -32,13 +32,13 @@
   * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix
   *
   * \param MatrixType the type of the dense matrix storing the coefficients
-  * \param TriangularPart can be either \c Lower or \c Upper
+  * \param TriangularPart can be either \c #Lower or \c #Upper
   *
   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
   * and most of the time this is the only way that it is used.
   *
-  * \sa class TriangularBase, MatrixBase::selfAdjointView()
+  * \sa class TriangularBase, MatrixBase::selfadjointView()
   */
 
 namespace internal {

Eigen/src/Core/SolveTriangular.h

@@ -74,26 +74,19 @@ struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,1>
     // FIXME find a way to allow an inner stride if packet_traits<Scalar>::size==1
 
     bool useRhsDirectly = Rhs::InnerStrideAtCompileTime==1 || rhs.innerStride()==1;
-    RhsScalar* actualRhs;
-    if(useRhsDirectly)
-    {
-      actualRhs = &rhs.coeffRef(0);
-    }
-    else
-    {
-      actualRhs = ei_aligned_stack_new(RhsScalar,rhs.size());
+
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhs,rhs.size(),
+                                                  (useRhsDirectly ? rhs.data() : 0));
+                                                  
+    if(!useRhsDirectly)
       MappedRhs(actualRhs,rhs.size()) = rhs;
-    }
 
     triangular_solve_vector<LhsScalar, RhsScalar, typename Lhs::Index, Side, Mode, LhsProductTraits::NeedToConjugate,
                             (int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor>
       ::run(actualLhs.cols(), actualLhs.data(), actualLhs.outerStride(), actualRhs);
 
     if(!useRhsDirectly)
-    {
       rhs = MappedRhs(actualRhs, rhs.size());
-      ei_aligned_stack_delete(RhsScalar, actualRhs, rhs.size());
-    }
   }
 };
 

Eigen/src/Core/StableNorm.h

@@ -56,6 +56,7 @@ template<typename Derived>
 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
 MatrixBase<Derived>::stableNorm() const
 {
+  using std::min;
   const Index blockSize = 4096;
   RealScalar scale = 0;
   RealScalar invScale = 1;
@@ -68,7 +69,7 @@ MatrixBase<Derived>::stableNorm() const
   if (bi>0)
     internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
   for (; bi<n; bi+=blockSize)
-    internal::stable_norm_kernel(this->segment(bi,std::min(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
+    internal::stable_norm_kernel(this->segment(bi,min(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
   return scale * internal::sqrt(ssq);
 }
 
@@ -85,6 +86,9 @@ template<typename Derived>
 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
 MatrixBase<Derived>::blueNorm() const
 {
+  using std::pow;
+  using std::min;
+  using std::max;
   static Index nmax = -1;
   static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
   if(nmax <= 0)
@@ -107,17 +111,17 @@ MatrixBase<Derived>::blueNorm() const
     rbig  = std::numeric_limits<RealScalar>::max();         // largest floating-point number
 
     iexp  = -((1-iemin)/2);
-    b1    = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));  // lower boundary of midrange
+    b1    = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));  // lower boundary of midrange
     iexp  = (iemax + 1 - it)/2;
-    b2    = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));   // upper boundary of midrange
+    b2    = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));   // upper boundary of midrange
 
     iexp  = (2-iemin)/2;
-    s1m   = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for lower range
+    s1m   = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for lower range
     iexp  = - ((iemax+it)/2);
-    s2m   = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for upper range
+    s2m   = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for upper range
 
     overfl  = rbig*s2m;             // overflow boundary for abig
-    eps     = RealScalar(std::pow(double(ibeta), 1-it));
+    eps     = RealScalar(pow(double(ibeta), 1-it));
     relerr  = internal::sqrt(eps);         // tolerance for neglecting asml
     abig    = RealScalar(1.0/eps - 1.0);
     if (RealScalar(nbig)>abig)  nmax = int(abig);  // largest safe n
@@ -163,8 +167,8 @@ MatrixBase<Derived>::blueNorm() const
   }
   else
     return internal::sqrt(amed);
-  asml = std::min(abig, amed);
-  abig = std::max(abig, amed);
+  asml = min(abig, amed);
+  abig = max(abig, amed);
   if(asml <= abig*relerr)
     return abig;
   else

Eigen/src/Core/TriangularMatrix.h

@@ -134,13 +134,13 @@ template<typename Derived> class TriangularBase : public EigenBase<Derived>
   * \brief Base class for triangular part in a matrix
   *
   * \param MatrixType the type of the object in which we are taking the triangular part
-  * \param Mode the kind of triangular matrix expression to construct. Can be Upper,
-  *             Lower, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field;
-  *             it must have either Upper or Lower, and additionnaly it may have either
-  *             UnitDiag or Selfadjoint.
+  * \param Mode the kind of triangular matrix expression to construct. Can be #Upper,
+  *             #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower.
+  *             This is in fact a bit field; it must have either #Upper or #Lower, 
+  *             and additionnaly it may have #UnitDiag or #ZeroDiag or neither.
   *
   * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular
-  * matrices one should speak ok "trapezoid" parts. This class is the return type
+  * matrices one should speak of "trapezoid" parts. This class is the return type
   * of MatrixBase::triangularView() and most of the time this is the only way it is used.
   *
   * \sa MatrixBase::triangularView()
@@ -448,6 +448,8 @@ struct triangular_assignment_selector
     col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
     row = (UnrollCount-1) % Derived1::RowsAtCompileTime
   };
+  
+  typedef typename Derived1::Scalar Scalar;
 
   inline static void run(Derived1 &dst, const Derived2 &src)
   {
@@ -466,9 +468,9 @@ struct triangular_assignment_selector
     else if(ClearOpposite)
     {
       if (Mode&UnitDiag && row==col)
-        dst.coeffRef(row, col) = 1;
+        dst.coeffRef(row, col) = Scalar(1);
       else
-        dst.coeffRef(row, col) = 0;
+        dst.coeffRef(row, col) = Scalar(0);
     }
   }
 };
@@ -484,6 +486,7 @@ template<typename Derived1, typename Derived2, bool ClearOpposite>
 struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite>
 {
   typedef typename Derived1::Index Index;
+  typedef typename Derived1::Scalar Scalar;
   inline static void run(Derived1 &dst, const Derived2 &src)
   {
     for(Index j = 0; j < dst.cols(); ++j)
@@ -493,7 +496,7 @@ struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearO
         dst.copyCoeff(i, j, src);
       if (ClearOpposite)
         for(Index i = maxi+1; i < dst.rows(); ++i)
-          dst.coeffRef(i, j) = 0;
+          dst.coeffRef(i, j) = Scalar(0);
     }
   }
 };
@@ -511,7 +514,7 @@ struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearO
       Index maxi = std::min(j, dst.rows());
       if (ClearOpposite)
         for(Index i = 0; i < maxi; ++i)
-          dst.coeffRef(i, j) = 0;
+          dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
     }
   }
 };
@@ -547,7 +550,7 @@ struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic
       Index maxi = std::min(j, dst.rows()-1);
       if (ClearOpposite)
         for(Index i = 0; i <= maxi; ++i)
-          dst.coeffRef(i, j) = 0;
+          dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
     }
   }
 };
@@ -756,8 +759,8 @@ typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Deriv
 /**
   * \returns an expression of a triangular view extracted from the current matrix
   *
-  * The parameter \a Mode can have the following values: \c Upper, \c StrictlyUpper, \c UnitUpper,
-  * \c Lower, \c StrictlyLower, \c UnitLower.
+  * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
+  * \c #Lower, \c #StrictlyLower, \c #UnitLower.
   *
   * Example: \include MatrixBase_extract.cpp
   * Output: \verbinclude MatrixBase_extract.out

Eigen/src/Core/VectorwiseOp.h

@@ -31,9 +31,9 @@
   *
   * \brief Generic expression of a partially reduxed matrix
   *
-  * \param MatrixType the type of the matrix we are applying the redux operation
-  * \param MemberOp type of the member functor
-  * \param Direction indicates the direction of the redux (Vertical or Horizontal)
+  * \tparam MatrixType the type of the matrix we are applying the redux operation
+  * \tparam MemberOp type of the member functor
+  * \tparam Direction indicates the direction of the redux (#Vertical or #Horizontal)
   *
   * This class represents an expression of a partial redux operator of a matrix.
   * It is the return type of some VectorwiseOp functions,
@@ -164,7 +164,7 @@ struct member_redux {
   * \brief Pseudo expression providing partial reduction operations
   *
   * \param ExpressionType the type of the object on which to do partial reductions
-  * \param Direction indicates the direction of the redux (Vertical or Horizontal)
+  * \param Direction indicates the direction of the redux (#Vertical or #Horizontal)
   *
   * This class represents a pseudo expression with partial reduction features.
   * It is the return type of DenseBase::colwise() and DenseBase::rowwise()

Eigen/src/Core/arch/AltiVec/Complex.h

@@ -48,6 +48,7 @@ template<> struct packet_traits<std::complex<float> >  : default_packet_traits
   typedef Packet2cf type;
   enum {
     Vectorizable = 1,
+    AlignedOnScalar = 1,
     size = 2,
 
     HasAdd    = 1,
@@ -69,19 +70,17 @@ template<> EIGEN_STRONG_INLINE Packet2cf pset1<Packet2cf>(const std::complex<flo
 {
   Packet2cf res;
   /* On AltiVec we cannot load 64-bit registers, so wa have to take care of alignment */
-  if ((ptrdiff_t)&from % 16 == 0) {
-    res.v = pload((const float *)&from);
-    res.v = vec_perm(res.v, res.v, p16uc_PSET_HI);
-  } else {
-    res.v = ploadu((const float *)&from);
-    res.v = vec_perm(res.v, res.v, p16uc_PSET_LO);
-  }
+  if((ptrdiff_t(&from) % 16) == 0)
+    res.v = pload<Packet4f>((const float *)&from);
+  else
+    res.v = ploadu<Packet4f>((const float *)&from);
+  res.v = vec_perm(res.v, res.v, p16uc_PSET_HI);
   return res;
 }
 
 template<> EIGEN_STRONG_INLINE Packet2cf padd<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(vec_add(a.v,b.v)); }
 template<> EIGEN_STRONG_INLINE Packet2cf psub<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(vec_sub(a.v,b.v)); }
-template<> EIGEN_STRONG_INLINE Packet2cf pnegate(const Packet2cf& a) { return Packet2cf(psub<Packet4f>(p4f_ZERO, a.v)); }
+template<> EIGEN_STRONG_INLINE Packet2cf pnegate(const Packet2cf& a) { return Packet2cf(pnegate(a.v)); }
 template<> EIGEN_STRONG_INLINE Packet2cf pconj(const Packet2cf& a) { return Packet2cf((Packet4f)vec_xor((Packet4ui)a.v, p4ui_CONJ_XOR)); }
 
 template<> EIGEN_STRONG_INLINE Packet2cf pmul<Packet2cf>(const Packet2cf& a, const Packet2cf& b)
@@ -108,8 +107,13 @@ template<> EIGEN_STRONG_INLINE Packet2cf por    <Packet2cf>(const Packet2cf& a,
 template<> EIGEN_STRONG_INLINE Packet2cf pxor   <Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(vec_xor(a.v,b.v)); }
 template<> EIGEN_STRONG_INLINE Packet2cf pandnot<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(vec_and(a.v, vec_nor(b.v,b.v))); }
 
-template<> EIGEN_STRONG_INLINE Packet2cf pload <std::complex<float> >(const std::complex<float>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload((const float*)from)); }
-template<> EIGEN_STRONG_INLINE Packet2cf ploadu<std::complex<float> >(const std::complex<float>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu((const float*)from)); }
+template<> EIGEN_STRONG_INLINE Packet2cf pload <Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload<Packet4f>((const float*)from)); }
+template<> EIGEN_STRONG_INLINE Packet2cf ploadu<Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu<Packet4f>((const float*)from)); }
+
+template<> EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex<float>*     from)
+{
+  return pset1<Packet2cf>(*from);
+}
 
 template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((float*)to, from.v); }
 template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((float*)to, from.v); }
@@ -136,7 +140,7 @@ template<> EIGEN_STRONG_INLINE std::complex<float> predux<Packet2cf>(const Packe
   Packet4f b;
   b = (Packet4f) vec_sld(a.v, a.v, 8);
   b = padd(a.v, b);
-  return pfirst(Packet2cf(sum));
+  return pfirst(Packet2cf(b));
 }
 
 template<> EIGEN_STRONG_INLINE Packet2cf preduxp<Packet2cf>(const Packet2cf* vecs)
@@ -180,7 +184,7 @@ template<> struct conj_helper<Packet2cf, Packet2cf, false,true>
 
   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
   {
-    return pmul(a, pconj(b));
+    return internal::pmul(a, pconj(b));
   }
 };
 
@@ -191,7 +195,7 @@ template<> struct conj_helper<Packet2cf, Packet2cf, true,false>
 
   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
   {
-    return pmul(pconj(a), b);
+    return internal::pmul(pconj(a), b);
   }
 };
 
@@ -202,7 +206,7 @@ template<> struct conj_helper<Packet2cf, Packet2cf, true,true>
 
   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
   {
-    return pconj(pmul(a, b));
+    return pconj(internal::pmul(a, b));
   }
 };
 
@@ -214,6 +218,11 @@ template<> EIGEN_STRONG_INLINE Packet2cf pdiv<Packet2cf>(const Packet2cf& a, con
   return Packet2cf(pdiv(res.v, vec_add(s,vec_perm(s, s, p16uc_COMPLEX_REV))));
 }
 
+template<> EIGEN_STRONG_INLINE Packet2cf pcplxflip<Packet2cf>(const Packet2cf& x)
+{
+  return Packet2cf(vec_perm(x.v, x.v, p16uc_COMPLEX_REV));
+}
+
 } // end namespace internal
 
 #endif // EIGEN_COMPLEX_ALTIVEC_H

Eigen/src/Core/arch/AltiVec/PacketMath.h

@@ -73,6 +73,7 @@ static Packet4f p4f_COUNTDOWN = { 3.0, 2.0, 1.0, 0.0 };
 static Packet4i p4i_COUNTDOWN = { 3, 2, 1, 0 };
 static Packet16uc p16uc_REVERSE = {12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3};
 static Packet16uc p16uc_FORWARD = vec_lvsl(0, (float*)0);
+static Packet16uc p16uc_DUPLICATE = {0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7};
 
 static _EIGEN_DECLARE_CONST_FAST_Packet4f(ZERO, 0);
 static _EIGEN_DECLARE_CONST_FAST_Packet4i(ZERO, 0);
@@ -292,6 +293,21 @@ template<> EIGEN_STRONG_INLINE Packet4i ploadu<Packet4i>(const int* from)
   return (Packet4i) vec_perm(MSQ, LSQ, mask);    // align the data
 }
 
+template<> EIGEN_STRONG_INLINE Packet4f ploaddup<Packet4f>(const float*   from)
+{
+  Packet4f p;
+  if((ptrdiff_t(&from) % 16) == 0)  p = pload<Packet4f>(from);
+  else                              p = ploadu<Packet4f>(from);
+  return vec_perm(p, p, p16uc_DUPLICATE);
+}
+template<> EIGEN_STRONG_INLINE Packet4i ploaddup<Packet4i>(const int*     from)
+{
+  Packet4i p;
+  if((ptrdiff_t(&from) % 16) == 0)  p = pload<Packet4i>(from);
+  else                              p = ploadu<Packet4i>(from);
+  return vec_perm(p, p, p16uc_DUPLICATE);
+}
+
 template<> EIGEN_STRONG_INLINE void pstore<float>(float*   to, const Packet4f& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); }
 template<> EIGEN_STRONG_INLINE void pstore<int>(int*       to, const Packet4i& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); }
 

Eigen/src/Core/arch/NEON/Complex.h

@@ -43,6 +43,7 @@ template<> struct packet_traits<std::complex<float> >  : default_packet_traits
   typedef Packet2cf type;
   enum {
     Vectorizable = 1,
+    AlignedOnScalar = 1,
     size = 2,
 
     HasAdd    = 1,
@@ -120,6 +121,8 @@ template<> EIGEN_STRONG_INLINE Packet2cf pandnot<Packet2cf>(const Packet2cf& a,
 template<> EIGEN_STRONG_INLINE Packet2cf pload<Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload<Packet4f>((const float*)from)); }
 template<> EIGEN_STRONG_INLINE Packet2cf ploadu<Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu<Packet4f>((const float*)from)); }
 
+template<> EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex<float>* from) { return pset1<Packet2cf>(*from); }
+
 template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((float*)to, from.v); }
 template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((float*)to, from.v); }
 
@@ -144,7 +147,7 @@ template<> EIGEN_STRONG_INLINE Packet2cf preverse(const Packet2cf& a)
   return Packet2cf(a_r128);
 }
 
-EIGEN_STRONG_INLINE Packet2cf pcplxflip(const Packet2cf& a)
+template<> EIGEN_STRONG_INLINE Packet2cf pcplxflip<Packet2cf>(const Packet2cf& a)
 {
   return Packet2cf(vrev64q_f32(a.v));
 }
@@ -220,7 +223,7 @@ template<> struct conj_helper<Packet2cf, Packet2cf, false,true>
 
   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
   {
-    return pmul(a, pconj(b));
+    return internal::pmul(a, pconj(b));
   }
 };
 
@@ -231,7 +234,7 @@ template<> struct conj_helper<Packet2cf, Packet2cf, true,false>
 
   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
   {
-    return pmul(pconj(a), b);
+    return internal::pmul(pconj(a), b);
   }
 };
 
@@ -242,7 +245,7 @@ template<> struct conj_helper<Packet2cf, Packet2cf, true,true>
 
   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
   {
-    return pconj(pmul(a, b));
+    return pconj(internal::pmul(a, b));
   }
 };
 

Eigen/src/Core/arch/NEON/PacketMath.h

@@ -41,7 +41,7 @@ namespace internal {
 
 typedef float32x4_t Packet4f;
 typedef int32x4_t   Packet4i;
-typedef uint32x4_t   Packet4ui;
+typedef uint32x4_t  Packet4ui;
 
 #define _EIGEN_DECLARE_CONST_Packet4f(NAME,X) \
   const Packet4f p4f_##NAME = pset1<Packet4f>(X)
@@ -84,8 +84,8 @@ template<> struct packet_traits<int>    : default_packet_traits
   };
 };
 
-#if (defined __GNUC__) && (!(EIGEN_GNUC_AT_LEAST(4,4)))
-// workaround gcc 4.2 and 4.3 compilatin issue
+#if EIGEN_GNUC_AT_MOST(4,4)
+// workaround gcc 4.2, 4.3 and 4.4 compilatin issue
 EIGEN_STRONG_INLINE float32x4_t vld1q_f32(const float* x) { return ::vld1q_f32((const float32_t*)x); }
 EIGEN_STRONG_INLINE float32x2_t vld1_f32 (const float* x) { return ::vld1_f32 ((const float32_t*)x); }
 EIGEN_STRONG_INLINE void        vst1q_f32(float* to, float32x4_t from) { ::vst1q_f32((float32_t*)to,from); }
@@ -100,12 +100,12 @@ template<> EIGEN_STRONG_INLINE Packet4i pset1<Packet4i>(const int&    from)   {
 
 template<> EIGEN_STRONG_INLINE Packet4f plset<float>(const float& a)
 {
-  Packet4f countdown = { 3, 2, 1, 0 };
+  Packet4f countdown = { 0, 1, 2, 3 };
   return vaddq_f32(pset1<Packet4f>(a), countdown);
 }
 template<> EIGEN_STRONG_INLINE Packet4i plset<int>(const int& a)
 {
-  Packet4i countdown = { 3, 2, 1, 0 };
+  Packet4i countdown = { 0, 1, 2, 3 };
   return vaddq_s32(pset1<Packet4i>(a), countdown);
 }
 
@@ -191,14 +191,14 @@ template<> EIGEN_STRONG_INLINE Packet4f ploaddup<Packet4f>(const float*   from)
 {
   float32x2_t lo, hi;
   lo = vdup_n_f32(*from);
-  hi = vdup_n_f32(*from);
+  hi = vdup_n_f32(*(from+1));
   return vcombine_f32(lo, hi);
 }
 template<> EIGEN_STRONG_INLINE Packet4i ploaddup<Packet4i>(const int*     from)
 {
   int32x2_t lo, hi;
   lo = vdup_n_s32(*from);
-  hi = vdup_n_s32(*from);
+  hi = vdup_n_s32(*(from+1));
   return vcombine_s32(lo, hi);
 }
 

Eigen/src/Core/arch/SSE/Complex.h

@@ -97,23 +97,30 @@ template<> EIGEN_STRONG_INLINE Packet2cf pandnot<Packet2cf>(const Packet2cf& a,
 template<> EIGEN_STRONG_INLINE Packet2cf pload <Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload<Packet4f>(&real_ref(*from))); }
 template<> EIGEN_STRONG_INLINE Packet2cf ploadu<Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu<Packet4f>(&real_ref(*from))); }
 
+template<> EIGEN_STRONG_INLINE Packet2cf pset1<Packet2cf>(const std::complex<float>&  from)
+{
+  Packet2cf res;
+  #if EIGEN_GNUC_AT_MOST(4,2)
+  // workaround annoying "may be used uninitialized in this function" warning with gcc 4.2
+  res.v = _mm_loadl_pi(_mm_set1_ps(0.0f), (const __m64*)&from);
+  #else
+  res.v = _mm_loadl_pi(res.v, (const __m64*)&from);
+  #endif
+  return Packet2cf(_mm_movelh_ps(res.v,res.v));
+}
+
+template<> EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex<float>* from) { return pset1<Packet2cf>(*from); }
+
 template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore(&real_ref(*to), from.v); }
 template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu(&real_ref(*to), from.v); }
 
 template<> EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float> *   addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }
 
-template<> EIGEN_STRONG_INLINE Packet2cf pset1<Packet2cf>(const std::complex<float>&  from)
-{
-  Packet2cf res;
-  res.v = _mm_loadl_pi(res.v, (const __m64*)&from);
-  return Packet2cf(_mm_movelh_ps(res.v,res.v));
-}
-
 template<> EIGEN_STRONG_INLINE std::complex<float>  pfirst<Packet2cf>(const Packet2cf& a)
 {
-  #if (defined __GNUC__) && (__GNUC__==4) && (__GNUC_MINOR__==2) && (__GNUC_PATCHLEVEL__<=3)
-  // workaround gcc 4.2.1 ICE (mac's gcc version) - I'm not sure how the 4.2.2 and 4.2.3 deal with it, but 4.2.4 works well.
-  // this is not performance wise ideal, but who cares...
+  #if EIGEN_GNUC_AT_MOST(4,3)
+  // Workaround gcc 4.2 ICE - this is not performance wise ideal, but who cares...
+  // This workaround also fix invalid code generation with gcc 4.3
   EIGEN_ALIGN16 std::complex<float> res[2];
   _mm_store_ps((float*)res, a.v);
   return res[0];
@@ -308,6 +315,8 @@ template<> EIGEN_STRONG_INLINE Packet1cd ploadu<Packet1cd>(const std::complex<do
 template<> EIGEN_STRONG_INLINE Packet1cd pset1<Packet1cd>(const std::complex<double>&  from)
 { /* here we really have to use unaligned loads :( */ return ploadu<Packet1cd>(&from); }
 
+template<> EIGEN_STRONG_INLINE Packet1cd ploaddup<Packet1cd>(const std::complex<double>* from) { return pset1<Packet1cd>(*from); }
+
 // FIXME force unaligned store, this is a temporary fix
 template<> EIGEN_STRONG_INLINE void pstore <std::complex<double> >(std::complex<double> *   to, const Packet1cd& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((double*)to, from.v); }
 template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<double> >(std::complex<double> *   to, const Packet1cd& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v); }

Eigen/src/Core/arch/SSE/PacketMath.h

@@ -110,22 +110,8 @@ template<> struct unpacket_traits<Packet4f> { typedef float  type; enum {size=4}
 template<> struct unpacket_traits<Packet2d> { typedef double type; enum {size=2}; };
 template<> struct unpacket_traits<Packet4i> { typedef int    type; enum {size=4}; };
 
-#ifdef __GNUC__
-// Sometimes GCC implements _mm_set1_p* using multiple moves,
-// that is inefficient :( (e.g., see gemm_pack_rhs)
-template<> EIGEN_STRONG_INLINE Packet4f pset1<Packet4f>(const float&  from) {
-  Packet4f res = _mm_set_ss(from);
-  return vec4f_swizzle1(res,0,0,0,0);
-}
-template<> EIGEN_STRONG_INLINE Packet2d pset1<Packet2d>(const double&  from) {
-  // NOTE the SSE3 intrinsic _mm_loaddup_pd is never faster but sometimes much slower
-  Packet2d res = _mm_set_sd(from);
-  return vec2d_swizzle1(res, 0, 0);
-}
-#else
 template<> EIGEN_STRONG_INLINE Packet4f pset1<Packet4f>(const float&  from) { return _mm_set1_ps(from); }
 template<> EIGEN_STRONG_INLINE Packet2d pset1<Packet2d>(const double& from) { return _mm_set1_pd(from); }
-#endif
 template<> EIGEN_STRONG_INLINE Packet4i pset1<Packet4i>(const int&    from) { return _mm_set1_epi32(from); }
 
 template<> EIGEN_STRONG_INLINE Packet4f plset<float>(const float& a) { return _mm_add_ps(pset1<Packet4f>(a), _mm_set_ps(3,2,1,0)); }
@@ -245,29 +231,52 @@ template<> EIGEN_STRONG_INLINE Packet4i pload<Packet4i>(const int*     from) { E
 // a correct instruction dependency.
 // TODO: do the same for MSVC (ICC is compatible)
 // NOTE: with the code below, MSVC's compiler crashes!
+
+#if defined(__GNUC__) && defined(__i386__)
+  // bug 195: gcc/i386 emits weird x87 fldl/fstpl instructions for _mm_load_sd
+  #define EIGEN_AVOID_CUSTOM_UNALIGNED_LOADS 1
+#elif defined(__clang__)
+  // bug 201: Segfaults in __mm_loadh_pd with clang 2.8
+  #define EIGEN_AVOID_CUSTOM_UNALIGNED_LOADS 1
+#else
+  #define EIGEN_AVOID_CUSTOM_UNALIGNED_LOADS 0
+#endif
+
 template<> EIGEN_STRONG_INLINE Packet4f ploadu<Packet4f>(const float* from)
 {
   EIGEN_DEBUG_UNALIGNED_LOAD
+#if EIGEN_AVOID_CUSTOM_UNALIGNED_LOADS
+  return _mm_loadu_ps(from);
+#else
   __m128d res;
   res =  _mm_load_sd((const double*)(from)) ;
   res =  _mm_loadh_pd(res, (const double*)(from+2)) ;
   return _mm_castpd_ps(res);
+#endif
 }
 template<> EIGEN_STRONG_INLINE Packet2d ploadu<Packet2d>(const double* from)
 {
   EIGEN_DEBUG_UNALIGNED_LOAD
+#if EIGEN_AVOID_CUSTOM_UNALIGNED_LOADS
+  return _mm_loadu_pd(from);
+#else
   __m128d res;
   res = _mm_load_sd(from) ;
   res = _mm_loadh_pd(res,from+1);
   return res;
+#endif
 }
 template<> EIGEN_STRONG_INLINE Packet4i ploadu<Packet4i>(const int* from)
 {
   EIGEN_DEBUG_UNALIGNED_LOAD
+#if EIGEN_AVOID_CUSTOM_UNALIGNED_LOADS
+  return _mm_loadu_si128(reinterpret_cast<const Packet4i*>(from));
+#else
   __m128d res;
   res =  _mm_load_sd((const double*)(from)) ;
   res =  _mm_loadh_pd(res, (const double*)(from+2)) ;
   return _mm_castpd_si128(res);
+#endif
 }
 #endif
 
@@ -296,6 +305,19 @@ template<> EIGEN_STRONG_INLINE void pstoreu<double>(double* to, const Packet2d&
 template<> EIGEN_STRONG_INLINE void pstoreu<float>(float*  to, const Packet4f& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, _mm_castps_pd(from)); }
 template<> EIGEN_STRONG_INLINE void pstoreu<int>(int*      to, const Packet4i& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, _mm_castsi128_pd(from)); }
 
+// some compilers might be tempted to perform multiple moves instead of using a vector path.
+template<> EIGEN_STRONG_INLINE void pstore1<Packet4f>(float* to, const float& a)
+{
+  Packet4f pa = _mm_set_ss(a);
+  pstore(to, vec4f_swizzle1(pa,0,0,0,0));
+}
+// some compilers might be tempted to perform multiple moves instead of using a vector path.
+template<> EIGEN_STRONG_INLINE void pstore1<Packet2d>(double* to, const double& a)
+{
+  Packet2d pa = _mm_set_sd(a);
+  pstore(to, vec2d_swizzle1(pa,0,0));
+}
+
 template<> EIGEN_STRONG_INLINE void prefetch<float>(const float*   addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }
 template<> EIGEN_STRONG_INLINE void prefetch<double>(const double* addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }
 template<> EIGEN_STRONG_INLINE void prefetch<int>(const int*       addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }

Eigen/src/Core/products/GeneralBlockPanelKernel.h

@@ -199,7 +199,7 @@ public:
   EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
   {
     for(DenseIndex k=0; k<n; k++)
-      pstore(&b[k*RhsPacketSize], pset1<RhsPacket>(rhs[k]));
+      pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
   }
 
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
@@ -270,7 +270,7 @@ public:
   EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
   {
     for(DenseIndex k=0; k<n; k++)
-      pstore(&b[k*RhsPacketSize], pset1<RhsPacket>(rhs[k]));
+      pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
   }
 
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
@@ -363,8 +363,8 @@ public:
     {
       if(Vectorizable)
       {
-        pstore((RealScalar*)&b[k*ResPacketSize*2+0], pset1<RealPacket>(real(rhs[k])));
-        pstore((RealScalar*)&b[k*ResPacketSize*2+ResPacketSize], pset1<RealPacket>(imag(rhs[k])));
+        pstore1<RealPacket>((RealScalar*)&b[k*ResPacketSize*2+0],             real(rhs[k]));
+        pstore1<RealPacket>((RealScalar*)&b[k*ResPacketSize*2+ResPacketSize], imag(rhs[k]));
       }
       else
         b[k] = rhs[k];
@@ -475,7 +475,7 @@ public:
   EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
   {
     for(DenseIndex k=0; k<n; k++)
-      pstore(&b[k*RhsPacketSize], pset1<RhsPacket>(rhs[k]));
+      pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
   }
 
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
@@ -1009,12 +1009,7 @@ EIGEN_ASM_COMMENT("mybegin4");
     for(Index j2=packet_cols; j2<cols; j2++)
     {
       // unpack B
-      {
-        traits.unpackRhs(depth, &blockB[j2*strideB+offsetB], unpackedB);
-//         const RhsScalar* blB = &blockB[j2*strideB+offsetB];
-//         for(Index k=0; k<depth; k++)
-//           pstore(&unpackedB[k*RhsPacketSize], pset1<RhsPacket>(blB[k]));
-      }
+      traits.unpackRhs(depth, &blockB[j2*strideB+offsetB], unpackedB);
 
       for(Index i=0; i<peeled_mc; i+=mr)
       {

Eigen/src/Core/products/GeneralMatrixMatrix.h

@@ -94,8 +94,9 @@ static void run(Index rows, Index cols, Index depth,
     
     std::size_t sizeA = kc*mc;
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
-    LhsScalar* blockA = ei_aligned_stack_new(LhsScalar, sizeA);
-    RhsScalar* w = ei_aligned_stack_new(RhsScalar, sizeW);
+    ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, 0);
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, w, sizeW, 0);
+    
     RhsScalar* blockB = blocking.blockB();
     eigen_internal_assert(blockB!=0);
 
@@ -154,9 +155,6 @@ static void run(Index rows, Index cols, Index depth,
         #pragma omp atomic
         --(info[j].users);
     }
-
-    ei_aligned_stack_delete(LhsScalar, blockA, kc*mc);
-    ei_aligned_stack_delete(RhsScalar, w, sizeW);
   }
   else
 #endif // EIGEN_HAS_OPENMP
@@ -167,9 +165,10 @@ static void run(Index rows, Index cols, Index depth,
     std::size_t sizeA = kc*mc;
     std::size_t sizeB = kc*cols;
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
-    LhsScalar *blockA = blocking.blockA()==0 ? ei_aligned_stack_new(LhsScalar, sizeA) : blocking.blockA();
-    RhsScalar *blockB = blocking.blockB()==0 ? ei_aligned_stack_new(RhsScalar, sizeB) : blocking.blockB();
-    RhsScalar *blockW = blocking.blockW()==0 ? ei_aligned_stack_new(RhsScalar, sizeW) : blocking.blockW();
+
+    ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, blocking.blockA());
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, blockB, sizeB, blocking.blockB());
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, blockW, sizeW, blocking.blockW());
 
     // For each horizontal panel of the rhs, and corresponding panel of the lhs...
     // (==GEMM_VAR1)
@@ -200,10 +199,6 @@ static void run(Index rows, Index cols, Index depth,
 
       }
     }
-
-    if(blocking.blockA()==0) ei_aligned_stack_delete(LhsScalar, blockA, sizeA);
-    if(blocking.blockB()==0) ei_aligned_stack_delete(RhsScalar, blockB, sizeB);
-    if(blocking.blockW()==0) ei_aligned_stack_delete(RhsScalar, blockW, sizeW);
   }
 }
 

Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h

@@ -83,10 +83,10 @@ struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,
     if(mc > Traits::nr)
       mc = (mc/Traits::nr)*Traits::nr;
 
-    LhsScalar* blockA = ei_aligned_stack_new(LhsScalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*size;
-    RhsScalar* allocatedBlockB = ei_aligned_stack_new(RhsScalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, allocatedBlockB, sizeB, 0);
     RhsScalar* blockB = allocatedBlockB + sizeW;
     
     gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
@@ -125,8 +125,6 @@ struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,
         }
       }
     }
-    ei_aligned_stack_delete(LhsScalar, blockA, kc*mc);
-    ei_aligned_stack_delete(RhsScalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/Parallelizer.h

@@ -30,7 +30,7 @@ namespace internal {
 /** \internal */
 inline void manage_multi_threading(Action action, int* v)
 {
-  static int m_maxThreads = -1;
+  static EIGEN_UNUSED int m_maxThreads = -1;
 
   if(action==SetAction)
   {

Eigen/src/Core/products/SelfadjointMatrixMatrix.h

@@ -263,10 +263,10 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
     // kc must smaller than mc
     kc = std::min(kc,mc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
@@ -313,9 +313,6 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
         gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 
@@ -343,11 +340,10 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLh
     Index mc = rows;  // cache block size along the M direction
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar>(kc, mc, nc);
-
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
@@ -369,9 +365,6 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLh
         gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/SelfadjointMatrixVector.h

@@ -62,14 +62,12 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
   // FIXME this copy is now handled outside product_selfadjoint_vector, so it could probably be removed.
   // if the rhs is not sequentially stored in memory we copy it to a temporary buffer,
   // this is because we need to extract packets
-  const Scalar* EIGEN_RESTRICT rhs = _rhs;
+  ei_declare_aligned_stack_constructed_variable(Scalar,rhs,size,rhsIncr==1 ? const_cast<Scalar*>(_rhs) : 0);  
   if (rhsIncr!=1)
   {
-    Scalar* r = ei_aligned_stack_new(Scalar, size);
     const Scalar* it = _rhs;
     for (Index i=0; i<size; ++i, it+=rhsIncr)
-      r[i] = *it;
-    rhs = r;
+      rhs[i] = *it;
   }
 
   Index bound = std::max(Index(0),size-8) & 0xfffffffe;
@@ -160,9 +158,6 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
     }
     res[j] += alpha * t2;
   }
-
-  if(rhsIncr!=1)
-    ei_aligned_stack_delete(Scalar, const_cast<Scalar*>(rhs), size);
 }
 
 } // end namespace internal 
@@ -211,40 +206,28 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
     
     internal::gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,!EvalToDest> static_dest;
     internal::gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!UseRhs> static_rhs;
+
+    ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
+                                                  EvalToDest ? dest.data() : static_dest.data());
+                                                  
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,rhs.size(),
+        UseRhs ? const_cast<RhsScalar*>(rhs.data()) : static_rhs.data());
     
-    bool freeDestPtr = false;
-    ResScalar* actualDestPtr;
-    if(EvalToDest)
-      actualDestPtr = dest.data();
-    else
+    if(!EvalToDest)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = dest.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualDestPtr=static_dest.data())==0)
-      {
-        freeDestPtr = true;
-        actualDestPtr = ei_aligned_stack_new(ResScalar,dest.size());
-      }
       MappedDest(actualDestPtr, dest.size()) = dest;
     }
       
-    bool freeRhsPtr = false;
-    RhsScalar* actualRhsPtr;
-    if(UseRhs)
-      actualRhsPtr = const_cast<RhsScalar*>(rhs.data());
-    else
+    if(!UseRhs)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = rhs.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualRhsPtr=static_rhs.data())==0)
-      {
-        freeRhsPtr = true;
-        actualRhsPtr = ei_aligned_stack_new(RhsScalar,rhs.size());
-      }
       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, rhs.size()) = rhs;
     }
       
@@ -259,11 +242,7 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
       );
     
     if(!EvalToDest)
-    {
       dest = MappedDest(actualDestPtr, dest.size());
-      if(freeDestPtr) ei_aligned_stack_delete(ResScalar, actualDestPtr, dest.size());
-    }
-    if(freeRhsPtr) ei_aligned_stack_delete(RhsScalar, actualRhsPtr, rhs.size());
   }
 };
 

Eigen/src/Core/products/SelfadjointProduct.h

@@ -81,27 +81,17 @@ struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
       UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
     };
     internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
-    
-    bool freeOtherPtr = false;
-    Scalar* actualOtherPtr;
-    if(UseOtherDirectly)
-      actualOtherPtr = const_cast<Scalar*>(actualOther.data());
-    else
-    {
-      if((actualOtherPtr=static_other.data())==0)
-      {
-        freeOtherPtr = true;
-        actualOtherPtr = ei_aligned_stack_new(Scalar,other.size());
-      }
+
+    ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
+      (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
+      
+    if(!UseOtherDirectly)
       Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
-    }
     
     selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
                               OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
                             (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
           ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
-    
-    if((!UseOtherDirectly) && freeOtherPtr) ei_aligned_stack_delete(Scalar, actualOtherPtr, other.size());
   }
 };
 

Eigen/src/Core/products/TriangularMatrixMatrix.h

@@ -96,33 +96,38 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
                                            LhsStorageOrder,ConjugateLhs,
                                            RhsStorageOrder,ConjugateRhs,ColMajor>
 {
+  
+  typedef gebp_traits<Scalar,Scalar> Traits;
+  enum {
+    SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
+    IsLower = (Mode&Lower) == Lower,
+    SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
+  };
 
   static EIGEN_DONT_INLINE void run(
-    Index rows, Index cols, Index depth,
+    Index _rows, Index _cols, Index _depth,
     const Scalar* _lhs, Index lhsStride,
     const Scalar* _rhs, Index rhsStride,
     Scalar* res,        Index resStride,
     Scalar alpha)
   {
+    // strip zeros
+    Index diagSize  = std::min(_rows,_depth);
+    Index rows      = IsLower ? _rows : diagSize;
+    Index depth     = IsLower ? diagSize : _depth;
+    Index cols      = _cols;
+    
     const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
     const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
 
-    typedef gebp_traits<Scalar,Scalar> Traits;
-    enum {
-      SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
-      IsLower = (Mode&Lower) == Lower,
-      SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
-    };
-
     Index kc = depth; // cache block size along the K direction
     Index mc = rows;  // cache block size along the M direction
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
-
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);    
     Scalar* blockB = allocatedBlockB + sizeW;
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,LhsStorageOrder> triangularBuffer;
@@ -153,10 +158,11 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
       pack_rhs(blockB, &rhs(actual_k2,0), rhsStride, actual_kc, cols);
 
       // the selected lhs's panel has to be split in three different parts:
-      //  1 - the part which is above the diagonal block => skip it
+      //  1 - the part which is zero => skip it
       //  2 - the diagonal block => special kernel
-      //  3 - the panel below the diagonal block => GEPP
-      // the block diagonal, if any
+      //  3 - the dense panel below (lower case) or above (upper case) the diagonal block => GEPP
+
+      // the block diagonal, if any:
       if(IsLower || actual_k2<rows)
       {
         // for each small vertical panels of lhs
@@ -194,7 +200,7 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
           }
         }
       }
-      // the part below the diagonal => GEPP
+      // the part below (lower case) or above (upper case) the diagonal => GEPP
       {
         Index start = IsLower ? k2 : 0;
         Index end   = IsLower ? rows : std::min(actual_k2,rows);
@@ -208,10 +214,6 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
         }
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
-//     delete[] allocatedBlockB;
   }
 };
 
@@ -223,33 +225,38 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,false,
                                            LhsStorageOrder,ConjugateLhs,
                                            RhsStorageOrder,ConjugateRhs,ColMajor>
 {
+  typedef gebp_traits<Scalar,Scalar> Traits;
+  enum {
+    SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
+    IsLower = (Mode&Lower) == Lower,
+    SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
+  };
 
   static EIGEN_DONT_INLINE void run(
-    Index rows, Index cols, Index depth,
+    Index _rows, Index _cols, Index _depth,
     const Scalar* _lhs, Index lhsStride,
     const Scalar* _rhs, Index rhsStride,
     Scalar* res,        Index resStride,
     Scalar alpha)
   {
+    // strip zeros
+    Index diagSize  = std::min(_cols,_depth);
+    Index rows      = _rows;
+    Index depth     = IsLower ? _depth : diagSize;
+    Index cols      = IsLower ? diagSize : _cols;
+    
     const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
     const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
 
-    typedef gebp_traits<Scalar,Scalar> Traits;
-    enum {
-      SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
-      IsLower = (Mode&Lower) == Lower,
-      SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
-    };
-
     Index kc = depth; // cache block size along the K direction
     Index mc = rows;  // cache block size along the M direction
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar,sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,RhsStorageOrder> triangularBuffer;
@@ -347,9 +354,6 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,false,
                     -1, -1, 0, 0, allocatedBlockB);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/TriangularMatrixVector.h

@@ -41,9 +41,6 @@ struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,C
   static EIGEN_DONT_INLINE  void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
                                      const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
   {
-    EIGEN_UNUSED_VARIABLE(resIncr);
-    eigen_assert(resIncr==1);
-    
     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
 
     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
@@ -95,9 +92,6 @@ struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,C
   static void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
                   const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
   {
-    eigen_assert(rhsIncr==1);
-    EIGEN_UNUSED_VARIABLE(rhsIncr);
-    
     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
 
     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
@@ -185,7 +179,7 @@ struct TriangularProduct<Mode,false,Lhs,true,Rhs,false>
   template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
   {
     eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
-    
+
     typedef TriangularProduct<(Mode & UnitDiag) | ((Mode & Lower) ? Upper : Lower),true,Transpose<const Rhs>,false,Transpose<const Lhs>,true> TriangularProductTranspose;
     Transpose<Dest> dstT(dst);
     internal::trmv_selector<(int(internal::traits<Rhs>::Flags)&RowMajorBit) ? ColMajor : RowMajor>::run(
@@ -235,23 +229,15 @@ template<> struct trmv_selector<ColMajor>
     
     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
 
-    ResScalar* actualDestPtr;
-    bool freeDestPtr = false;
-    if (evalToDest)
-    {
-      actualDestPtr = dest.data();
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
+                                                  evalToDest ? dest.data() : static_dest.data());
+
+    if(!evalToDest)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = dest.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualDestPtr = static_dest.data())==0)
-      {
-        freeDestPtr = true;
-        actualDestPtr = ei_aligned_stack_new(ResScalar,dest.size());
-      }
       if(!alphaIsCompatible)
       {
         MappedDest(actualDestPtr, dest.size()).setZero();
@@ -277,7 +263,6 @@ template<> struct trmv_selector<ColMajor>
         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
       else
         dest = MappedDest(actualDestPtr, dest.size());
-      if(freeDestPtr) ei_aligned_stack_delete(ResScalar, actualDestPtr, dest.size());
     }
   }
 };
@@ -310,23 +295,15 @@ template<> struct trmv_selector<RowMajor>
 
     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
 
-    RhsScalar* actualRhsPtr;
-    bool freeRhsPtr = false;
-    if (DirectlyUseRhs)
-    {
-      actualRhsPtr = const_cast<RhsScalar*>(actualRhs.data());
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
+        DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
+
+    if(!DirectlyUseRhs)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = actualRhs.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualRhsPtr = static_rhs.data())==0)
-      {
-        freeRhsPtr = true;
-        actualRhsPtr = ei_aligned_stack_new(RhsScalar, actualRhs.size());
-      }
       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
     }
     
@@ -340,8 +317,6 @@ template<> struct trmv_selector<RowMajor>
             actualRhsPtr,1,
             dest.data(),dest.innerStride(),
             actualAlpha);
-
-    if((!DirectlyUseRhs) && freeRhsPtr) ei_aligned_stack_delete(RhsScalar, actualRhsPtr, prod.rhs().size());
   }
 };
 

Eigen/src/Core/products/TriangularSolverMatrix.h

@@ -70,10 +70,10 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     conj_if<Conjugate> conj;
@@ -174,9 +174,6 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
         }
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 
@@ -209,10 +206,10 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
     Index nc = rows;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*size;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     conj_if<Conjugate> conj;
@@ -314,9 +311,6 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
                       -1, -1, 0, 0, allocatedBlockB);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/util/Constants.h

@@ -161,23 +161,72 @@ const unsigned int HereditaryBits = RowMajorBit
                                   | EvalBeforeNestingBit
                                   | EvalBeforeAssigningBit;
 
-// Possible values for the Mode parameter of triangularView()
-enum {
-  Lower=0x1, Upper=0x2, UnitDiag=0x4, ZeroDiag=0x8,
-  UnitLower=UnitDiag|Lower, UnitUpper=UnitDiag|Upper,
-  StrictlyLower=ZeroDiag|Lower, StrictlyUpper=ZeroDiag|Upper,
-  SelfAdjoint=0x10};
+/** \defgroup enums Enumerations
+  * \ingroup Core_Module
+  *
+  * Various enumerations used in %Eigen. Many of these are used as template parameters.
+  */
+
+/** \ingroup enums
+  * Enum containing possible values for the \p Mode parameter of 
+  * MatrixBase::selfadjointView() and MatrixBase::triangularView(). */
+enum {
+  /** View matrix as a lower triangular matrix. */
+  Lower=0x1,                      
+  /** View matrix as an upper triangular matrix. */
+  Upper=0x2,                      
+  /** %Matrix has ones on the diagonal; to be used in combination with #Lower or #Upper. */
+  UnitDiag=0x4, 
+  /** %Matrix has zeros on the diagonal; to be used in combination with #Lower or #Upper. */
+  ZeroDiag=0x8,
+  /** View matrix as a lower triangular matrix with ones on the diagonal. */
+  UnitLower=UnitDiag|Lower, 
+  /** View matrix as an upper triangular matrix with ones on the diagonal. */
+  UnitUpper=UnitDiag|Upper,
+  /** View matrix as a lower triangular matrix with zeros on the diagonal. */
+  StrictlyLower=ZeroDiag|Lower, 
+  /** View matrix as an upper triangular matrix with zeros on the diagonal. */
+  StrictlyUpper=ZeroDiag|Upper,
+  /** Used in BandMatrix and SelfAdjointView to indicate that the matrix is self-adjoint. */
+  SelfAdjoint=0x10
+};
+
+/** \ingroup enums
+  * Enum for indicating whether an object is aligned or not. */
+enum { 
+  /** Object is not correctly aligned for vectorization. */
+  Unaligned=0, 
+  /** Object is aligned for vectorization. */
+  Aligned=1 
+};
 
-enum { Unaligned=0, Aligned=1 };
 enum { ConditionalJumpCost = 5 };
 
+/** \ingroup enums
+ * Enum used by DenseBase::corner() in Eigen2 compatibility mode. */
 // FIXME after the corner() API change, this was not needed anymore, except by AlignedBox
 // TODO: find out what to do with that. Adapt the AlignedBox API ?
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 
-enum DirectionType { Vertical, Horizontal, BothDirections };
+/** \ingroup enums
+  * Enum containing possible values for the \p Direction parameter of
+  * Reverse, PartialReduxExpr and VectorwiseOp. */
+enum DirectionType { 
+  /** For Reverse, all columns are reversed; 
+    * for PartialReduxExpr and VectorwiseOp, act on columns. */
+  Vertical, 
+  /** For Reverse, all rows are reversed; 
+    * for PartialReduxExpr and VectorwiseOp, act on rows. */
+  Horizontal, 
+  /** For Reverse, both rows and columns are reversed; 
+    * not used for PartialReduxExpr and VectorwiseOp. */
+  BothDirections 
+};
+
 enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct };
 
+/** \internal \ingroup enums
+  * Enum to specify how to traverse the entries of a matrix. */
 enum {
   /** \internal Default traversal, no vectorization, no index-based access */
   DefaultTraversal,
@@ -196,14 +245,25 @@ enum {
   InvalidTraversal
 };
 
+/** \internal \ingroup enums
+  * Enum to specify whether to unroll loops when traversing over the entries of a matrix. */
 enum {
+  /** \internal Do not unroll loops. */
   NoUnrolling,
+  /** \internal Unroll only the inner loop, but not the outer loop. */
   InnerUnrolling,
+  /** \internal Unroll both the inner and the outer loop. If there is only one loop, 
+    * because linear traversal is used, then unroll that loop. */
   CompleteUnrolling
 };
 
+/** \ingroup enums
+  * Enum containing possible values for the \p _Options template parameter of
+  * Matrix, Array and BandMatrix. */
 enum {
+  /** Storage order is column major (see \ref TopicStorageOrders). */
   ColMajor = 0,
+  /** Storage order is row major (see \ref TopicStorageOrders). */
   RowMajor = 0x1,  // it is only a coincidence that this is equal to RowMajorBit -- don't rely on that
   /** \internal Align the matrix itself if it is vectorizable fixed-size */
   AutoAlign = 0,
@@ -211,11 +271,13 @@ enum {
   DontAlign = 0x2
 };
 
-/** \brief Enum for specifying whether to apply or solve on the left or right. 
-  */
+/** \ingroup enums
+  * Enum for specifying whether to apply or solve on the left or right. */
 enum {
-  OnTheLeft = 1,  /**< \brief Apply transformation on the left. */
-  OnTheRight = 2  /**< \brief Apply transformation on the right. */
+  /** Apply transformation on the left. */
+  OnTheLeft = 1,  
+  /** Apply transformation on the right. */
+  OnTheRight = 2  
 };
 
 /* the following could as well be written:
@@ -239,53 +301,108 @@ namespace {
   EIGEN_UNUSED Default_t Default;
 }
 
+/** \internal \ingroup enums
+  * Used in AmbiVector. */
 enum {
   IsDense         = 0,
   IsSparse
 };
 
+/** \ingroup enums
+  * Used as template parameter in DenseCoeffBase and MapBase to indicate 
+  * which accessors should be provided. */
 enum AccessorLevels {
-  ReadOnlyAccessors, WriteAccessors, DirectAccessors, DirectWriteAccessors
+  /** Read-only access via a member function. */
+  ReadOnlyAccessors, 
+  /** Read/write access via member functions. */
+  WriteAccessors, 
+  /** Direct read-only access to the coefficients. */
+  DirectAccessors, 
+  /** Direct read/write access to the coefficients. */
+  DirectWriteAccessors
 };
 
+/** \ingroup enums
+  * Enum with options to give to various decompositions. */
 enum DecompositionOptions {
-  Pivoting            = 0x01, // LDLT,
-  NoPivoting          = 0x02, // LDLT,
-  ComputeFullU        = 0x04, // SVD,
-  ComputeThinU        = 0x08, // SVD,
-  ComputeFullV        = 0x10, // SVD,
-  ComputeThinV        = 0x20, // SVD,
-  EigenvaluesOnly     = 0x40, // all eigen solvers
-  ComputeEigenvectors = 0x80, // all eigen solvers
+  /** \internal Not used (meant for LDLT?). */
+  Pivoting            = 0x01, 
+  /** \internal Not used (meant for LDLT?). */
+  NoPivoting          = 0x02, 
+  /** Used in JacobiSVD to indicate that the square matrix U is to be computed. */
+  ComputeFullU        = 0x04,
+  /** Used in JacobiSVD to indicate that the thin matrix U is to be computed. */
+  ComputeThinU        = 0x08,
+  /** Used in JacobiSVD to indicate that the square matrix V is to be computed. */
+  ComputeFullV        = 0x10,
+  /** Used in JacobiSVD to indicate that the thin matrix V is to be computed. */
+  ComputeThinV        = 0x20,
+  /** Used in SelfAdjointEigenSolver and GeneralizedSelfAdjointEigenSolver to specify
+    * that only the eigenvalues are to be computed and not the eigenvectors. */
+  EigenvaluesOnly     = 0x40,
+  /** Used in SelfAdjointEigenSolver and GeneralizedSelfAdjointEigenSolver to specify
+    * that both the eigenvalues and the eigenvectors are to be computed. */
+  ComputeEigenvectors = 0x80,
+  /** \internal */
   EigVecMask = EigenvaluesOnly | ComputeEigenvectors,
+  /** Used in GeneralizedSelfAdjointEigenSolver to indicate that it should
+    * solve the generalized eigenproblem \f$ Ax = \lambda B x \f$. */
   Ax_lBx              = 0x100,
+  /** Used in GeneralizedSelfAdjointEigenSolver to indicate that it should
+    * solve the generalized eigenproblem \f$ ABx = \lambda x \f$. */
   ABx_lx              = 0x200,
+  /** Used in GeneralizedSelfAdjointEigenSolver to indicate that it should
+    * solve the generalized eigenproblem \f$ BAx = \lambda x \f$. */
   BAx_lx              = 0x400,
+  /** \internal */
   GenEigMask = Ax_lBx | ABx_lx | BAx_lx
 };
 
+/** \ingroup enums
+  * Possible values for the \p QRPreconditioner template parameter of JacobiSVD. */
 enum QRPreconditioners {
+  /** Do not specify what is to be done if the SVD of a non-square matrix is asked for. */
   NoQRPreconditioner,
+  /** Use a QR decomposition without pivoting as the first step. */
   HouseholderQRPreconditioner,
+  /** Use a QR decomposition with column pivoting as the first step. */
   ColPivHouseholderQRPreconditioner,
+  /** Use a QR decomposition with full pivoting as the first step. */
   FullPivHouseholderQRPreconditioner
 };
 
-/** \brief Enum for reporting the status of a computation.
-  */
+#ifdef Success
+#error The preprocessor symbol 'Success' is defined, possibly by the X11 header file X.h
+#endif
+
+/** \ingroups enums
+  * Enum for reporting the status of a computation. */
 enum ComputationInfo {
-  Success = 0,        /**< \brief Computation was successful. */
-  NumericalIssue = 1, /**< \brief The provided data did not satisfy the prerequisites. */
-  NoConvergence = 2   /**< \brief Iterative procedure did not converge. */
+  /** Computation was successful. */
+  Success = 0,        
+  /** The provided data did not satisfy the prerequisites. */
+  NumericalIssue = 1, 
+  /** Iterative procedure did not converge. */
+  NoConvergence = 2
 };
 
+/** \ingroup enums
+  * Enum used to specify how a particular transformation is stored in a matrix.
+  * \sa Transform, Hyperplane::transform(). */
 enum TransformTraits {
+  /** Transformation is an isometry. */
   Isometry      = 0x1,
+  /** Transformation is an affine transformation stored as a (Dim+1)^2 matrix whose last row is 
+    * assumed to be [0 ... 0 1]. */
   Affine        = 0x2,
+  /** Transformation is an affine transformation stored as a (Dim) x (Dim+1) matrix. */
   AffineCompact = 0x10 | Affine,
+  /** Transformation is a general projective transformation stored as a (Dim+1)^2 matrix. */
   Projective    = 0x20
 };
 
+/** \internal \ingroup enums
+  * Enum used to choose between implementation depending on the computer architecture. */
 namespace Architecture
 {
   enum Type {
@@ -302,8 +419,12 @@ namespace Architecture
   };
 }
 
+/** \internal \ingroup enums
+  * Enum used as template parameter in GeneralProduct. */
 enum { CoeffBasedProductMode, LazyCoeffBasedProductMode, OuterProduct, InnerProduct, GemvProduct, GemmProduct };
 
+/** \internal \ingroup enums
+  * Enum used in experimental parallel implementation. */
 enum Action {GetAction, SetAction};
 
 /** The type used to identify a dense storage. */

Eigen/src/Core/util/DisableMSVCWarnings.h

@@ -1,17 +0,0 @@
-
-#ifdef _MSC_VER
-  // 4100 - unreferenced formal parameter (occurred e.g. in aligned_allocator::destroy(pointer p))
-  // 4101 - unreferenced local variable
-  // 4127 - conditional expression is constant
-  // 4181 - qualifier applied to reference type ignored
-  // 4211 - nonstandard extension used : redefined extern to static
-  // 4244 - 'argument' : conversion from 'type1' to 'type2', possible loss of data
-  // 4273 - QtAlignedMalloc, inconsistent DLL linkage
-  // 4324 - structure was padded due to declspec(align())
-  // 4512 - assignment operator could not be generated
-  // 4522 - 'class' : multiple assignment operators specified
-  // 4700 - uninitialized local variable 'xyz' used
-  // 4717 - 'function' : recursive on all control paths, function will cause runtime stack overflow
-  #pragma warning( push )
-  #pragma warning( disable : 4100 4101 4127 4181 4211 4244 4273 4324 4512 4522 4700 4717 )
-#endif

Eigen/src/Core/util/DisableStupidWarnings.h

@@ -0,0 +1,42 @@
+#ifndef EIGEN_WARNINGS_DISABLED
+#define EIGEN_WARNINGS_DISABLED
+
+#ifdef _MSC_VER
+  // 4100 - unreferenced formal parameter (occurred e.g. in aligned_allocator::destroy(pointer p))
+  // 4101 - unreferenced local variable
+  // 4127 - conditional expression is constant
+  // 4181 - qualifier applied to reference type ignored
+  // 4211 - nonstandard extension used : redefined extern to static
+  // 4244 - 'argument' : conversion from 'type1' to 'type2', possible loss of data
+  // 4273 - QtAlignedMalloc, inconsistent DLL linkage
+  // 4324 - structure was padded due to declspec(align())
+  // 4512 - assignment operator could not be generated
+  // 4522 - 'class' : multiple assignment operators specified
+  // 4700 - uninitialized local variable 'xyz' used
+  // 4717 - 'function' : recursive on all control paths, function will cause runtime stack overflow
+  #ifndef EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
+    #pragma warning( push )
+  #endif
+  #pragma warning( disable : 4100 4101 4127 4181 4211 4244 4273 4324 4512 4522 4700 4717 )
+#elif defined __INTEL_COMPILER
+  // 2196 - routine is both "inline" and "noinline" ("noinline" assumed)
+  //        ICC 12 generates this warning even without any inline keyword, when defining class methods 'inline' i.e. inside of class body
+  // 2536 - type qualifiers are meaningless here
+  //        ICC 12 generates this warning when a function return type is const qualified, even if that type is a template-parameter-dependent
+  //        typedef that may be a reference type.
+  // 279  - controlling expression is constant
+  //        ICC 12 generates this warning on assert(constant_expression_depending_on_template_params) and frankly this is a legitimate use case.
+  #ifndef EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
+    #pragma warning push
+  #endif
+  #pragma warning disable 2196 2536 279
+#elif defined __clang__
+  // -Wconstant-logical-operand - warning: use of logical && with constant operand; switch to bitwise & or remove constant
+  //     this is really a stupid warning as it warns on compile-time expressions involving enums
+  #ifndef EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
+    #pragma clang diagnostic push
+  #endif
+  #pragma clang diagnostic ignored "-Wconstant-logical-operand"
+#endif
+
+#endif // not EIGEN_WARNINGS_DISABLED

Eigen/src/Core/util/EnableMSVCWarnings.h

@@ -1,4 +0,0 @@
-
-#ifdef _MSC_VER
-  #pragma warning( pop )
-#endif

Eigen/src/Core/util/ForwardDeclarations.h

@@ -179,6 +179,9 @@ template<typename Scalar> struct scalar_exp_op;
 template<typename Scalar> struct scalar_log_op;
 template<typename Scalar> struct scalar_cos_op;
 template<typename Scalar> struct scalar_sin_op;
+template<typename Scalar> struct scalar_acos_op;
+template<typename Scalar> struct scalar_asin_op;
+template<typename Scalar> struct scalar_tan_op;
 template<typename Scalar> struct scalar_pow_op;
 template<typename Scalar> struct scalar_inverse_op;
 template<typename Scalar> struct scalar_square_op;

Eigen/src/Core/util/Macros.h

@@ -26,18 +26,31 @@
 #ifndef EIGEN_MACROS_H
 #define EIGEN_MACROS_H
 
-#define EIGEN_WORLD_VERSION 2
-#define EIGEN_MAJOR_VERSION 93
-#define EIGEN_MINOR_VERSION 0
+#define EIGEN_WORLD_VERSION 3
+#define EIGEN_MAJOR_VERSION 0
+#define EIGEN_MINOR_VERSION 1
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
                                                                  EIGEN_MINOR_VERSION>=z))))
 #ifdef __GNUC__
-  #define EIGEN_GNUC_AT_LEAST(x,y) ((__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x)
+  #define EIGEN_GNUC_AT_LEAST(x,y) ((__GNUC__==x && __GNUC_MINOR__>=y) || __GNUC__>x)
 #else
   #define EIGEN_GNUC_AT_LEAST(x,y) 0
 #endif
+ 
+#ifdef __GNUC__
+  #define EIGEN_GNUC_AT_MOST(x,y) ((__GNUC__==x && __GNUC_MINOR__<=y) || __GNUC__<x)
+#else
+  #define EIGEN_GNUC_AT_MOST(x,y) 0
+#endif
+
+#if EIGEN_GNUC_AT_MOST(4,3)
+  // see bug 89
+  #define EIGEN_SAFE_TO_USE_STANDARD_ASSERT_MACRO 0
+#else
+  #define EIGEN_SAFE_TO_USE_STANDARD_ASSERT_MACRO 1
+#endif
 
 #if defined(__GNUC__) && (__GNUC__ <= 3)
 #define EIGEN_GCC3_OR_OLDER 1
@@ -109,31 +122,13 @@
 
 #define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl;
 
-#ifdef NDEBUG
-# ifndef EIGEN_NO_DEBUG
-#  define EIGEN_NO_DEBUG
-# endif
-#endif
+// concatenate two tokens
+#define EIGEN_CAT2(a,b) a ## b
+#define EIGEN_CAT(a,b) EIGEN_CAT2(a,b)
 
-#ifndef eigen_assert
-#ifdef EIGEN_NO_DEBUG
-#define eigen_assert(x)
-#else
-#define eigen_assert(x) assert(x)
-#endif
-#endif
-
-#ifdef EIGEN_INTERNAL_DEBUGGING
-#define eigen_internal_assert(x) eigen_assert(x)
-#else
-#define eigen_internal_assert(x)
-#endif
-
-#ifdef EIGEN_NO_DEBUG
-#define EIGEN_ONLY_USED_FOR_DEBUG(x) (void)x
-#else
-#define EIGEN_ONLY_USED_FOR_DEBUG(x)
-#endif
+// convert a token to a string
+#define EIGEN_MAKESTRING2(a) #a
+#define EIGEN_MAKESTRING(a) EIGEN_MAKESTRING2(a)
 
 // EIGEN_ALWAYS_INLINE_ATTRIB should be use in the declaration of function
 // which should be inlined even in debug mode.
@@ -147,14 +142,14 @@
 #define EIGEN_ALWAYS_INLINE_ATTRIB
 #endif
 
-#if EIGEN_GNUC_AT_LEAST(4,1)
+#if EIGEN_GNUC_AT_LEAST(4,1) && !defined(__clang__) && !defined(__INTEL_COMPILER)
 #define EIGEN_FLATTEN_ATTRIB __attribute__((flatten))
 #else
 #define EIGEN_FLATTEN_ATTRIB
 #endif
 
 // EIGEN_FORCE_INLINE means "inline as much as possible"
-#if (defined _MSC_VER) || (defined __intel_compiler)
+#if (defined _MSC_VER) || (defined __INTEL_COMPILER)
 #define EIGEN_STRONG_INLINE __forceinline
 #else
 #define EIGEN_STRONG_INLINE inline
@@ -175,6 +170,67 @@
 #define EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
 #define EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS inline
 
+#ifdef NDEBUG
+# ifndef EIGEN_NO_DEBUG
+#  define EIGEN_NO_DEBUG
+# endif
+#endif
+
+// eigen_plain_assert is where we implement the workaround for the assert() bug in GCC <= 4.3, see bug 89
+#ifdef EIGEN_NO_DEBUG
+  #define eigen_plain_assert(x)
+#else
+  #if EIGEN_SAFE_TO_USE_STANDARD_ASSERT_MACRO
+    namespace Eigen {
+    namespace internal {
+    inline bool copy_bool(bool b) { return b; }
+    }
+    }
+    #define eigen_plain_assert(x) assert(x)
+  #else
+    // work around bug 89
+    #include <cstdlib>   // for abort
+    #include <iostream>  // for std::cerr
+
+    namespace Eigen {
+    namespace internal {
+    // trivial function copying a bool. Must be EIGEN_DONT_INLINE, so we implement it after including Eigen headers.
+    // see bug 89.
+    namespace {
+    EIGEN_DONT_INLINE bool copy_bool(bool b) { return b; }
+    }
+    inline void assert_fail(const char *condition, const char *function, const char *file, int line)
+    {
+      std::cerr << "assertion failed: " << condition << " in function " << function << " at " << file << ":" << line << std::endl;
+      abort();
+    }
+    }
+    }
+    #define eigen_plain_assert(x) \
+      do { \
+        if(!Eigen::internal::copy_bool(x)) \
+          Eigen::internal::assert_fail(EIGEN_MAKESTRING(x), __PRETTY_FUNCTION__, __FILE__, __LINE__); \
+      } while(false)
+  #endif
+#endif
+
+// eigen_assert can be overridden
+#ifndef eigen_assert
+#define eigen_assert(x) eigen_plain_assert(x)
+#endif
+
+#ifdef EIGEN_INTERNAL_DEBUGGING
+#define eigen_internal_assert(x) eigen_assert(x)
+#else
+#define eigen_internal_assert(x)
+#endif
+
+#ifdef EIGEN_NO_DEBUG
+#define EIGEN_ONLY_USED_FOR_DEBUG(x) (void)x
+#else
+#define EIGEN_ONLY_USED_FOR_DEBUG(x)
+#endif
+
 #if (defined __GNUC__)
 #define EIGEN_DEPRECATED __attribute__((deprecated))
 #elif (defined _MSC_VER)
@@ -205,9 +261,7 @@
  * If we made alignment depend on whether or not EIGEN_VECTORIZE is defined, it would be impossible to link
  * vectorized and non-vectorized code.
  */
-#if !EIGEN_ALIGN_STATICALLY
-  #define EIGEN_ALIGN_TO_BOUNDARY(n)
-#elif (defined __GNUC__) || (defined __PGI) || (defined __IBMCPP__)
+#if (defined __GNUC__) || (defined __PGI) || (defined __IBMCPP__)
   #define EIGEN_ALIGN_TO_BOUNDARY(n) __attribute__((aligned(n)))
 #elif (defined _MSC_VER)
   #define EIGEN_ALIGN_TO_BOUNDARY(n) __declspec(align(n))
@@ -220,6 +274,14 @@
 
 #define EIGEN_ALIGN16 EIGEN_ALIGN_TO_BOUNDARY(16)
 
+#if EIGEN_ALIGN_STATICALLY
+#define EIGEN_USER_ALIGN_TO_BOUNDARY(n) EIGEN_ALIGN_TO_BOUNDARY(n)
+#define EIGEN_USER_ALIGN16 EIGEN_ALIGN16
+#else
+#define EIGEN_USER_ALIGN_TO_BOUNDARY(n)
+#define EIGEN_USER_ALIGN16
+#endif
+
 #ifdef EIGEN_DONT_USE_RESTRICT_KEYWORD
   #define EIGEN_RESTRICT
 #endif
@@ -244,14 +306,6 @@
 // just an empty macro !
 #define EIGEN_EMPTY
 
-// concatenate two tokens
-#define EIGEN_CAT2(a,b) a ## b
-#define EIGEN_CAT(a,b) EIGEN_CAT2(a,b)
-
-// convert a token to a string
-#define EIGEN_MAKESTRING2(a) #a
-#define EIGEN_MAKESTRING(a) EIGEN_MAKESTRING2(a)
-
 #if defined(_MSC_VER) && (!defined(__INTEL_COMPILER))
 #define EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Derived) \
   using Base::operator =;

Eigen/src/Core/util/Memory.h

@@ -167,14 +167,36 @@ inline void* generic_aligned_realloc(void* ptr, size_t size, size_t old_size)
 *** Implementation of portable aligned versions of malloc/free/realloc     ***
 *****************************************************************************/
 
+#ifdef EIGEN_NO_MALLOC
+inline void check_that_malloc_is_allowed()
+{
+  eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
+}
+#elif defined EIGEN_RUNTIME_NO_MALLOC
+inline bool is_malloc_allowed_impl(bool update, bool new_value = false)
+{
+  static bool value = true;
+  if (update == 1)
+    value = new_value;
+  return value;
+}
+inline bool is_malloc_allowed() { return is_malloc_allowed_impl(false); }
+inline bool set_is_malloc_allowed(bool new_value) { return is_malloc_allowed_impl(true, new_value); }
+inline void check_that_malloc_is_allowed()
+{
+  eigen_assert(is_malloc_allowed() && "heap allocation is forbidden (EIGEN_RUNTIME_NO_MALLOC is defined and g_is_malloc_allowed is false)");
+}
+#else 
+inline void check_that_malloc_is_allowed()
+{}
+#endif
+
 /** \internal Allocates \a size bytes. The returned pointer is guaranteed to have 16 bytes alignment.
   * On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown.
   */
 inline void* aligned_malloc(size_t size)
 {
-  #ifdef EIGEN_NO_MALLOC
-    eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
-  #endif
+  check_that_malloc_is_allowed();
 
   void *result;
   #if !EIGEN_ALIGN
@@ -268,9 +290,7 @@ template<bool Align> inline void* conditional_aligned_malloc(size_t size)
 
 template<> inline void* conditional_aligned_malloc<false>(size_t size)
 {
-  #ifdef EIGEN_NO_MALLOC
-    eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
-  #endif
+  check_that_malloc_is_allowed();
 
   void *result = std::malloc(size);
   #ifdef EIGEN_EXCEPTIONS
@@ -448,36 +468,87 @@ inline static Index first_aligned(const Scalar* array, Index size)
 *** Implementation of runtime stack allocation (falling back to malloc)    ***
 *****************************************************************************/
 
-/** \internal
-  * Allocates an aligned buffer of SIZE bytes on the stack if SIZE is smaller than
-  * EIGEN_STACK_ALLOCATION_LIMIT, and if stack allocation is supported by the platform
-  * (currently, this is Linux only). Otherwise the memory is allocated on the heap.
-  * Data allocated with ei_aligned_stack_alloc \b must be freed by calling
-  * ei_aligned_stack_free(PTR,SIZE).
-  * \code
-  * float * data = ei_aligned_stack_alloc(float,array.size());
-  * // ...
-  * ei_aligned_stack_free(data,float,array.size());
-  * \endcode
-  */
-#if (defined __linux__)
-  #define ei_aligned_stack_alloc(SIZE) (SIZE<=EIGEN_STACK_ALLOCATION_LIMIT) \
-                                    ? alloca(SIZE) \
-                                    : Eigen::internal::aligned_malloc(SIZE)
-  #define ei_aligned_stack_free(PTR,SIZE) if(SIZE>EIGEN_STACK_ALLOCATION_LIMIT) Eigen::internal::aligned_free(PTR)
-#elif defined(_MSC_VER)
-  #define ei_aligned_stack_alloc(SIZE) (SIZE<=EIGEN_STACK_ALLOCATION_LIMIT) \
-                                    ? _alloca(SIZE) \
-                                    : Eigen::internal::aligned_malloc(SIZE)
-  #define ei_aligned_stack_free(PTR,SIZE) if(SIZE>EIGEN_STACK_ALLOCATION_LIMIT) Eigen::internal::aligned_free(PTR)
-#else
-  #define ei_aligned_stack_alloc(SIZE) Eigen::internal::aligned_malloc(SIZE)
-  #define ei_aligned_stack_free(PTR,SIZE) Eigen::internal::aligned_free(PTR)
+// you can overwrite Eigen's default behavior regarding alloca by defining EIGEN_ALLOCA
+// to the appropriate stack allocation function
+#ifndef EIGEN_ALLOCA
+  #if (defined __linux__)
+    #define EIGEN_ALLOCA alloca
+  #elif defined(_MSC_VER)
+    #define EIGEN_ALLOCA _alloca
+  #endif
 #endif
 
-#define ei_aligned_stack_new(TYPE,SIZE) Eigen::internal::construct_elements_of_array(reinterpret_cast<TYPE*>(ei_aligned_stack_alloc(sizeof(TYPE)*SIZE)), SIZE)
-#define ei_aligned_stack_delete(TYPE,PTR,SIZE) do {Eigen::internal::destruct_elements_of_array<TYPE>(PTR, SIZE); \
-                                                   ei_aligned_stack_free(PTR,sizeof(TYPE)*SIZE);} while(0)
+namespace internal {
+
+// This helper class construct the allocated memory, and takes care of destructing and freeing the handled data
+// at destruction time. In practice this helper class is mainly useful to avoid memory leak in case of exceptions.
+template<typename T> class aligned_stack_memory_handler
+{
+  public:
+    /* Creates a stack_memory_handler responsible for the buffer \a ptr of size \a size.
+     * Note that \a ptr can be 0 regardless of the other parameters.
+     * This constructor takes care of constructing/initializing the elements of the buffer if required by the scalar type T (see NumTraits<T>::RequireInitialization).
+     * In this case, the buffer elements will also be destructed when this handler will be destructed.
+     * Finally, if \a dealloc is true, then the pointer \a ptr is freed.
+     **/
+    aligned_stack_memory_handler(T* ptr, size_t size, bool dealloc)
+      : m_ptr(ptr), m_size(size), m_deallocate(dealloc)
+    {
+      if(NumTraits<T>::RequireInitialization)
+        Eigen::internal::construct_elements_of_array(m_ptr, size);
+    }
+    ~aligned_stack_memory_handler()
+    {
+      if(NumTraits<T>::RequireInitialization)
+        Eigen::internal::destruct_elements_of_array<T>(m_ptr, m_size);
+      if(m_deallocate)
+        Eigen::internal::aligned_free(m_ptr);
+    }
+  protected:
+    T* m_ptr;
+    size_t m_size;
+    bool m_deallocate;
+};
+
+}
+
+/** \internal
+  * Declares, allocates and construct an aligned buffer named NAME of SIZE elements of type TYPE on the stack
+  * if SIZE is smaller than EIGEN_STACK_ALLOCATION_LIMIT, and if stack allocation is supported by the platform
+  * (currently, this is Linux and Visual Studio only). Otherwise the memory is allocated on the heap.
+  * The allocated buffer is automatically deleted when exiting the scope of this declaration.
+  * If BUFFER is non nul, then the declared variable is simply an alias for BUFFER, and no allocation/deletion occurs.
+  * Here is an example:
+  * \code
+  * {
+  *   ei_declare_aligned_stack_constructed_variable(float,data,size,0);
+  *   // use data[0] to data[size-1]
+  * }
+  * \endcode
+  * The underlying stack allocation function can controlled with the EIGEN_ALLOCA preprocessor token.
+  */
+#ifdef EIGEN_ALLOCA
+
+  #ifdef __arm__
+    #define EIGEN_ALIGNED_ALLOCA(SIZE) reinterpret_cast<void*>((reinterpret_cast<size_t>(EIGEN_ALLOCA(SIZE+16)) & ~(size_t(15))) + 16)
+  #else
+    #define EIGEN_ALIGNED_ALLOCA EIGEN_ALLOCA
+  #endif
+
+  #define ei_declare_aligned_stack_constructed_variable(TYPE,NAME,SIZE,BUFFER) \
+    TYPE* NAME = (BUFFER)!=0 ? (BUFFER) \
+               : reinterpret_cast<TYPE*>( \
+                      (sizeof(TYPE)*SIZE<=EIGEN_STACK_ALLOCATION_LIMIT) ? EIGEN_ALIGNED_ALLOCA(sizeof(TYPE)*SIZE) \
+                    : Eigen::internal::aligned_malloc(sizeof(TYPE)*SIZE) );  \
+    Eigen::internal::aligned_stack_memory_handler<TYPE> EIGEN_CAT(NAME,_stack_memory_destructor)((BUFFER)==0 ? NAME : 0,SIZE,sizeof(TYPE)*SIZE>EIGEN_STACK_ALLOCATION_LIMIT)
+
+#else
+
+  #define ei_declare_aligned_stack_constructed_variable(TYPE,NAME,SIZE,BUFFER) \
+    TYPE* NAME = (BUFFER)!=0 ? BUFFER : reinterpret_cast<TYPE*>(Eigen::internal::aligned_malloc(sizeof(TYPE)*SIZE));    \
+    Eigen::internal::aligned_stack_memory_handler<TYPE> EIGEN_CAT(NAME,_stack_memory_destructor)((BUFFER)==0 ? NAME : 0,SIZE,true)
+    
+#endif
 
 
 /*****************************************************************************

Eigen/src/Core/util/ReenableStupidWarnings.h

@@ -0,0 +1,14 @@
+#ifdef EIGEN_WARNINGS_DISABLED
+#undef EIGEN_WARNINGS_DISABLED
+
+#ifndef EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
+  #ifdef _MSC_VER
+    #pragma warning( pop )
+  #elif defined __INTEL_COMPILER
+    #pragma warning pop
+  #elif defined __clang__
+    #pragma clang diagnostic pop
+  #endif
+#endif
+
+#endif // EIGEN_WARNINGS_DISABLED

Eigen/src/Eigen2Support/CMakeLists.txt

@@ -4,3 +4,5 @@ INSTALL(FILES
   ${Eigen_Eigen2Support_SRCS}
   DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support COMPONENT Devel
   )
+
+ADD_SUBDIRECTORY(Geometry)

Eigen/src/Eigen2Support/Cwise.h

@@ -55,6 +55,9 @@
   * Example: \include MatrixBase_cwise_const.cpp
   * Output: \verbinclude MatrixBase_cwise_const.out
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_CWISE_PLUGIN.
+  *
   * \sa MatrixBase::cwise() const, MatrixBase::cwise()
   */
 template<typename ExpressionType> class Cwise

Eigen/src/Eigen2Support/Geometry/CMakeLists.txt

@@ -1,6 +1,6 @@
-FILE(GLOB Eigen_Geometry_SRCS "*.h")
+FILE(GLOB Eigen_Eigen2Support_Geometry_SRCS "*.h")
 
 INSTALL(FILES
-  ${Eigen_Geometry_SRCS}
+  ${Eigen_Eigen2Support_Geometry_SRCS}
   DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support/Geometry
   )

Eigen/src/Eigenvalues/ComplexEigenSolver.h

@@ -186,7 +186,8 @@ template<typename _MatrixType> class ComplexEigenSolver
       * This function returns a column vector containing the
       * eigenvalues. Eigenvalues are repeated according to their
       * algebraic multiplicity, so there are as many eigenvalues as
-      * rows in the matrix.
+      * rows in the matrix. The eigenvalues are not sorted in any particular
+      * order.
       *
       * Example: \include ComplexEigenSolver_eigenvalues.cpp
       * Output: \verbinclude ComplexEigenSolver_eigenvalues.out

Eigen/src/Eigenvalues/EigenSolver.h

@@ -228,6 +228,7 @@ template<typename _MatrixType> class EigenSolver
       * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2
       * blocks of the form
       * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$.
+      * These blocks are not sorted in any particular order.
       * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by
       * pseudoEigenvectors() satisfy \f$ AV = VD \f$.
       *
@@ -244,7 +245,8 @@ template<typename _MatrixType> class EigenSolver
       * compute(const MatrixType&, bool) has been called before.
       *
       * The eigenvalues are repeated according to their algebraic multiplicity,
-      * so there are as many eigenvalues as rows in the matrix.
+      * so there are as many eigenvalues as rows in the matrix. The eigenvalues 
+      * are not sorted in any particular order.
       *
       * Example: \include EigenSolver_eigenvalues.cpp
       * Output: \verbinclude EigenSolver_eigenvalues.out
@@ -341,6 +343,7 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
     {
       // we have a real eigen value
       matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
+      matV.col(j).normalize();
     }
     else
     {
@@ -447,7 +450,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
     Scalar q = m_eivalues.coeff(n).imag();
 
     // Scalar vector
-    if (q == 0)
+    if (q == Scalar(0))
     {
       Scalar lastr=0, lastw=0;
       Index l = n;
@@ -488,12 +491,12 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
 
           // Overflow control
           Scalar t = internal::abs(m_matT.coeff(i,n));
-          if ((eps * t) * t > 1)
+          if ((eps * t) * t > Scalar(1))
             m_matT.col(n).tail(size-i) /= t;
         }
       }
     }
-    else if (q < 0) // Complex vector
+    else if (q < Scalar(0) && n > 0) // Complex vector
     {
       Scalar lastra=0, lastsa=0, lastw=0;
       Index l = n-1;
@@ -527,7 +530,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
         else
         {
           l = i;
-          if (m_eivalues.coeff(i).imag() == 0)
+          if (m_eivalues.coeff(i).imag() == RealScalar(0))
           {
             std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
             m_matT.coeffRef(i,n-1) = internal::real(cc);
@@ -560,13 +563,18 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
           }
 
           // Overflow control
-          Scalar t = std::max(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
-          if ((eps * t) * t > 1)
+          using std::max;
+          Scalar t = max(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
+          if ((eps * t) * t > Scalar(1))
             m_matT.block(i, n-1, size-i, 2) /= t;
 
         }
       }
     }
+    else
+    {
+      eigen_assert("Internal bug in EigenSolver"); // this should not happen
+    }
   }
 
   // Back transformation to get eigenvectors of original matrix

Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h

@@ -70,13 +70,9 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
     /** \brief Default constructor for fixed-size matrices.
       *
       * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via compute(const MatrixType&, bool) or
-      * compute(const MatrixType&, const MatrixType&, bool). This constructor
+      * perform decompositions via compute(). This constructor
       * can only be used if \p _MatrixType is a fixed-size matrix; use
-      * SelfAdjointEigenSolver(Index) for dynamic-size matrices.
-      *
-      * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
-      * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
+      * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
       */
     GeneralizedSelfAdjointEigenSolver() : Base() {}
 
@@ -86,12 +82,11 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
       * eigenvalues and eigenvectors will be computed.
       *
       * This constructor is useful for dynamic-size matrices, when the user
-      * intends to perform decompositions via compute(const MatrixType&, bool)
-      * or compute(const MatrixType&, const MatrixType&, bool). The \p size
+      * intends to perform decompositions via compute(). The \p size
       * parameter is only used as a hint. It is not an error to give a wrong
       * \p size, but it may impair performance.
       *
-      * \sa compute(const MatrixType&, bool) for an example
+      * \sa compute() for an example
       */
     GeneralizedSelfAdjointEigenSolver(Index size)
         : Base(size)
@@ -103,8 +98,8 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
       *                   Only the lower triangular part of the matrix is referenced.
       * \param[in]  matB  Positive-definite matrix in matrix pencil.
       *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}.
-      *                     Default is ComputeEigenvectors|Ax_lBx.
+      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
+      *                     Default is #ComputeEigenvectors|#Ax_lBx.
       *
       * This constructor calls compute(const MatrixType&, const MatrixType&, int)
       * to compute the eigenvalues and (if requested) the eigenvectors of the
@@ -136,8 +131,8 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
       *                   Only the lower triangular part of the matrix is referenced.
       * \param[in]  matB  Positive-definite matrix in matrix pencil.
       *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}.
-      *                     Default is ComputeEigenvectors|Ax_lBx.
+      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
+      *                     Default is #ComputeEigenvectors|#Ax_lBx.
       *
       * \returns    Reference to \c *this
       *

Eigen/src/Eigenvalues/RealSchur.h

@@ -324,11 +324,11 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scal
   m_matT.coeffRef(iu,iu) += exshift;
   m_matT.coeffRef(iu-1,iu-1) += exshift;
 
-  if (q >= 0) // Two real eigenvalues
+  if (q >= Scalar(0)) // Two real eigenvalues
   {
     Scalar z = internal::sqrt(internal::abs(q));
     JacobiRotation<Scalar> rot;
-    if (p >= 0)
+    if (p >= Scalar(0))
       rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
     else
       rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
@@ -369,7 +369,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
   {
     Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
     s = s * s + shiftInfo.coeff(2);
-    if (s > 0)
+    if (s > Scalar(0))
     {
       s = internal::sqrt(s);
       if (shiftInfo.coeff(1) < shiftInfo.coeff(0))

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -62,12 +62,12 @@ class GeneralizedSelfAdjointEigenSolver;
   *
   * Call the function compute() to compute the eigenvalues and eigenvectors of
   * a given matrix. Alternatively, you can use the
-  * SelfAdjointEigenSolver(const MatrixType&, bool) constructor which computes
+  * SelfAdjointEigenSolver(const MatrixType&, int) constructor which computes
   * the eigenvalues and eigenvectors at construction time. Once the eigenvalue
   * and eigenvectors are computed, they can be retrieved with the eigenvalues()
   * and eigenvectors() functions.
   *
-  * The documentation for SelfAdjointEigenSolver(const MatrixType&, bool)
+  * The documentation for SelfAdjointEigenSolver(const MatrixType&, int)
   * contains an example of the typical use of this class.
   *
   * To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and
@@ -110,8 +110,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     /** \brief Default constructor for fixed-size matrices.
       *
       * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via compute(const MatrixType&, bool) or
-      * compute(const MatrixType&, const MatrixType&, bool). This constructor
+      * perform decompositions via compute(). This constructor
       * can only be used if \p _MatrixType is a fixed-size matrix; use
       * SelfAdjointEigenSolver(Index) for dynamic-size matrices.
       *
@@ -131,12 +130,11 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       * eigenvalues and eigenvectors will be computed.
       *
       * This constructor is useful for dynamic-size matrices, when the user
-      * intends to perform decompositions via compute(const MatrixType&, bool)
-      * or compute(const MatrixType&, const MatrixType&, bool). The \p size
+      * intends to perform decompositions via compute(). The \p size
       * parameter is only used as a hint. It is not an error to give a wrong
       * \p size, but it may impair performance.
       *
-      * \sa compute(const MatrixType&, bool) for an example
+      * \sa compute() for an example
       */
     SelfAdjointEigenSolver(Index size)
         : m_eivec(size, size),
@@ -149,17 +147,16 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       *
       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
       *    be computed. Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
+      * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
       *
-      * This constructor calls compute(const MatrixType&, bool) to compute the
+      * This constructor calls compute(const MatrixType&, int) to compute the
       * eigenvalues of the matrix \p matrix. The eigenvectors are computed if
-      * \p options equals ComputeEigenvectors.
+      * \p options equals #ComputeEigenvectors.
       *
       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp
       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out
       *
-      * \sa compute(const MatrixType&, bool),
-      *     SelfAdjointEigenSolver(const MatrixType&, const MatrixType&, bool)
+      * \sa compute(const MatrixType&, int)
       */
     SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors)
       : m_eivec(matrix.rows(), matrix.cols()),
@@ -174,11 +171,11 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       *
       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
       *    be computed. Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
+      * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
       * \returns    Reference to \c *this
       *
       * This function computes the eigenvalues of \p matrix.  The eigenvalues()
-      * function can be used to retrieve them.  If \p options equals ComputeEigenvectors,
+      * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
       * then the eigenvectors are also computed and can be retrieved by
       * calling eigenvectors().
       *
@@ -198,11 +195,11 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       * Example: \include SelfAdjointEigenSolver_compute_MatrixType.cpp
       * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType.out
       *
-      * \sa SelfAdjointEigenSolver(const MatrixType&, bool)
+      * \sa SelfAdjointEigenSolver(const MatrixType&, int)
       */
     SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
 
-    /** \brief Returns the eigenvectors of given matrix (pencil).
+    /** \brief Returns the eigenvectors of given matrix.
       *
       * \returns  A const reference to the matrix whose columns are the eigenvectors.
       *
@@ -227,14 +224,15 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       return m_eivec;
     }
 
-    /** \brief Returns the eigenvalues of given matrix (pencil).
+    /** \brief Returns the eigenvalues of given matrix.
       *
       * \returns A const reference to the column vector containing the eigenvalues.
       *
       * \pre The eigenvalues have been computed before.
       *
       * The eigenvalues are repeated according to their algebraic multiplicity,
-      * so there are as many eigenvalues as rows in the matrix.
+      * so there are as many eigenvalues as rows in the matrix. The eigenvalues
+      * are sorted in increasing order.
       *
       * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp
       * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out
@@ -402,10 +400,11 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
 
   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
   RealScalar scale = matrix.cwiseAbs().maxCoeff();
+  if(scale==Scalar(0)) scale = 1;
   mat = matrix / scale;
   m_subdiag.resize(n-1);
   internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
-
+  
   Index end = n-1;
   Index start = 0;
   Index iter = 0; // number of iterations we are working on one element
@@ -458,7 +457,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
       }
     }
   }
-
+  
   // scale back the eigen values
   m_eivalues *= scale;
 
@@ -471,12 +470,17 @@ namespace internal {
 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
 {
+  // NOTE this version avoids over & underflow, however since the matrix is prescaled, overflow cannot occur,
+  // and underflows should be meaningless anyway. So I don't any reason to enable this version, but I keep
+  // it here for reference:
+//   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
+//   RealScalar e = subdiag[end-1];
+//   RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
   RealScalar e2 = abs2(subdiag[end-1]);
   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
   RealScalar x = diag[start] - mu;
   RealScalar z = subdiag[start];
-
   for (Index k = start; k < end; ++k)
   {
     JacobiRotation<RealScalar> rot;
@@ -489,6 +493,7 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
     diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
     diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
     subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
+    
 
     if (k > start)
       subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
@@ -500,7 +505,7 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
       z = -rot.s() * subdiag[k+1];
       subdiag[k + 1] = rot.c() * subdiag[k+1];
     }
-
+    
     // apply the givens rotation to the unit matrix Q = Q * G
     if (matrixQ)
     {

Eigen/src/Geometry/AngleAxis.h

@@ -171,6 +171,9 @@ template<typename Scalar>
 template<typename QuatDerived>
 AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived>& q)
 {
+  using std::acos;
+  using std::min;
+  using std::max;
   Scalar n2 = q.vec().squaredNorm();
   if (n2 < NumTraits<Scalar>::dummy_precision()*NumTraits<Scalar>::dummy_precision())
   {
@@ -179,7 +182,7 @@ AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived
   }
   else
   {
-    m_angle = Scalar(2)*std::acos(std::min(std::max(Scalar(-1),q.w()),Scalar(1)));
+    m_angle = Scalar(2)*acos(min(max(Scalar(-1),q.w()),Scalar(1)));
     m_axis = q.vec() / internal::sqrt(n2);
   }
   return *this;

Eigen/src/Geometry/Homogeneous.h

@@ -232,13 +232,15 @@ template<typename MatrixType,typename Lhs>
 struct traits<homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs> >
 {
   typedef typename take_matrix_for_product<Lhs>::type LhsMatrixType;
+  typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
+  typedef typename remove_all<LhsMatrixType>::type LhsMatrixTypeCleaned;
   typedef typename make_proper_matrix_type<
-                 typename traits<MatrixType>::Scalar,
-                 LhsMatrixType::RowsAtCompileTime,
-                 MatrixType::ColsAtCompileTime,
-                 MatrixType::PlainObject::Options,
-                 LhsMatrixType::MaxRowsAtCompileTime,
-                 MatrixType::MaxColsAtCompileTime>::type ReturnType;
+                 typename traits<MatrixTypeCleaned>::Scalar,
+                 LhsMatrixTypeCleaned::RowsAtCompileTime,
+                 MatrixTypeCleaned::ColsAtCompileTime,
+                 MatrixTypeCleaned::PlainObject::Options,
+                 LhsMatrixTypeCleaned::MaxRowsAtCompileTime,
+                 MatrixTypeCleaned::MaxColsAtCompileTime>::type ReturnType;
 };
 
 template<typename MatrixType,typename Lhs>
@@ -246,7 +248,8 @@ struct homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs>
   : public ReturnByValue<homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs> >
 {
   typedef typename traits<homogeneous_left_product_impl>::LhsMatrixType LhsMatrixType;
-  typedef typename remove_all<typename LhsMatrixType::Nested>::type LhsMatrixTypeNested;
+  typedef typename remove_all<LhsMatrixType>::type LhsMatrixTypeCleaned;
+  typedef typename remove_all<typename LhsMatrixTypeCleaned::Nested>::type LhsMatrixTypeNested;
   typedef typename MatrixType::Index Index;
   homogeneous_left_product_impl(const Lhs& lhs, const MatrixType& rhs)
     : m_lhs(take_matrix_for_product<Lhs>::run(lhs)),
@@ -267,7 +270,7 @@ struct homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs>
             .template replicate<MatrixType::ColsAtCompileTime>(m_rhs.cols());
   }
 
-  const typename LhsMatrixType::Nested m_lhs;
+  const typename LhsMatrixTypeCleaned::Nested m_lhs;
   const typename MatrixType::Nested m_rhs;
 };
 

Eigen/src/Geometry/Hyperplane.h

@@ -213,8 +213,8 @@ public:
   /** Applies the transformation matrix \a mat to \c *this and returns a reference to \c *this.
     *
     * \param mat the Dim x Dim transformation matrix
-    * \param traits specifies whether the matrix \a mat represents an Isometry
-    *               or a more generic Affine transformation. The default is Affine.
+    * \param traits specifies whether the matrix \a mat represents an #Isometry
+    *               or a more generic #Affine transformation. The default is #Affine.
     */
   template<typename XprType>
   inline Hyperplane& transform(const MatrixBase<XprType>& mat, TransformTraits traits = Affine)
@@ -233,8 +233,8 @@ public:
   /** Applies the transformation \a t to \c *this and returns a reference to \c *this.
     *
     * \param t the transformation of dimension Dim
-    * \param traits specifies whether the transformation \a t represents an Isometry
-    *               or a more generic Affine transformation. The default is Affine.
+    * \param traits specifies whether the transformation \a t represents an #Isometry
+    *               or a more generic #Affine transformation. The default is #Affine.
     *               Other kind of transformations are not supported.
     */
   template<int TrOptions>

Eigen/src/Geometry/Quaternion.h

@@ -49,6 +49,9 @@ public:
   typedef typename internal::traits<Derived>::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef typename internal::traits<Derived>::Coefficients Coefficients;
+  enum {
+    Flags = Eigen::internal::traits<Derived>::Flags
+  };
 
  // typedef typename Matrix<Scalar,4,1> Coefficients;
   /** the type of a 3D vector */
@@ -222,7 +225,8 @@ struct traits<Quaternion<_Scalar,_Options> >
   typedef _Scalar Scalar;
   typedef Matrix<_Scalar,4,1,_Options> Coefficients;
   enum{
-    PacketAccess = _Options & DontAlign ? Unaligned : Aligned
+    IsAligned = bool(EIGEN_ALIGN) && ((int(_Options)&Aligned)==Aligned),
+    Flags = IsAligned ? (AlignedBit | LvalueBit) : LvalueBit
   };
 };
 }
@@ -294,33 +298,53 @@ typedef Quaternion<double> Quaterniond;
 ***************************************************************************/
 
 namespace internal {
-template<typename _Scalar, int _PacketAccess>
-struct traits<Map<Quaternion<_Scalar>, _PacketAccess> >:
-traits<Quaternion<_Scalar> >
-{
-  typedef _Scalar Scalar;
-  typedef Map<Matrix<_Scalar,4,1>, _PacketAccess> Coefficients;
-  enum {
-    PacketAccess = _PacketAccess
+  template<typename _Scalar, int _Options>
+  struct traits<Map<Quaternion<_Scalar>, _Options> >:
+  traits<Quaternion<_Scalar, _Options> >
+  {
+    typedef _Scalar Scalar;
+    typedef Map<Matrix<_Scalar,4,1>, _Options> Coefficients;
+
+    typedef traits<Quaternion<_Scalar, _Options> > TraitsBase;
+    enum {
+      IsAligned = TraitsBase::IsAligned,
+
+      Flags = TraitsBase::Flags
+    };
+  };
+}
+
+namespace internal {
+  template<typename _Scalar, int _Options>
+  struct traits<Map<const Quaternion<_Scalar>, _Options> >:
+  traits<Quaternion<_Scalar> >
+  {
+    typedef _Scalar Scalar;
+    typedef Map<const Matrix<_Scalar,4,1>, _Options> Coefficients;
+
+    typedef traits<Quaternion<_Scalar, _Options> > TraitsBase;
+    enum {
+      IsAligned = TraitsBase::IsAligned,
+      Flags = TraitsBase::Flags & ~LvalueBit
+    };
   };
-};
 }
 
 /** \brief Quaternion expression mapping a constant memory buffer
   *
   * \param _Scalar the type of the Quaternion coefficients
-  * \param PacketAccess see class Map
+  * \param _Options see class Map
   *
   * This is a specialization of class Map for Quaternion. This class allows to view
   * a 4 scalar memory buffer as an Eigen's Quaternion object.
   *
   * \sa class Map, class Quaternion, class QuaternionBase
   */
-template<typename _Scalar, int PacketAccess>
-class Map<const Quaternion<_Scalar>, PacketAccess >
-  : public QuaternionBase<Map<const Quaternion<_Scalar>, PacketAccess> >
+template<typename _Scalar, int _Options>
+class Map<const Quaternion<_Scalar>, _Options >
+  : public QuaternionBase<Map<const Quaternion<_Scalar>, _Options> >
 {
-    typedef QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> > Base;
+    typedef QuaternionBase<Map<const Quaternion<_Scalar>, _Options> > Base;
 
   public:
     typedef _Scalar Scalar;
@@ -333,7 +357,7 @@ class Map<const Quaternion<_Scalar>, PacketAccess >
       * The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order:
       * \code *coeffs == {x, y, z, w} \endcode
       *
-      * If the template parameter PacketAccess is set to Aligned, then the pointer coeffs must be aligned. */
+      * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */
     EIGEN_STRONG_INLINE Map(const Scalar* coeffs) : m_coeffs(coeffs) {}
 
     inline const Coefficients& coeffs() const { return m_coeffs;}
@@ -345,18 +369,18 @@ class Map<const Quaternion<_Scalar>, PacketAccess >
 /** \brief Expression of a quaternion from a memory buffer
   *
   * \param _Scalar the type of the Quaternion coefficients
-  * \param PacketAccess see class Map
+  * \param _Options see class Map
   *
   * This is a specialization of class Map for Quaternion. This class allows to view
   * a 4 scalar memory buffer as an Eigen's  Quaternion object.
   *
   * \sa class Map, class Quaternion, class QuaternionBase
   */
-template<typename _Scalar, int PacketAccess>
-class Map<Quaternion<_Scalar>, PacketAccess >
-  : public QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> >
+template<typename _Scalar, int _Options>
+class Map<Quaternion<_Scalar>, _Options >
+  : public QuaternionBase<Map<Quaternion<_Scalar>, _Options> >
 {
-    typedef QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> > Base;
+    typedef QuaternionBase<Map<Quaternion<_Scalar>, _Options> > Base;
 
   public:
     typedef _Scalar Scalar;
@@ -369,7 +393,7 @@ class Map<Quaternion<_Scalar>, PacketAccess >
       * The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order:
       * \code *coeffs == {x, y, z, w} \endcode
       *
-      * If the template parameter PacketAccess is set to Aligned, then the pointer coeffs must be aligned. */
+      * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */
     EIGEN_STRONG_INLINE Map(Scalar* coeffs) : m_coeffs(coeffs) {}
 
     inline Coefficients& coeffs() { return m_coeffs; }
@@ -399,7 +423,7 @@ typedef Map<Quaternion<double>, Aligned>  QuaternionMapAlignedd;
 // Generic Quaternion * Quaternion product
 // This product can be specialized for a given architecture via the Arch template argument.
 namespace internal {
-template<int Arch, class Derived1, class Derived2, typename Scalar, int PacketAccess> struct quat_product
+template<int Arch, class Derived1, class Derived2, typename Scalar, int _Options> struct quat_product
 {
   EIGEN_STRONG_INLINE static Quaternion<Scalar> run(const QuaternionBase<Derived1>& a, const QuaternionBase<Derived2>& b){
     return Quaternion<Scalar>
@@ -423,7 +447,7 @@ QuaternionBase<Derived>::operator* (const QuaternionBase<OtherDerived>& other) c
    YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
   return internal::quat_product<Architecture::Target, Derived, OtherDerived,
                          typename internal::traits<Derived>::Scalar,
-                         internal::traits<Derived>::PacketAccess && internal::traits<OtherDerived>::PacketAccess>::run(*this, other);
+                         internal::traits<Derived>::IsAligned && internal::traits<OtherDerived>::IsAligned>::run(*this, other);
 }
 
 /** \sa operator*(Quaternion) */
@@ -551,6 +575,7 @@ template<class Derived>
 template<typename Derived1, typename Derived2>
 inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b)
 {
+  using std::max;
   Vector3 v0 = a.normalized();
   Vector3 v1 = b.normalized();
   Scalar c = v1.dot(v0);
@@ -565,7 +590,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri
   //    which yields a singular value problem
   if (c < Scalar(-1)+NumTraits<Scalar>::dummy_precision())
   {
-    c = std::max<Scalar>(c,-1);
+    c = max<Scalar>(c,-1);
     Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
     JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
     Vector3 axis = svd.matrixV().col(2);
@@ -625,10 +650,11 @@ template <class OtherDerived>
 inline typename internal::traits<Derived>::Scalar
 QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& other) const
 {
+  using std::acos;
   double d = internal::abs(this->dot(other));
   if (d>=1.0)
     return Scalar(0);
-  return static_cast<Scalar>(2 * std::acos(d));
+  return static_cast<Scalar>(2 * acos(d));
 }
 
 /** \returns the spherical linear interpolation between the two quaternions
@@ -639,6 +665,7 @@ template <class OtherDerived>
 Quaternion<typename internal::traits<Derived>::Scalar>
 QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const
 {
+  using std::acos;
   static const Scalar one = Scalar(1) - NumTraits<Scalar>::epsilon();
   Scalar d = this->dot(other);
   Scalar absD = internal::abs(d);
@@ -654,7 +681,7 @@ QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& oth
   else
   {
     // theta is the angle between the 2 quaternions
-    Scalar theta = std::acos(absD);
+    Scalar theta = acos(absD);
     Scalar sinTheta = internal::sin(theta);
 
     scale0 = internal::sin( ( Scalar(1) - t ) * theta) / sinTheta;

Eigen/src/Geometry/Transform.h

@@ -43,7 +43,9 @@ struct transform_traits
 
 template< typename TransformType,
           typename MatrixType,
-          bool IsProjective = transform_traits<TransformType>::IsProjective>
+          int Case = transform_traits<TransformType>::IsProjective ? 0
+                   : int(MatrixType::RowsAtCompileTime) == int(transform_traits<TransformType>::HDim) ? 1
+                   : 2>
 struct transform_right_product_impl;
 
 template< typename Other,
@@ -81,15 +83,16 @@ template<typename TransformType> struct transform_take_affine_part;
   *
   * \brief Represents an homogeneous transformation in a N dimensional space
   *
-  * \param _Scalar the scalar type, i.e., the type of the coefficients
-  * \param _Dim the dimension of the space
-  * \param _Mode the type of the transformation. Can be:
-  *              - Affine: the transformation is stored as a (Dim+1)^2 matrix,
-  *                        where the last row is assumed to be [0 ... 0 1].
-  *              - AffineCompact: the transformation is stored as a (Dim)x(Dim+1) matrix.
-  *              - Projective: the transformation is stored as a (Dim+1)^2 matrix
-  *                            without any assumption.
-  * \param _Options can be \b AutoAlign or \b DontAlign. Default is \b AutoAlign
+  * \tparam _Scalar the scalar type, i.e., the type of the coefficients
+  * \tparam _Dim the dimension of the space
+  * \tparam _Mode the type of the transformation. Can be:
+  *              - #Affine: the transformation is stored as a (Dim+1)^2 matrix,
+  *                         where the last row is assumed to be [0 ... 0 1].
+  *              - #AffineCompact: the transformation is stored as a (Dim)x(Dim+1) matrix.
+  *              - #Projective: the transformation is stored as a (Dim+1)^2 matrix
+  *                             without any assumption.
+  * \tparam _Options has the same meaning as in class Matrix. It allows to specify DontAlign and/or RowMajor.
+  *                  These Options are passed directly to the underlying matrix type.
   *
   * The homography is internally represented and stored by a matrix which
   * is available through the matrix() method. To understand the behavior of
@@ -177,6 +180,9 @@ template<typename TransformType> struct transform_take_affine_part;
   * Conversion methods from/to Qt's QMatrix and QTransform are available if the
   * preprocessor token EIGEN_QT_SUPPORT is defined.
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_TRANSFORM_PLUGIN.
+  *
   * \sa class Matrix, class Quaternion
   */
 template<typename _Scalar, int _Dim, int _Mode, int _Options>
@@ -195,11 +201,11 @@ public:
   typedef _Scalar Scalar;
   typedef DenseIndex Index;
   /** type of the matrix used to represent the transformation */
-  typedef Matrix<Scalar,Rows,HDim,Options&DontAlign> MatrixType;
+  typedef typename internal::make_proper_matrix_type<Scalar,Rows,HDim,Options>::type MatrixType;
   /** constified MatrixType */
   typedef const MatrixType ConstMatrixType;
   /** type of the matrix used to represent the linear part of the transformation */
-  typedef Matrix<Scalar,Dim,Dim> LinearMatrixType;
+  typedef Matrix<Scalar,Dim,Dim,Options> LinearMatrixType;
   /** type of read/write reference to the linear part of the transformation */
   typedef Block<MatrixType,Dim,Dim> LinearPart;
   /** type of read reference to the linear part of the transformation */
@@ -517,7 +523,7 @@ public:
   template<typename Derived>
   inline Transform operator*(const RotationBase<Derived,Dim>& r) const;
 
-  LinearMatrixType rotation() const;
+  const LinearMatrixType rotation() const;
   template<typename RotationMatrixType, typename ScalingMatrixType>
   void computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const;
   template<typename ScalingMatrixType, typename RotationMatrixType>
@@ -604,7 +610,7 @@ protected:
   #ifndef EIGEN_PARSED_BY_DOXYGEN
     EIGEN_STRONG_INLINE static void check_template_params()
     {
-      EIGEN_STATIC_ASSERT((Options & (DontAlign)) == Options, INVALID_MATRIX_TEMPLATE_PARAMETERS)
+      EIGEN_STATIC_ASSERT((Options & (DontAlign|RowMajor)) == Options, INVALID_MATRIX_TEMPLATE_PARAMETERS)
     }
   #endif
 
@@ -666,7 +672,7 @@ Transform<Scalar,Dim,Mode,Options>::Transform(const QMatrix& other)
   *
   * This function is available only if the token EIGEN_QT_SUPPORT is defined.
   */
-template<typename Scalar, int Dim, int Mode,int Otpions>
+template<typename Scalar, int Dim, int Mode,int Options>
 Transform<Scalar,Dim,Mode,Options>& Transform<Scalar,Dim,Mode,Options>::operator=(const QMatrix& other)
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
@@ -712,9 +718,13 @@ Transform<Scalar,Dim,Mode,Options>& Transform<Scalar,Dim,Mode,Options>::operator
 {
   check_template_params();
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  m_matrix << other.m11(), other.m21(), other.dx(),
-              other.m12(), other.m22(), other.dy(),
-              other.m13(), other.m23(), other.m33();
+  if (Mode == int(AffineCompact))
+    m_matrix << other.m11(), other.m21(), other.dx(),
+                other.m12(), other.m22(), other.dy();
+  else
+    m_matrix << other.m11(), other.m21(), other.dx(),
+                other.m12(), other.m22(), other.dy(),
+                other.m13(), other.m23(), other.m33();
   return *this;
 }
 
@@ -726,9 +736,14 @@ template<typename Scalar, int Dim, int Mode, int Options>
 QTransform Transform<Scalar,Dim,Mode,Options>::toQTransform(void) const
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  return QTransform(matrix.coeff(0,0), matrix.coeff(1,0), matrix.coeff(2,0)
-                    matrix.coeff(0,1), matrix.coeff(1,1), matrix.coeff(2,1)
-                    matrix.coeff(0,2), matrix.coeff(1,2), matrix.coeff(2,2));
+  if (Mode == int(AffineCompact))
+    return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0),
+                      m_matrix.coeff(0,1), m_matrix.coeff(1,1),
+                      m_matrix.coeff(0,2), m_matrix.coeff(1,2));
+  else
+    return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0), m_matrix.coeff(2,0),
+                      m_matrix.coeff(0,1), m_matrix.coeff(1,1), m_matrix.coeff(2,1),
+                      m_matrix.coeff(0,2), m_matrix.coeff(1,2), m_matrix.coeff(2,2));
 }
 #endif
 
@@ -964,7 +979,7 @@ inline Transform<Scalar,Dim,Mode,Options> Transform<Scalar,Dim,Mode,Options>::op
   * \sa computeRotationScaling(), computeScalingRotation(), class SVD
   */
 template<typename Scalar, int Dim, int Mode, int Options>
-typename Transform<Scalar,Dim,Mode,Options>::LinearMatrixType
+const typename Transform<Scalar,Dim,Mode,Options>::LinearMatrixType
 Transform<Scalar,Dim,Mode,Options>::rotation() const
 {
   LinearMatrixType result;
@@ -1076,10 +1091,10 @@ struct projective_transform_inverse<TransformType, Projective>
   *
   * \param hint allows to optimize the inversion process when the transformation
   * is known to be not a general transformation (optional). The possible values are:
-  *  - Projective if the transformation is not necessarily affine, i.e., if the
+  *  - #Projective if the transformation is not necessarily affine, i.e., if the
   *    last row is not guaranteed to be [0 ... 0 1]
-  *  - Affine if the last row can be assumed to be [0 ... 0 1]
-  *  - Isometry if the transformation is only a concatenations of translations
+  *  - #Affine if the last row can be assumed to be [0 ... 0 1]
+  *  - #Isometry if the transformation is only a concatenations of translations
   *    and rotations.
   *  The default is the template class parameter \c Mode.
   *
@@ -1136,9 +1151,9 @@ template<typename TransformType> struct transform_take_affine_part {
   { return m.template block<TransformType::Dim,TransformType::HDim>(0,0); }
 };
 
-template<typename Scalar, int Dim>
-struct transform_take_affine_part<Transform<Scalar,Dim,AffineCompact> > {
-  typedef typename Transform<Scalar,Dim,AffineCompact>::MatrixType MatrixType;
+template<typename Scalar, int Dim, int Options>
+struct transform_take_affine_part<Transform<Scalar,Dim,AffineCompact, Options> > {
+  typedef typename Transform<Scalar,Dim,AffineCompact,Options>::MatrixType MatrixType;
   static inline MatrixType& run(MatrixType& m) { return m; }
   static inline const MatrixType& run(const MatrixType& m) { return m; }
 };
@@ -1178,7 +1193,7 @@ struct transform_construct_from_matrix<Other, Mode,Options,Dim,HDim, HDim,HDim>
 template<typename Other, int Options, int Dim, int HDim>
 struct transform_construct_from_matrix<Other, AffineCompact,Options,Dim,HDim, HDim,HDim>
 {
-  static inline void run(Transform<typename Other::Scalar,Dim,AffineCompact> *transform, const Other& other)
+  static inline void run(Transform<typename Other::Scalar,Dim,AffineCompact,Options> *transform, const Other& other)
   { transform->matrix() = other.template block<Dim,HDim>(0,0); }
 };
 
@@ -1200,7 +1215,7 @@ struct transform_product_result
 };
 
 template< typename TransformType, typename MatrixType >
-struct transform_right_product_impl< TransformType, MatrixType, true >
+struct transform_right_product_impl< TransformType, MatrixType, 0 >
 {
   typedef typename MatrixType::PlainObject ResultType;
 
@@ -1211,7 +1226,7 @@ struct transform_right_product_impl< TransformType, MatrixType, true >
 };
 
 template< typename TransformType, typename MatrixType >
-struct transform_right_product_impl< TransformType, MatrixType, false >
+struct transform_right_product_impl< TransformType, MatrixType, 1 >
 {
   enum { 
     Dim = TransformType::Dim, 
@@ -1224,20 +1239,39 @@ struct transform_right_product_impl< TransformType, MatrixType, false >
 
   EIGEN_STRONG_INLINE static ResultType run(const TransformType& T, const MatrixType& other)
   {
-    EIGEN_STATIC_ASSERT(OtherRows==Dim || OtherRows==HDim, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES);
+    EIGEN_STATIC_ASSERT(OtherRows==HDim, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES);
 
     typedef Block<ResultType, Dim, OtherCols> TopLeftLhs;
-    typedef Block<const MatrixType, Dim, OtherCols> TopLeftRhs;
 
     ResultType res(other.rows(),other.cols());
+    TopLeftLhs(res, 0, 0, Dim, other.cols()).noalias() = T.affine() * other;
+    res.row(OtherRows-1) = other.row(OtherRows-1);
+    
+    return res;
+  }
+};
 
-    TopLeftLhs(res, 0, 0, Dim, other.cols()) =
-      ( T.linear() * TopLeftRhs(other, 0, 0, Dim, other.cols()) ).colwise() +
-        T.translation();
+template< typename TransformType, typename MatrixType >
+struct transform_right_product_impl< TransformType, MatrixType, 2 >
+{
+  enum { 
+    Dim = TransformType::Dim, 
+    HDim = TransformType::HDim,
+    OtherRows = MatrixType::RowsAtCompileTime,
+    OtherCols = MatrixType::ColsAtCompileTime
+  };
 
-    // we need to take .rows() because OtherRows might be Dim or HDim
-    if (OtherRows==HDim)
-      res.row(other.rows()) = other.row(other.rows());
+  typedef typename MatrixType::PlainObject ResultType;
+
+  EIGEN_STRONG_INLINE static ResultType run(const TransformType& T, const MatrixType& other)
+  {
+    EIGEN_STATIC_ASSERT(OtherRows==Dim, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES);
+
+    typedef Block<ResultType, Dim, OtherCols> TopLeftLhs;
+
+    ResultType res(other.rows(),other.cols());
+    TopLeftLhs(res, 0, 0, Dim, other.cols()).noalias() = T.linear() * other;
+    TopLeftLhs(res, 0, 0, Dim, other.cols()).colwise() += T.translation();
 
     return res;
   }

Eigen/src/Householder/Householder.h

@@ -72,8 +72,9 @@ void MatrixBase<Derived>::makeHouseholder(
 
   if(tailSqNorm == RealScalar(0) && internal::imag(c0)==RealScalar(0))
   {
-    tau = 0;
+    tau = RealScalar(0);
     beta = internal::real(c0);
+    essential.setZero();
   }
   else
   {

Eigen/src/Jacobi/Jacobi.h

@@ -104,9 +104,9 @@ bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z)
   else
   {
     RealScalar tau = (x-z)/(RealScalar(2)*internal::abs(y));
-    RealScalar w = internal::sqrt(internal::abs2(tau) + 1);
+    RealScalar w = internal::sqrt(internal::abs2(tau) + RealScalar(1));
     RealScalar t;
-    if(tau>0)
+    if(tau>RealScalar(0))
     {
       t = RealScalar(1) / (tau + w);
     }
@@ -114,8 +114,8 @@ bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z)
     {
       t = RealScalar(1) / (tau - w);
     }
-    RealScalar sign_t = t > 0 ? 1 : -1;
-    RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+1);
+    RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
+    RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+RealScalar(1));
     m_s = - sign_t * (internal::conj(y) / internal::abs(y)) * internal::abs(t) * n;
     m_c = n;
     return true;
@@ -221,15 +221,15 @@ template<typename Scalar>
 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type)
 {
 
-  if(q==0)
+  if(q==Scalar(0))
   {
     m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
-    m_s = 0;
+    m_s = Scalar(0);
     if(r) *r = internal::abs(p);
   }
-  else if(p==0)
+  else if(p==Scalar(0))
   {
-    m_c = 0;
+    m_c = Scalar(0);
     m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
     if(r) *r = internal::abs(q);
   }

Eigen/src/LU/PartialPivLU.h

@@ -268,7 +268,7 @@ struct partial_lu_impl
 
       row_transpositions[k] = row_of_biggest_in_col;
 
-      if(biggest_in_corner != 0)
+      if(biggest_in_corner != RealScalar(0))
       {
         if(k != row_of_biggest_in_col)
         {

Eigen/src/QR/ColPivHouseholderQR.h

@@ -330,12 +330,12 @@ template<typename _MatrixType> class ColPivHouseholderQR
       */
     inline Index nonzeroPivots() const
     {
-      eigen_assert(m_isInitialized && "LU is not initialized.");
+      eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_nonzero_pivots;
     }
 
     /** \returns the absolute value of the biggest pivot, i.e. the biggest
-      *          diagonal coefficient of U.
+      *          diagonal coefficient of R.
       */
     RealScalar maxPivot() const { return m_maxpivot; }
 
@@ -387,7 +387,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   for(Index k = 0; k < cols; ++k)
     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
 
-  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / rows;
+  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
@@ -413,7 +413,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
     // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
     // repetitions of the unit test, with the result of solve() filled with large values of the order
     // of 1/(size*epsilon).
-    if(biggest_col_sq_norm < threshold_helper * (rows-k))
+    if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
     {
       m_nonzero_pivots = k;
       m_hCoeffs.tail(size-k).setZero();

Eigen/src/SVD/JacobiSVD.h

@@ -271,13 +271,13 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
   RealScalar d = m.coeff(1,0) - m.coeff(0,1);
   if(t == RealScalar(0))
   {
-    rot1.c() = 0;
-    rot1.s() = d > 0 ? 1 : -1;
+    rot1.c() = RealScalar(0);
+    rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
   }
   else
   {
     RealScalar u = d / t;
-    rot1.c() = RealScalar(1) / sqrt(1 + abs2(u));
+    rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u));
     rot1.s() = rot1.c() * u;
   }
   m.applyOnTheLeft(0,1,rot1);
@@ -292,7 +292,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
   *
   * \class JacobiSVD
   *
-  * \brief Two-sided Jacobi SVD decomposition of a square matrix
+  * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix
   *
   * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
   * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally
@@ -376,7 +376,12 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
       * The default constructor is useful in cases in which the user intends to
       * perform decompositions via JacobiSVD::compute(const MatrixType&).
       */
-    JacobiSVD() : m_isInitialized(false) {}
+    JacobiSVD()
+      : m_isInitialized(false),
+        m_isAllocated(false),
+        m_computationOptions(0),
+        m_rows(-1), m_cols(-1)
+    {}
 
 
     /** \brief Default Constructor with memory preallocation
@@ -386,6 +391,10 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
       * \sa JacobiSVD()
       */
     JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
+      : m_isInitialized(false),
+        m_isAllocated(false),
+        m_computationOptions(0),
+        m_rows(-1), m_cols(-1)
     {
       allocate(rows, cols, computationOptions);
     }
@@ -394,33 +403,48 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
      *
      * \param matrix the matrix to decompose
      * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
-     *                           By default, none is computed. This is a bit-field, the possible bits are ComputeFullU, ComputeThinU,
-     *                           ComputeFullV, ComputeThinV.
+     *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
+     *                           #ComputeFullV, #ComputeThinV.
      *
      * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
      * available with the (non-default) FullPivHouseholderQR preconditioner.
      */
     JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
+      : m_isInitialized(false),
+        m_isAllocated(false),
+        m_computationOptions(0),
+        m_rows(-1), m_cols(-1)
     {
       compute(matrix, computationOptions);
     }
 
-    /** \brief Method performing the decomposition of given matrix.
+    /** \brief Method performing the decomposition of given matrix using custom options.
      *
      * \param matrix the matrix to decompose
      * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
-     *                           By default, none is computed. This is a bit-field, the possible bits are ComputeFullU, ComputeThinU,
-     *                           ComputeFullV, ComputeThinV.
+     *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
+     *                           #ComputeFullV, #ComputeThinV.
      *
      * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
      * available with the (non-default) FullPivHouseholderQR preconditioner.
      */
-    JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0);
+    JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
+
+    /** \brief Method performing the decomposition of given matrix using current options.
+     *
+     * \param matrix the matrix to decompose
+     *
+     * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int).
+     */
+    JacobiSVD& compute(const MatrixType& matrix)
+    {
+      return compute(matrix, m_computationOptions);
+    }
 
     /** \returns the \a U matrix.
      *
      * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
-     * the U matrix is n-by-n if you asked for ComputeFullU, and is n-by-m if you asked for ComputeThinU.
+     * the U matrix is n-by-n if you asked for #ComputeFullU, and is n-by-m if you asked for #ComputeThinU.
      *
      * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed.
      *
@@ -436,7 +460,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     /** \returns the \a V matrix.
      *
      * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
-     * the V matrix is p-by-p if you asked for ComputeFullV, and is p-by-m if you asked for ComputeThinV.
+     * the V matrix is p-by-p if you asked for #ComputeFullV, and is p-by-m if you asked for ComputeThinV.
      *
      * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed.
      *
@@ -452,7 +476,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     /** \returns the vector of singular values.
      *
      * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the
-     * returned vector has size \a m.
+     * returned vector has size \a m.  Singular values are always sorted in decreasing order.
      */
     const SingularValuesType& singularValues() const
     {
@@ -494,16 +518,17 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     inline Index cols() const { return m_cols; }
 
   private:
-    void allocate(Index rows, Index cols, unsigned int computationOptions = 0);
+    void allocate(Index rows, Index cols, unsigned int computationOptions);
 
   protected:
     MatrixUType m_matrixU;
     MatrixVType m_matrixV;
     SingularValuesType m_singularValues;
     WorkMatrixType m_workMatrix;
-    bool m_isInitialized;
+    bool m_isInitialized, m_isAllocated;
     bool m_computeFullU, m_computeThinU;
     bool m_computeFullV, m_computeThinV;
+    unsigned int m_computationOptions;
     Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
 
     template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
@@ -515,9 +540,21 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
 template<typename MatrixType, int QRPreconditioner>
 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
 {
+  eigen_assert(rows >= 0 && cols >= 0);
+
+  if (m_isAllocated &&
+      rows == m_rows &&
+      cols == m_cols &&
+      computationOptions == m_computationOptions)
+  {
+    return;
+  }
+
   m_rows = rows;
   m_cols = cols;
   m_isInitialized = false;
+  m_isAllocated = true;
+  m_computationOptions = computationOptions;
   m_computeFullU = (computationOptions & ComputeFullU) != 0;
   m_computeThinU = (computationOptions & ComputeThinU) != 0;
   m_computeFullV = (computationOptions & ComputeFullV) != 0;
@@ -581,8 +618,9 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
         // if this 2x2 sub-matrix is not diagonal already...
         // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
         // keep us iterating forever.
-        if(std::max(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p)))
-            > std::max(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision)
+        using std::max;
+        if(max(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p)))
+            > max(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision)
         {
           finished = false;
 

Eigen/src/Sparse/SparseMatrix.h

@@ -34,13 +34,15 @@
   * This class implements a sparse matrix using the very common compressed row/column storage
   * scheme.
   *
-  * \param _Scalar the scalar type, i.e. the type of the coefficients
-  * \param _Options Union of bit flags controlling the storage scheme. Currently the only possibility
+  * \tparam _Scalar the scalar type, i.e. the type of the coefficients
+  * \tparam _Options Union of bit flags controlling the storage scheme. Currently the only possibility
   *                 is RowMajor. The default is 0 which means column-major.
-  * \param _Index the type of the indices. Default is \c int.
+  * \tparam _Index the type of the indices. Default is \c int.
   *
   * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme.
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_SPARSEMATRIX_PLUGIN.
   */
 
 namespace internal {

Eigen/src/Sparse/SparseMatrixBase.h

@@ -31,10 +31,10 @@
   *
   * \brief Base class of any sparse matrices or sparse expressions
   *
-  * \param Derived
-  *
-  *
+  * \tparam Derived
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_SPARSEMATRIXBASE_PLUGIN.
   */
 template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
 {

Eigen/src/Sparse/SparseSelfAdjointView.h

@@ -31,13 +31,13 @@
   * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
   *
   * \param MatrixType the type of the dense matrix storing the coefficients
-  * \param UpLo can be either \c Lower or \c Upper
+  * \param UpLo can be either \c #Lower or \c #Upper
   *
   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
   * and most of the time this is the only way that it is used.
   *
-  * \sa SparseMatrixBase::selfAdjointView()
+  * \sa SparseMatrixBase::selfadjointView()
   */
 template<typename Lhs, typename Rhs, int UpLo>
 class SparseSelfAdjointTimeDenseProduct;

Eigen/src/Sparse/SparseSparseProduct.h

@@ -133,7 +133,12 @@ static void sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
   float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
 
-  res.resize(rows, cols);
+  // mimics a resizeByInnerOuter:
+  if(ResultType::IsRowMajor)
+    res.resize(cols, rows);
+  else
+    res.resize(rows, cols);
+
   res.reserve(Index(ratioRes*rows*cols));
   for (Index j=0; j<cols; ++j)
   {

Eigen/src/Sparse/SparseVector.h

@@ -29,10 +29,12 @@
   *
   * \brief a sparse vector class
   *
-  * \param _Scalar the scalar type, i.e. the type of the coefficients
+  * \tparam _Scalar the scalar type, i.e. the type of the coefficients
   *
   * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme.
   *
+  * This class can be extended with the help of the plugin mechanism described on the page
+  * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_SPARSEVECTOR_PLUGIN.
   */
 
 namespace internal {

Eigen/src/plugins/ArrayCwiseUnaryOps.h

@@ -70,7 +70,7 @@ sqrt() const
   * Example: \include Cwise_cos.cpp
   * Output: \verbinclude Cwise_cos.out
   *
-  * \sa sin(), exp()
+  * \sa sin(), acos()
   */
 inline const CwiseUnaryOp<internal::scalar_cos_op<Scalar>, const Derived>
 cos() const
@@ -84,7 +84,7 @@ cos() const
   * Example: \include Cwise_sin.cpp
   * Output: \verbinclude Cwise_sin.out
   *
-  * \sa cos(), exp()
+  * \sa cos(), asin()
   */
 inline const CwiseUnaryOp<internal::scalar_sin_op<Scalar>, const Derived>
 sin() const
@@ -92,6 +92,31 @@ sin() const
   return derived();
 }
 
+/** \returns an expression of the coefficient-wise arc cosine of *this.
+  *
+  * Example: \include Cwise_acos.cpp
+  * Output: \verbinclude Cwise_acos.out
+  *
+  * \sa cos(), asin()
+  */
+inline const CwiseUnaryOp<internal::scalar_acos_op<Scalar>, const Derived>
+acos() const
+{
+  return derived();
+}
+
+/** \returns an expression of the coefficient-wise arc sine of *this.
+  *
+  * Example: \include Cwise_asin.cpp
+  * Output: \verbinclude Cwise_asin.out
+  *
+  * \sa sin(), acos()
+  */
+inline const CwiseUnaryOp<internal::scalar_asin_op<Scalar>, const Derived>
+asin() const
+{
+  return derived();
+}
 
 /** \returns an expression of the coefficient-wise tan of *this.
   *

bench/BenchTimer.h

@@ -119,7 +119,7 @@ public:
 
   inline double getCpuTime()
   {
-#ifdef WIN32
+#ifdef _WIN32
     LARGE_INTEGER query_ticks;
     QueryPerformanceCounter(&query_ticks);
     return query_ticks.QuadPart/m_frequency;
@@ -132,7 +132,7 @@ public:
 
   inline double getRealTime()
   {
-#ifdef WIN32
+#ifdef _WIN32
     SYSTEMTIME st;
     GetSystemTime(&st);
     return (double)st.wSecond + 1.e-3 * (double)st.wMilliseconds;

bench/btl/libs/C_BLAS/blas.h

@@ -46,6 +46,11 @@ double  BLASFUNC(xdotu)  (int *, double  *, int *, double  *, int *);
 double  BLASFUNC(xdotc)  (int *, double  *, int *, double  *, int *);
 #endif
 
+int  BLASFUNC(cdotuw)  (int *, float  *, int *, float  *, int *, float*);
+int  BLASFUNC(cdotcw)  (int *, float  *, int *, float  *, int *, float*);
+int  BLASFUNC(zdotuw)  (int *, double  *, int *, double  *, int *, double*);
+int  BLASFUNC(zdotcw)  (int *, double  *, int *, double  *, int *, double*);
+
 int    BLASFUNC(saxpy) (int *, float  *, float  *, int *, float  *, int *);
 int    BLASFUNC(daxpy) (int *, double *, double *, int *, double *, int *);
 int    BLASFUNC(qaxpy) (int *, double *, double *, int *, double *, int *);

blas/CMakeLists.txt

@@ -18,6 +18,7 @@ if(CMAKE_Fortran_COMPILER_WORKS)
 add_custom_target(blas)
 
 set(EigenBlas_SRCS single.cpp double.cpp complex_single.cpp complex_double.cpp xerbla.cpp
+    complexdots.f
     srotm.f srotmg.f drotm.f drotmg.f
     lsame.f   chpr2.f  ctbsv.f  dspmv.f  dtbmv.f  dtpsv.f ssbmv.f  sspr.f   stpmv.f     zhpr2.f  ztbsv.f chbmv.f  chpr.f   ctpmv.f     dspr2.f dtbsv.f  sspmv.f  stbmv.f  stpsv.f   zhbmv.f  zhpr.f   ztpmv.f chpmv.f  ctbmv.f  ctpsv.f    dsbmv.f  dspr.f   dtpmv.f   sspr2.f  stbsv.f  zhpmv.f  ztbmv.f  ztpsv.f
 )
@@ -41,7 +42,7 @@ install(TARGETS eigen_blas
 if(EIGEN_LEAVE_TEST_IN_ALL_TARGET)
   add_subdirectory(testing) # can't do EXCLUDE_FROM_ALL here, breaks CTest
 else()
-  add_subdirectory(testing)
+  add_subdirectory(testing EXCLUDE_FROM_ALL)
 endif()
 
 endif(CMAKE_Fortran_COMPILER_WORKS)

blas/common.h

@@ -105,6 +105,7 @@ enum
   Conj = IsComplex
 };
 
+typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> PlainMatrixType;
 typedef Map<Matrix<Scalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > MatrixType;
 typedef Map<Matrix<Scalar,Dynamic,1>, 0, InnerStride<Dynamic> > StridedVectorType;
 typedef Map<Matrix<Scalar,Dynamic,1> > CompactVectorType;

blas/complexdots.f

@@ -0,0 +1,43 @@
+      COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      COMPLEX CX(*),CY(*)
+      COMPLEX RES
+      EXTERNAL CDOTCW
+      
+      CALL CDOTCW(N,CX,INCX,CY,INCY,RES)
+      CDOTC = RES
+      RETURN
+      END
+      
+      COMPLEX FUNCTION CDOTU(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      COMPLEX CX(*),CY(*)
+      COMPLEX RES
+      EXTERNAL CDOTUW
+      
+      CALL CDOTUW(N,CX,INCX,CY,INCY,RES)
+      CDOTU = RES
+      RETURN
+      END
+      
+      DOUBLE COMPLEX FUNCTION ZDOTC(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      DOUBLE COMPLEX CX(*),CY(*)
+      DOUBLE COMPLEX RES
+      EXTERNAL ZDOTCW
+      
+      CALL ZDOTCW(N,CX,INCX,CY,INCY,RES)
+      ZDOTC = RES
+      RETURN
+      END
+      
+      DOUBLE COMPLEX FUNCTION ZDOTU(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      DOUBLE COMPLEX CX(*),CY(*)
+      DOUBLE COMPLEX RES
+      EXTERNAL ZDOTUW
+      
+      CALL ZDOTUW(N,CX,INCX,CY,INCY,RES)
+      ZDOTU = RES
+      RETURN
+      END

blas/level1_cplx_impl.h

@@ -52,7 +52,7 @@ RealScalar EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),asum_)(int *n,
 }
 
 // computes a dot product of a conjugated vector with another vector.
-Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+int EIGEN_BLAS_FUNC(dotcw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres)
 {
 //   std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n";
 
@@ -60,18 +60,18 @@ Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py,
 
   Scalar* x = reinterpret_cast<Scalar*>(px);
   Scalar* y = reinterpret_cast<Scalar*>(py);
+  Scalar* res = reinterpret_cast<Scalar*>(pres);
 
-  Scalar res;
-  if(*incx==1 && *incy==1)    res = (vector(x,*n).dot(vector(y,*n)));
-  else if(*incx>0 && *incy>0) res = (vector(x,*n,*incx).dot(vector(y,*n,*incy)));
-  else if(*incx<0 && *incy>0) res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,*incy)));
-  else if(*incx>0 && *incy<0) res = (vector(x,*n,*incx).dot(vector(y,*n,-*incy).reverse()));
-  else if(*incx<0 && *incy<0) res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,-*incy).reverse()));
-  return res;
+  if(*incx==1 && *incy==1)    *res = (vector(x,*n).dot(vector(y,*n)));
+  else if(*incx>0 && *incy>0) *res = (vector(x,*n,*incx).dot(vector(y,*n,*incy)));
+  else if(*incx<0 && *incy>0) *res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,*incy)));
+  else if(*incx>0 && *incy<0) *res = (vector(x,*n,*incx).dot(vector(y,*n,-*incy).reverse()));
+  else if(*incx<0 && *incy<0) *res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,-*incy).reverse()));
+  return 0;
 }
 
 // computes a vector-vector dot product without complex conjugation.
-Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+int EIGEN_BLAS_FUNC(dotuw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres)
 {
 //   std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n";
 
@@ -79,13 +79,14 @@ Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py,
 
   Scalar* x = reinterpret_cast<Scalar*>(px);
   Scalar* y = reinterpret_cast<Scalar*>(py);
-  Scalar res;
-  if(*incx==1 && *incy==1)    res = (vector(x,*n).cwiseProduct(vector(y,*n))).sum();
-  else if(*incx>0 && *incy>0) res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum();
-  else if(*incx<0 && *incy>0) res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum();
-  else if(*incx>0 && *incy<0) res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
-  else if(*incx<0 && *incy<0) res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
-  return res;
+  Scalar* res = reinterpret_cast<Scalar*>(pres);
+
+  if(*incx==1 && *incy==1)    *res = (vector(x,*n).cwiseProduct(vector(y,*n))).sum();
+  else if(*incx>0 && *incy>0) *res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum();
+  else if(*incx<0 && *incy>0) *res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum();
+  else if(*incx>0 && *incy<0) *res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
+  else if(*incx<0 && *incy<0) *res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
+  return 0;
 }
 
 RealScalar EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),nrm2_)(int *n, RealScalar *px, int *incx)

blas/level3_impl.h

@@ -343,8 +343,6 @@ int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palp
   if(info)
     return xerbla_(SCALAR_SUFFIX_UP"SYRK ",&info,6);
 
-  int code = OP(*op) | (UPLO(*uplo) << 2);
-
   if(beta!=Scalar(1))
   {
     if(UPLO(*uplo)==UP)
@@ -372,6 +370,7 @@ int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palp
       matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
   }
   #else
+  int code = OP(*op) | (UPLO(*uplo) << 2);
   func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha);
   #endif
 

blas/testing/CMakeLists.txt

@@ -15,6 +15,7 @@ macro(ei_add_blas_test testname)
   target_link_libraries(${targetname} ${EXTERNAL_LIBS})
 
   add_test(${testname} "${Eigen_SOURCE_DIR}/blas/testing/runblastest.sh" "${testname}" "${Eigen_SOURCE_DIR}/blas/testing/${testname}.dat")
+  add_dependencies(buildtests ${targetname})
   
 endmacro(ei_add_blas_test)
 

cmake/FindCholmod.cmake

@@ -1,3 +1,5 @@
+# Cholmod lib usually requires linking to a blas and lapack library.
+# It is up to the user of this module to find a BLAS and link to it.
 
 if (CHOLMOD_INCLUDES AND CHOLMOD_LIBRARIES)
   set(CHOLMOD_FIND_QUIETLY TRUE)