bump to 3.1.0-rc2

function which must be called at the initialization time of any multi-threaded
application calling Eigen from multiple threads.

CMakeLists.txt

@@ -338,12 +338,9 @@ if(EIGEN_BUILD_BTL)
   add_subdirectory(bench/btl EXCLUDE_FROM_ALL)
 endif(EIGEN_BUILD_BTL)
 
-if(TEST_REAL_CASES)
-  if(NOT WIN32)
-    add_subdirectory(bench/spbench EXCLUDE_FROM_ALL)
-    set(ENV(EIGEN_MATRIX_DIR) ${TEST_REAL_CASES})
-  endif(NOT WIN32)
-endif(TEST_REAL_CASES)
+if(NOT WIN32)
+  add_subdirectory(bench/spbench EXCLUDE_FROM_ALL)
+endif(NOT WIN32)
 
 ei_testing_print_summary()
 

Eigen/CholmodSupport

@@ -26,6 +26,11 @@ extern "C" {
   * \code
   * #include <Eigen/CholmodSupport>
   * \endcode
+  *
+  * In order to use this module, the cholmod headers must be accessible from the include paths, and your binary must be linked to the cholmod library and its dependencies.
+  * The dependencies depend on how cholmod has been compiled.
+  * For a cmake based project, you can use our FindCholmod.cmake module to help you in this task.
+  *
   */
 
 #include "src/misc/Solve.h"

Eigen/Core

@@ -329,12 +329,12 @@ using std::ptrdiff_t;
 #include "src/Core/GeneralProduct.h"
 #include "src/Core/TriangularMatrix.h"
 #include "src/Core/SelfAdjointView.h"
-#include "src/Core/SolveTriangular.h"
+#include "src/Core/products/GeneralBlockPanelKernel.h"
 #include "src/Core/products/Parallelizer.h"
 #include "src/Core/products/CoeffBasedProduct.h"
-#include "src/Core/products/GeneralBlockPanelKernel.h"
 #include "src/Core/products/GeneralMatrixVector.h"
 #include "src/Core/products/GeneralMatrixMatrix.h"
+#include "src/Core/SolveTriangular.h"
 #include "src/Core/products/GeneralMatrixMatrixTriangular.h"
 #include "src/Core/products/SelfadjointMatrixVector.h"
 #include "src/Core/products/SelfadjointMatrixMatrix.h"

Eigen/PaStiXSupport

@@ -17,10 +17,22 @@ extern "C" {
 
 /** \ingroup Support_modules
   * \defgroup PaStiXSupport_Module PaStiXSupport module
-  *
+  * 
+  * This module provides an interface to the <a href="http://pastix.gforge.inria.fr/">PaSTiX</a> library.
+  * PaSTiX is a general \b supernodal, \b parallel and \b opensource sparse solver.
+  * It provides the two following main factorization classes:
+  * - class PastixLLT : a supernodal, parallel LLt Cholesky factorization.
+  * - class PastixLDLT: a supernodal, parallel LDLt Cholesky factorization.
+  * - class PastixLU : a supernodal, parallel LU factorization (optimized for a symmetric pattern).
+  * 
   * \code
   * #include <Eigen/PaStiXSupport>
   * \endcode
+  *
+  * In order to use this module, the PaSTiX headers must be accessible from the include paths, and your binary must be linked to the PaSTiX library and its dependencies.
+  * The dependencies depend on how PaSTiX has been compiled.
+  * For a cmake based project, you can use our FindPaSTiX.cmake module to help you in this task.
+  *
   */
 
 #include "src/misc/Solve.h"

Eigen/PardisoSupport

@@ -12,11 +12,15 @@
 /** \ingroup Support_modules
   * \defgroup PardisoSupport_Module PardisoSupport module
   *
-  * This module brings support for the Intel(R) MKL PARDISO direct sparse solvers
+  * This module brings support for the Intel(R) MKL PARDISO direct sparse solvers.
   *
   * \code
   * #include <Eigen/PardisoSupport>
   * \endcode
+  *
+  * In order to use this module, the MKL headers must be accessible from the include paths, and your binary must be linked to the MKL library and its dependencies.
+  * See this \ref TopicUsingIntelMKL "page" for more information on MKL-Eigen integration.
+  * 
   */
 
 #include "src/PardisoSupport/PardisoSupport.h"

Eigen/SuperLUSupport

@@ -31,11 +31,21 @@ namespace Eigen { struct SluMatrix; }
 /** \ingroup Support_modules
   * \defgroup SuperLUSupport_Module SuperLUSupport module
   *
+  * This module provides an interface to the <a href="http://crd-legacy.lbl.gov/~xiaoye/SuperLU/">SuperLU</a> library.
+  * It provides the following factorization class:
+  * - class SuperLU: a supernodal sequential LU factorization.
+  * - class SuperILU: a supernodal sequential incomplete LU factorization (to be used as a preconditioner for iterative methods).
+  *
   * \warning When including this module, you have to use SUPERLU_EMPTY instead of EMPTY which is no longer defined because it is too polluting.
   *
   * \code
   * #include <Eigen/SuperLUSupport>
   * \endcode
+  *
+  * In order to use this module, the superlu headers must be accessible from the include paths, and your binary must be linked to the superlu library and its dependencies.
+  * The dependencies depend on how superlu has been compiled.
+  * For a cmake based project, you can use our FindSuperLU.cmake module to help you in this task.
+  *
   */
 
 #include "src/misc/Solve.h"

Eigen/UmfPackSupport

@@ -12,12 +12,18 @@ extern "C" {
 /** \ingroup Support_modules
   * \defgroup UmfPackSupport_Module UmfPackSupport module
   *
-  *
-  *
+  * This module provides an interface to the UmfPack library which is part of the <a href="http://www.cise.ufl.edu/research/sparse/SuiteSparse/">suitesparse</a> package.
+  * It provides the following factorization class:
+  * - class UmfPackLU: a multifrontal sequential LU factorization.
   *
   * \code
   * #include <Eigen/UmfPackSupport>
   * \endcode
+  *
+  * In order to use this module, the umfpack headers must be accessible from the include paths, and your binary must be linked to the umfpack library and its dependencies.
+  * The dependencies depend on how umfpack has been compiled.
+  * For a cmake based project, you can use our FindUmfPack.cmake module to help you in this task.
+  *
   */
 
 #include "src/misc/Solve.h"

Eigen/src/Cholesky/LDLT.h

@@ -479,7 +479,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Deri
     for (Index i = 0; i < size; i++)
       m_transpositions.coeffRef(i) = i;
     m_temporary.resize(size);
-    m_sign = sigma;
+    m_sign = sigma>=0 ? 1 : -1;
     m_isInitialized = true;
   }
 

Eigen/src/Core/Assign_MKL.h

@@ -209,10 +209,13 @@ EIGEN_MKL_VML_DECLARE_UNARY_CALLS_LA(sqrt, Sqrt)
 
 EIGEN_MKL_VML_DECLARE_UNARY_CALLS_REAL(square, Sqr)
 
+// The vm*powx functions are not avaibale in the windows version of MKL.
+#ifdef _WIN32
 EIGEN_MKL_VML_DECLARE_POW_CALL(pow, vmspowx_, float, float)
 EIGEN_MKL_VML_DECLARE_POW_CALL(pow, vmdpowx_, double, double)
 EIGEN_MKL_VML_DECLARE_POW_CALL(pow, vmcpowx_, scomplex, MKL_Complex8)
 EIGEN_MKL_VML_DECLARE_POW_CALL(pow, vmzpowx_, dcomplex, MKL_Complex16)
+#endif
 
 } // end namespace internal
 

Eigen/src/Core/Functors.h

@@ -295,7 +295,7 @@ struct functor_traits<scalar_opposite_op<Scalar> >
 template<typename Scalar> struct scalar_abs_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_abs_op)
   typedef typename NumTraits<Scalar>::Real result_type;
-  EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return abs(a); }
+  EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return internal::abs(a); }
   template<typename Packet>
   EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
   { return internal::pabs(a); }
@@ -317,7 +317,7 @@ struct functor_traits<scalar_abs_op<Scalar> >
 template<typename Scalar> struct scalar_abs2_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_abs2_op)
   typedef typename NumTraits<Scalar>::Real result_type;
-  EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return abs2(a); }
+  EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return internal::abs2(a); }
   template<typename Packet>
   EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
   { return internal::pmul(a,a); }
@@ -333,7 +333,7 @@ struct functor_traits<scalar_abs2_op<Scalar> >
   */
 template<typename Scalar> struct scalar_conjugate_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_conjugate_op)
-  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { return conj(a); }
+  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { return internal::conj(a); }
   template<typename Packet>
   EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const { return internal::pconj(a); }
 };
@@ -370,7 +370,7 @@ template<typename Scalar>
 struct scalar_real_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_real_op)
   typedef typename NumTraits<Scalar>::Real result_type;
-  EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return real(a); }
+  EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return internal::real(a); }
 };
 template<typename Scalar>
 struct functor_traits<scalar_real_op<Scalar> >
@@ -385,7 +385,7 @@ template<typename Scalar>
 struct scalar_imag_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_imag_op)
   typedef typename NumTraits<Scalar>::Real result_type;
-  EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return imag(a); }
+  EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return internal::imag(a); }
 };
 template<typename Scalar>
 struct functor_traits<scalar_imag_op<Scalar> >
@@ -400,7 +400,7 @@ template<typename Scalar>
 struct scalar_real_ref_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_real_ref_op)
   typedef typename NumTraits<Scalar>::Real result_type;
-  EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return real_ref(*const_cast<Scalar*>(&a)); }
+  EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return internal::real_ref(*const_cast<Scalar*>(&a)); }
 };
 template<typename Scalar>
 struct functor_traits<scalar_real_ref_op<Scalar> >
@@ -415,7 +415,7 @@ template<typename Scalar>
 struct scalar_imag_ref_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_imag_ref_op)
   typedef typename NumTraits<Scalar>::Real result_type;
-  EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return imag_ref(*const_cast<Scalar*>(&a)); }
+  EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return internal::imag_ref(*const_cast<Scalar*>(&a)); }
 };
 template<typename Scalar>
 struct functor_traits<scalar_imag_ref_op<Scalar> >
@@ -429,7 +429,7 @@ struct functor_traits<scalar_imag_ref_op<Scalar> >
   */
 template<typename Scalar> struct scalar_exp_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_exp_op)
-  inline const Scalar operator() (const Scalar& a) const { return exp(a); }
+  inline const Scalar operator() (const Scalar& a) const { return internal::exp(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::pexp(a); }
 };
@@ -445,7 +445,7 @@ struct functor_traits<scalar_exp_op<Scalar> >
   */
 template<typename Scalar> struct scalar_log_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_log_op)
-  inline const Scalar operator() (const Scalar& a) const { return log(a); }
+  inline const Scalar operator() (const Scalar& a) const { return internal::log(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::plog(a); }
 };
@@ -703,7 +703,7 @@ struct functor_traits<scalar_add_op<Scalar> >
   */
 template<typename Scalar> struct scalar_sqrt_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_sqrt_op)
-  inline const Scalar operator() (const Scalar& a) const { return sqrt(a); }
+  inline const Scalar operator() (const Scalar& a) const { return internal::sqrt(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::psqrt(a); }
 };
@@ -721,7 +721,7 @@ struct functor_traits<scalar_sqrt_op<Scalar> >
   */
 template<typename Scalar> struct scalar_cos_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_cos_op)
-  inline Scalar operator() (const Scalar& a) const { return cos(a); }
+  inline Scalar operator() (const Scalar& a) const { return internal::cos(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::pcos(a); }
 };
@@ -740,7 +740,7 @@ struct functor_traits<scalar_cos_op<Scalar> >
   */
 template<typename Scalar> struct scalar_sin_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_sin_op)
-  inline const Scalar operator() (const Scalar& a) const { return sin(a); }
+  inline const Scalar operator() (const Scalar& a) const { return internal::sin(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::psin(a); }
 };
@@ -760,7 +760,7 @@ struct functor_traits<scalar_sin_op<Scalar> >
   */
 template<typename Scalar> struct scalar_tan_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_tan_op)
-  inline const Scalar operator() (const Scalar& a) const { return tan(a); }
+  inline const Scalar operator() (const Scalar& a) const { return internal::tan(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::ptan(a); }
 };
@@ -779,7 +779,7 @@ struct functor_traits<scalar_tan_op<Scalar> >
   */
 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); }
+  inline const Scalar operator() (const Scalar& a) const { return internal::acos(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::pacos(a); }
 };
@@ -798,7 +798,7 @@ struct functor_traits<scalar_acos_op<Scalar> >
   */
 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); }
+  inline const Scalar operator() (const Scalar& a) const { return internal::asin(a); }
   typedef typename packet_traits<Scalar>::type Packet;
   inline Packet packetOp(const Packet& a) const { return internal::pasin(a); }
 };

Eigen/src/Core/SolveTriangular.h

@@ -100,12 +100,22 @@ struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,Dynamic>
   typedef typename Rhs::Index Index;
   typedef blas_traits<Lhs> LhsProductTraits;
   typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType;
+
   static void run(const Lhs& lhs, Rhs& rhs)
   {
     typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsProductTraits::extract(lhs);
+
+    const Index size = lhs.rows();
+    const Index othersize = Side==OnTheLeft? rhs.cols() : rhs.rows();
+
+    typedef internal::gemm_blocking_space<(Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor,Scalar,Scalar,
+              Rhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxRowsAtCompileTime,4> BlockingType;
+
+    BlockingType blocking(rhs.rows(), rhs.cols(), size);
+
     triangular_solve_matrix<Scalar,Index,Side,Mode,LhsProductTraits::NeedToConjugate,(int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor,
                                (Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor>
-      ::run(lhs.rows(), Side==OnTheLeft? rhs.cols() : rhs.rows(), &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.outerStride());
+      ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.outerStride(), blocking);
   }
 };
 

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

@@ -123,7 +123,7 @@ Packet4f pexp<Packet4f>(const Packet4f& _x)
   _EIGEN_DECLARE_CONST_Packet4i(0x7f, 0x7f);
 
 
-  _EIGEN_DECLARE_CONST_Packet4f(exp_hi, 88.3762626647949f);
+  _EIGEN_DECLARE_CONST_Packet4f(exp_hi,  88.3762626647950f);
   _EIGEN_DECLARE_CONST_Packet4f(exp_lo, -88.3762626647949f);
 
   _EIGEN_DECLARE_CONST_Packet4f(cephes_LOG2EF, 1.44269504088896341f);
@@ -170,7 +170,7 @@ Packet4f pexp<Packet4f>(const Packet4f& _x)
   y = pmadd(y, z, x);
   y = padd(y, p4f_1);
 
-  /* build 2^n */
+  // build 2^n
   emm0 = _mm_cvttps_epi32(fx);
   emm0 = _mm_add_epi32(emm0, p4i_0x7f);
   emm0 = _mm_slli_epi32(emm0, 23);

Eigen/src/Core/products/GeneralBlockPanelKernel.h

@@ -26,7 +26,7 @@
 #define EIGEN_GENERAL_BLOCK_PANEL_H
 
 namespace Eigen { 
-
+  
 namespace internal {
 
 template<typename _LhsScalar, typename _RhsScalar, bool _ConjLhs=false, bool _ConjRhs=false>
@@ -42,9 +42,14 @@ inline std::ptrdiff_t manage_caching_sizes_helper(std::ptrdiff_t a, std::ptrdiff
 /** \internal */
 inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1=0, std::ptrdiff_t* l2=0)
 {
-  static std::ptrdiff_t m_l1CacheSize = manage_caching_sizes_helper(queryL1CacheSize(),8 * 1024);
-  static std::ptrdiff_t m_l2CacheSize = manage_caching_sizes_helper(queryTopLevelCacheSize(),1*1024*1024);
-
+  static std::ptrdiff_t m_l1CacheSize = 0;
+  static std::ptrdiff_t m_l2CacheSize = 0;
+  if(m_l2CacheSize==0)
+  {
+    m_l1CacheSize = manage_caching_sizes_helper(queryL1CacheSize(),8 * 1024);
+    m_l2CacheSize = manage_caching_sizes_helper(queryTopLevelCacheSize(),1*1024*1024);
+  }
+  
   if(action==SetAction)
   {
     // set the cpu cache size and cache all block sizes from a global cache size in byte

Eigen/src/Core/products/GeneralMatrixMatrix.h

@@ -79,7 +79,7 @@ static void run(Index rows, Index cols, Index depth,
 
   typedef gebp_traits<LhsScalar,RhsScalar> Traits;
 
-  Index kc = blocking.kc();                 // cache block size along the K direction
+  Index kc = blocking.kc();                   // cache block size along the K direction
   Index mc = (std::min)(rows,blocking.mc());  // cache block size along the M direction
   //Index nc = blocking.nc(); // cache block size along the N direction
 
@@ -249,7 +249,7 @@ struct gemm_functor
     BlockingType& m_blocking;
 };
 
-template<int StorageOrder, typename LhsScalar, typename RhsScalar, int MaxRows, int MaxCols, int MaxDepth,
+template<int StorageOrder, typename LhsScalar, typename RhsScalar, int MaxRows, int MaxCols, int MaxDepth, int KcFactor=1,
 bool FiniteAtCompileTime = MaxRows!=Dynamic && MaxCols!=Dynamic && MaxDepth != Dynamic> class gemm_blocking_space;
 
 template<typename _LhsScalar, typename _RhsScalar>
@@ -282,8 +282,8 @@ class level3_blocking
     inline RhsScalar* blockW() { return m_blockW; }
 };
 
-template<int StorageOrder, typename _LhsScalar, typename _RhsScalar, int MaxRows, int MaxCols, int MaxDepth>
-class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, true>
+template<int StorageOrder, typename _LhsScalar, typename _RhsScalar, int MaxRows, int MaxCols, int MaxDepth, int KcFactor>
+class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, KcFactor, true>
   : public level3_blocking<
       typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
       typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
@@ -324,8 +324,8 @@ class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, M
     inline void allocateAll() {}
 };
 
-template<int StorageOrder, typename _LhsScalar, typename _RhsScalar, int MaxRows, int MaxCols, int MaxDepth>
-class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, false>
+template<int StorageOrder, typename _LhsScalar, typename _RhsScalar, int MaxRows, int MaxCols, int MaxDepth, int KcFactor>
+class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, KcFactor, false>
   : public level3_blocking<
       typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
       typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
@@ -349,7 +349,7 @@ class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, M
       this->m_nc = Transpose ? rows : cols;
       this->m_kc = depth;
 
-      computeProductBlockingSizes<LhsScalar,RhsScalar>(this->m_kc, this->m_mc, this->m_nc);
+      computeProductBlockingSizes<LhsScalar,RhsScalar,KcFactor>(this->m_kc, this->m_mc, this->m_nc);
       m_sizeA = this->m_mc * this->m_kc;
       m_sizeB = this->m_kc * this->m_nc;
       m_sizeW = this->m_kc*Traits::WorkSpaceFactor;

Eigen/src/Core/products/Parallelizer.h

@@ -57,12 +57,23 @@ inline void manage_multi_threading(Action action, int* v)
   }
 }
 
+}
+
+/** Must be call first when calling Eigen from multiple threads */
+inline void initParallel()
+{
+  int nbt;
+  internal::manage_multi_threading(GetAction, &nbt);
+  std::ptrdiff_t l1, l2;
+  internal::manage_caching_sizes(GetAction, &l1, &l2);
+}
+
 /** \returns the max number of threads reserved for Eigen
   * \sa setNbThreads */
 inline int nbThreads()
 {
   int ret;
-  manage_multi_threading(GetAction, &ret);
+  internal::manage_multi_threading(GetAction, &ret);
   return ret;
 }
 
@@ -70,9 +81,11 @@ inline int nbThreads()
   * \sa nbThreads */
 inline void setNbThreads(int v)
 {
-  manage_multi_threading(SetAction, &v);
+  internal::manage_multi_threading(SetAction, &v);
 }
 
+namespace internal {
+
 template<typename Index> struct GemmParallelInfo
 {
   GemmParallelInfo() : sync(-1), users(0), rhs_start(0), rhs_length(0) {}
@@ -121,6 +134,7 @@ void parallelize_gemm(const Functor& func, Index rows, Index cols, bool transpos
   if(threads==1)
     return func(0,rows, 0,cols);
 
+  Eigen::initParallel();
   func.initParallelSession();
 
   if(transpose)

Eigen/src/Core/products/TriangularSolverMatrix.h

@@ -36,14 +36,15 @@ struct triangular_solve_matrix<Scalar,Index,Side,Mode,Conjugate,TriStorageOrder,
   static EIGEN_DONT_INLINE void run(
     Index size, Index cols,
     const Scalar*  tri, Index triStride,
-    Scalar* _other, Index otherStride)
+    Scalar* _other, Index otherStride,
+    level3_blocking<Scalar,Scalar>& blocking)
   {
     triangular_solve_matrix<
       Scalar, Index, Side==OnTheLeft?OnTheRight:OnTheLeft,
       (Mode&UnitDiag) | ((Mode&Upper) ? Lower : Upper),
       NumTraits<Scalar>::IsComplex && Conjugate,
       TriStorageOrder==RowMajor ? ColMajor : RowMajor, ColMajor>
-      ::run(size, cols, tri, triStride, _other, otherStride);
+      ::run(size, cols, tri, triStride, _other, otherStride, blocking);
   }
 };
 
@@ -55,7 +56,8 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
   static EIGEN_DONT_INLINE void run(
     Index size, Index otherSize,
     const Scalar* _tri, Index triStride,
-    Scalar* _other, Index otherStride)
+    Scalar* _other, Index otherStride,
+    level3_blocking<Scalar,Scalar>& blocking)
   {
     Index cols = otherSize;
     const_blas_data_mapper<Scalar, Index, TriStorageOrder> tri(_tri,triStride);
@@ -67,17 +69,16 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
       IsLower = (Mode&Lower) == Lower
     };
 
-    Index kc = size; // cache block size along the K direction
-    Index mc = size;  // cache block size along the M direction
-    Index nc = cols;  // cache block size along the N direction
-    computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
+    Index kc = blocking.kc();                   // cache block size along the K direction
+    Index mc = (std::min)(size,blocking.mc());  // cache block size along the M direction
 
+    std::size_t sizeA = kc*mc;
+    std::size_t sizeB = kc*cols;
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
-    std::size_t sizeB = sizeW + kc*cols;
-    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;
-    Scalar* blockW = allocatedBlockB;
+
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, sizeA, blocking.blockA());
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockB, sizeB, blocking.blockB());
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockW, sizeW, blocking.blockW());
 
     conj_if<Conjugate> conj;
     gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, Conjugate, false> gebp_kernel;
@@ -181,7 +182,7 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
           {
             pack_lhs(blockA, &tri(i2, IsLower ? k2 : k2-kc), triStride, actual_kc, actual_mc);
 
-            gebp_kernel(_other+i2, otherStride, blockA, blockB, actual_mc, actual_kc, cols, Scalar(-1));
+            gebp_kernel(_other+i2, otherStride, blockA, blockB, actual_mc, actual_kc, cols, Scalar(-1), -1, -1, 0, 0, blockW);
           }
         }
       }
@@ -197,7 +198,8 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
   static EIGEN_DONT_INLINE void run(
     Index size, Index otherSize,
     const Scalar* _tri, Index triStride,
-    Scalar* _other, Index otherStride)
+    Scalar* _other, Index otherStride,
+    level3_blocking<Scalar,Scalar>& blocking)
   {
     Index rows = otherSize;
     const_blas_data_mapper<Scalar, Index, TriStorageOrder> rhs(_tri,triStride);
@@ -210,19 +212,16 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
       IsLower = (Mode&Lower) == Lower
     };
 
-//     Index kc = std::min<Index>(Traits::Max_kc/4,size); // cache block size along the K direction
-//     Index mc = std::min<Index>(Traits::Max_mc,size);   // cache block size along the M direction
-    // check that !!!!
-    Index kc = size; // cache block size along the K direction
-    Index mc = size;  // cache block size along the M direction
-    Index nc = rows;  // cache block size along the N direction
-    computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
+    Index kc = blocking.kc();                   // cache block size along the K direction
+    Index mc = (std::min)(rows,blocking.mc());  // cache block size along the M direction
 
+    std::size_t sizeA = kc*mc;
+    std::size_t sizeB = kc*size;
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
-    std::size_t sizeB = sizeW + kc*size;
-    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;
+
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, sizeA, blocking.blockA());
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockB, sizeB, blocking.blockB());
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockW, sizeW, blocking.blockW());
 
     conj_if<Conjugate> conj;
     gebp_kernel<Scalar,Scalar, Index, Traits::mr, Traits::nr, false, Conjugate> gebp_kernel;
@@ -289,7 +288,7 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
                           Scalar(-1),
                           actual_kc, actual_kc, // strides
                           panelOffset, panelOffset, // offsets
-                          allocatedBlockB);  // workspace
+                          blockW);  // workspace
             }
 
             // unblocked triangular solve
@@ -320,7 +319,7 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
         if (rs>0)
           gebp_kernel(_other+i2+startPanel*otherStride, otherStride, blockA, geb,
                       actual_mc, actual_kc, rs, Scalar(-1),
-                      -1, -1, 0, 0, allocatedBlockB);
+                      -1, -1, 0, 0, blockW);
       }
     }
   }

Eigen/src/Core/products/TriangularSolverMatrix_MKL.h

@@ -51,7 +51,7 @@ struct triangular_solve_matrix<EIGTYPE,Index,OnTheLeft,Mode,Conjugate,TriStorage
   static EIGEN_DONT_INLINE void run( \
       Index size, Index otherSize, \
       const EIGTYPE* _tri, Index triStride, \
-      EIGTYPE* _other, Index otherStride) \
+      EIGTYPE* _other, Index otherStride, level3_blocking<EIGTYPE,EIGTYPE>& /*blocking*/) \
   { \
    MKL_INT m = size, n = otherSize, lda, ldb; \
    char side = 'L', uplo, diag='N', transa; \
@@ -106,7 +106,7 @@ struct triangular_solve_matrix<EIGTYPE,Index,OnTheRight,Mode,Conjugate,TriStorag
   static EIGEN_DONT_INLINE void run( \
       Index size, Index otherSize, \
       const EIGTYPE* _tri, Index triStride, \
-      EIGTYPE* _other, Index otherStride) \
+      EIGTYPE* _other, Index otherStride, level3_blocking<EIGTYPE,EIGTYPE>& /*blocking*/) \
   { \
    MKL_INT m = otherSize, n = size, lda, ldb; \
    char side = 'R', uplo, diag='N', transa; \

Eigen/src/Core/util/Macros.h

@@ -28,7 +28,7 @@
 
 #define EIGEN_WORLD_VERSION 3
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 93
+#define EIGEN_MINOR_VERSION 95
 
 #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/src/Eigenvalues/RealSchur.h

@@ -239,38 +239,40 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
   Scalar exshift(0); // sum of exceptional shifts
   Scalar norm = computeNormOfT();
 
-  while (iu >= 0)
+  if(norm!=0)
   {
-    Index il = findSmallSubdiagEntry(iu, norm);
+    while (iu >= 0)
+    {
+      Index il = findSmallSubdiagEntry(iu, norm);
 
-    // Check for convergence
-    if (il == iu) // One root found
-    {
-      m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
-      if (iu > 0) 
-        m_matT.coeffRef(iu, iu-1) = Scalar(0);
-      iu--;
-      iter = 0;
+      // Check for convergence
+      if (il == iu) // One root found
+      {
+        m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
+        if (iu > 0)
+          m_matT.coeffRef(iu, iu-1) = Scalar(0);
+        iu--;
+        iter = 0;
+      }
+      else if (il == iu-1) // Two roots found
+      {
+        splitOffTwoRows(iu, computeU, exshift);
+        iu -= 2;
+        iter = 0;
+      }
+      else // No convergence yet
+      {
+        // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
+        Vector3s firstHouseholderVector(0,0,0), shiftInfo;
+        computeShift(iu, iter, exshift, shiftInfo);
+        iter = iter + 1;
+        if (iter > m_maxIterations) break;
+        Index im;
+        initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
+        performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
+      }
     }
-    else if (il == iu-1) // Two roots found
-    {
-      splitOffTwoRows(iu, computeU, exshift);
-      iu -= 2;
-      iter = 0;
-    }
-    else // No convergence yet
-    {
-      // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
-      Vector3s firstHouseholderVector(0,0,0), shiftInfo;
-      computeShift(iu, iter, exshift, shiftInfo);
-      iter = iter + 1; 
-      if (iter > m_maxIterations) break;
-      Index im;
-      initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
-      performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
-    }
-  } 
-
+  }
   if(iter <= m_maxIterations) 
     m_info = Success;
   else

Eigen/src/Geometry/Transform.h

@@ -39,7 +39,7 @@ struct transform_traits
     Dim = Transform::Dim,
     HDim = Transform::HDim,
     Mode = Transform::Mode,
-    IsProjective = (Mode==Projective)
+    IsProjective = (int(Mode)==int(Projective))
   };
 };
 
@@ -468,15 +468,40 @@ public:
   {
     return internal::transform_transform_product_impl<Transform,Transform>::run(*this,other);
   }
-
+  
+  #ifdef __INTEL_COMPILER
+private:
+  // this intermediate structure permits to workaround a bug in ICC 11:
+  //   error: template instantiation resulted in unexpected function type of "Eigen::Transform<double, 3, 32, 0>
+  //             (const Eigen::Transform<double, 3, 2, 0> &) const"
+  //  (the meaning of a name may have changed since the template declaration -- the type of the template is:
+  // "Eigen::internal::transform_transform_product_impl<Eigen::Transform<double, 3, 32, 0>,
+  //     Eigen::Transform<double, 3, Mode, Options>, <expression>>::ResultType (const Eigen::Transform<double, 3, Mode, Options> &) const")
+  // 
+  template<int OtherMode,int OtherOptions> struct icc_11_workaround
+  {
+    typedef internal::transform_transform_product_impl<Transform,Transform<Scalar,Dim,OtherMode,OtherOptions> > ProductType;
+    typedef typename ProductType::ResultType ResultType;
+  };
+  
+public:
   /** Concatenates two different transformations */
   template<int OtherMode,int OtherOptions>
-  inline const typename internal::transform_transform_product_impl<
-    Transform,Transform<Scalar,Dim,OtherMode,OtherOptions> >::ResultType
+  inline typename icc_11_workaround<OtherMode,OtherOptions>::ResultType
+    operator * (const Transform<Scalar,Dim,OtherMode,OtherOptions>& other) const
+  {
+    typedef typename icc_11_workaround<OtherMode,OtherOptions>::ProductType ProductType;
+    return ProductType::run(*this,other);
+  }
+  #else
+  /** Concatenates two different transformations */
+  template<int OtherMode,int OtherOptions>
+  inline typename internal::transform_transform_product_impl<Transform,Transform<Scalar,Dim,OtherMode,OtherOptions> >::ResultType
     operator * (const Transform<Scalar,Dim,OtherMode,OtherOptions>& other) const
   {
     return internal::transform_transform_product_impl<Transform,Transform<Scalar,Dim,OtherMode,OtherOptions> >::run(*this,other);
   }
+  #endif
 
   /** \sa MatrixBase::setIdentity() */
   void setIdentity() { m_matrix.setIdentity(); }

Eigen/src/OrderingMethods/Amd.h

@@ -106,6 +106,7 @@ Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Ind
 template<typename Scalar, typename Index>
 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm)
 {
+  using std::sqrt;
   typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
   
   int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
@@ -114,7 +115,7 @@ void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, Permutation
   unsigned int h;
   
   Index n = C.cols();
-  dense = std::max<Index> (16, 10 * sqrt ((double) n));   /* find dense threshold */
+  dense = std::max<Index> (16, Index(10 * sqrt(double(n))));   /* find dense threshold */
   dense = std::min<Index> (n-2, dense);
   
   Index cnz = C.nonZeros();

Eigen/src/PaStiXSupport/PaStiXSupport.h

@@ -35,7 +35,6 @@ namespace Eigen {
   *
   * \sa TutorialSparseDirectSolvers
   */
-
 template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
 template<typename _MatrixType, int Options> class PastixLLT;
 template<typename _MatrixType, int Options> class PastixLDLT;
@@ -75,32 +74,34 @@ namespace internal
   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, float *vals, int *perm, int * invp, float *x, int nbrhs, int *iparm, double *dparm)
   {
     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
-    if (nbrhs == 0) x = NULL;
+    if (nbrhs == 0) {x = NULL; nbrhs=1;}
     s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
   }
   
   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, double *vals, int *perm, int * invp, double *x, int nbrhs, int *iparm, double *dparm)
   {
     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
-    if (nbrhs == 0) x = NULL;
+    if (nbrhs == 0) {x = NULL; nbrhs=1;}
     d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
   }
   
   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<float> *vals, int *perm, int * invp, std::complex<float> *x, int nbrhs, int *iparm, double *dparm)
   {
+    if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
+    if (nbrhs == 0) {x = NULL; nbrhs=1;}
     c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm); 
   }
   
   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm)
   {
     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
-    if (nbrhs == 0) x = NULL;
+    if (nbrhs == 0) {x = NULL; nbrhs=1;}
     z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm); 
   }
 
   // Convert the matrix  to Fortran-style Numbering
   template <typename MatrixType>
-  void EigenToFortranNumbering (MatrixType& mat)
+  void c_to_fortran_numbering (MatrixType& mat)
   {
     if ( !(mat.outerIndexPtr()[0]) ) 
     { 
@@ -114,7 +115,7 @@ namespace internal
   
   // Convert to C-style Numbering
   template <typename MatrixType>
-  void EigenToCNumbering (MatrixType& mat)
+  void fortran_to_c_numbering (MatrixType& mat)
   {
     // Check the Numbering
     if ( mat.outerIndexPtr()[0] == 1 ) 
@@ -126,38 +127,12 @@ namespace internal
         --mat.innerIndexPtr()[i];
     }
   }
-  
-  // Symmetrize the graph of the input matrix 
-  // In : The Input matrix to symmetrize the pattern
-  // Out : The output matrix
-  // StrMatTrans : The structural pattern of the transpose of In; It is 
-  // used to optimize the future symmetrization with the same matrix pattern
-  // WARNING It is assumed here that successive calls to this routine are done 
-  // with matrices having the same pattern.
-  template <typename MatrixType>
-  void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans, bool& hasTranspose)
-  {
-    eigen_assert(In.cols()==In.rows() && " Can only symmetrize the graph of a square matrix");
-    if (!hasTranspose)
-    { //First call to this routine, need to compute the structural pattern of In^T
-      StrMatTrans = In.transpose();
-      // Set the elements of the matrix to zero 
-      for (int i = 0; i < StrMatTrans.rows(); i++) 
-      {
-        for (typename MatrixType::InnerIterator it(StrMatTrans, i); it; ++it)
-          it.valueRef() = 0.0;
-      }
-      hasTranspose = true;
-    }
-    Out = (StrMatTrans + In).eval(); 
-  }
-
 }
 
 // This is the base class to interface with PaStiX functions. 
 // Users should not used this class directly. 
 template <class Derived>
-class PastixBase
+class PastixBase : internal::noncopyable
 {
   public:
     typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
@@ -166,29 +141,19 @@ class PastixBase
     typedef typename MatrixType::RealScalar RealScalar;
     typedef typename MatrixType::Index Index;
     typedef Matrix<Scalar,Dynamic,1> Vector;
+    typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix;
     
   public:
     
-    PastixBase():m_initisOk(false),m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) 
+    PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0)
     {
-      m_pastixdata = 0;
-      m_hasTranspose = false;
-      PastixInit();
+      init();
     }
     
     ~PastixBase() 
     {
-       PastixDestroy();
+      clean();
     }
-    
-    // Initialize the Pastix data structure, check the matrix
-    void PastixInit(); 
-    
-    // Compute the ordering and the symbolic factorization
-    Derived& analyzePattern (MatrixType& mat);
-    
-    // Compute the numerical factorization
-    Derived& factorize (MatrixType& mat);
 
     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
       *
@@ -269,7 +234,6 @@ class PastixBase
     /** Return a reference to a particular index parameter of the DPARM vector 
      * \sa dparm()
      */
-    
     double& dparm(int idxparam)
     {
       return m_dparm(idxparam);
@@ -307,17 +271,27 @@ class PastixBase
     }
     
   protected:
+
+    // Initialize the Pastix data structure, check the matrix
+    void init(); 
+    
+    // Compute the ordering and the symbolic factorization
+    void analyzePattern(ColSpMatrix& mat);
+    
+    // Compute the numerical factorization
+    void factorize(ColSpMatrix& mat);
+    
     // Free all the data allocated by Pastix
-    void PastixDestroy()
+    void clean()
     {
       eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 
       m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
       m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
-      internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
-               m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
+      internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
+                             m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
     }
     
-    Derived& compute (MatrixType& mat);
+    void compute(ColSpMatrix& mat);
     
     int m_initisOk; 
     int m_analysisIsOk;
@@ -325,22 +299,12 @@ class PastixBase
     bool m_isInitialized;
     mutable ComputationInfo m_info; 
     mutable pastix_data_t *m_pastixdata; // Data structure for pastix
-    mutable SparseMatrix<Scalar, ColMajor> m_mat_null; // An input  null matrix
-    mutable Matrix<Scalar, Dynamic,1> m_vec_null; // An input null vector
-    mutable SparseMatrix<Scalar, ColMajor> m_StrMatTrans; // The transpose pattern of the input matrix
-    mutable bool m_hasTranspose; // The transpose of the current matrix has already been computed
     mutable int m_comm; // The MPI communicator identifier
-    mutable Matrix<Index,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
+    mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
     mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
     mutable Matrix<Index,Dynamic,1> m_perm;  // Permutation vector
     mutable Matrix<Index,Dynamic,1> m_invp;  // Inverse permutation vector
-    mutable int m_ordering; // ordering method to use
-    mutable int m_amalgamation; // level of amalgamation
     mutable int m_size; // Size of the matrix 
-    
-  private:
-    PastixBase(PastixBase& ) {}
-    
 }; 
 
  /** Initialize the PaStiX data structure. 
@@ -348,29 +312,29 @@ class PastixBase
    * \sa iparm() dparm()
    */
 template <class Derived>
-void PastixBase<Derived>::PastixInit()
+void PastixBase<Derived>::init()
 {
   m_size = 0; 
-  m_iparm.resize(IPARM_SIZE);
-  m_dparm.resize(DPARM_SIZE);
+  m_iparm.setZero(IPARM_SIZE);
+  m_dparm.setZero(DPARM_SIZE);
   
   m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
-  if(m_pastixdata)
-  { // This trick is used to reset the Pastix internal data between successive
-    // calls with (structural) different matrices
-    PastixDestroy();
-    m_pastixdata = 0;
-    m_iparm(IPARM_MODIFY_PARAMETER) = API_YES;
-    m_hasTranspose = false; 
-  }
+  pastix(&m_pastixdata, MPI_COMM_WORLD,
+         0, 0, 0, 0,
+         0, 0, 0, 1, m_iparm.data(), m_dparm.data());
+  
+  m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
+  m_iparm[IPARM_VERBOSE]             = 2;
+  m_iparm[IPARM_ORDERING]            = API_ORDER_SCOTCH;
+  m_iparm[IPARM_INCOMPLETE]          = API_NO;
+  m_iparm[IPARM_OOC_LIMIT]           = 2000;
+  m_iparm[IPARM_RHS_MAKING]          = API_RHS_B;
+  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
   
   m_iparm(IPARM_START_TASK) = API_TASK_INIT;
   m_iparm(IPARM_END_TASK) = API_TASK_INIT;
-  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
-  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
-               m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
-  
-  m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
+  internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
+                         0, 0, 0, 0, m_iparm.data(), m_dparm.data());
   
   // Check the returned error
   if(m_iparm(IPARM_ERROR_NUMBER)) {
@@ -384,82 +348,74 @@ void PastixBase<Derived>::PastixInit()
 }
 
 template <class Derived>
-Derived& PastixBase<Derived>::compute(MatrixType& mat)
+void PastixBase<Derived>::compute(ColSpMatrix& mat)
 {
   eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
-  typedef typename MatrixType::Scalar Scalar;
   
-  // Save the size of the current matrix 
-  m_size = mat.rows();
-  // Convert the matrix in fortran-style numbering 
-  internal::EigenToFortranNumbering(mat);
-  analyzePattern(mat);
+  analyzePattern(mat);  
   factorize(mat);
+  
   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
-  if (m_factorizationIsOk) m_isInitialized = true;
-  
-  //Convert back the matrix -- Is it really necessary here 
-  internal::EigenToCNumbering(mat);
-  
-  return derived();
+  m_isInitialized = m_factorizationIsOk;
 }
 
 
 template <class Derived>
-Derived& PastixBase<Derived>::analyzePattern(MatrixType& mat)
-{
-  eigen_assert(m_initisOk && "PastixInit should be called first to set the default  parameters");
+void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
+{                         
+  eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
+  
+  // clean previous calls
+  if(m_size>0)
+    clean();
+  
   m_size = mat.rows();
   m_perm.resize(m_size);
   m_invp.resize(m_size);
   
-  // Convert the matrix in fortran-style numbering 
-  internal::EigenToFortranNumbering(mat);
-  
   m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
   m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
-  
   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
-               mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
+               mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
   
   // Check the returned error
-  if(m_iparm(IPARM_ERROR_NUMBER)) {
+  if(m_iparm(IPARM_ERROR_NUMBER))
+  {
     m_info = NumericalIssue;
     m_analysisIsOk = false;
   }
-  else { 
+  else
+  { 
     m_info = Success;
     m_analysisIsOk = true;
   }
-  return derived();
 }
 
 template <class Derived>
-Derived& PastixBase<Derived>::factorize(MatrixType& mat)
+void PastixBase<Derived>::factorize(ColSpMatrix& mat)
 {
+//   if(&m_cpyMat != &mat) m_cpyMat = mat;
   eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
   m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
   m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
   m_size = mat.rows();
   
-  // Convert the matrix in fortran-style numbering 
-  internal::EigenToFortranNumbering(mat);
-  
   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
-               mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
+               mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
   
   // Check the returned error
-  if(m_iparm(IPARM_ERROR_NUMBER)) {
+  if(m_iparm(IPARM_ERROR_NUMBER))
+  {
     m_info = NumericalIssue;
     m_factorizationIsOk = false;
     m_isInitialized = false;
   }
-  else { 
+  else
+  {
     m_info = Success;
     m_factorizationIsOk = true;
     m_isInitialized = true;
   }
-  return derived();
 }
 
 /* Solve the system */
@@ -475,20 +431,17 @@ bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) co
   x = b; /* on return, x is overwritten by the computed solution */
   
   for (int i = 0; i < b.cols(); i++){
-    m_iparm(IPARM_START_TASK) = API_TASK_SOLVE;
-    m_iparm(IPARM_END_TASK) = API_TASK_REFINE;
-    m_iparm(IPARM_RHS_MAKING) = API_RHS_B;
-    internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
-               m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
+    m_iparm[IPARM_START_TASK]          = API_TASK_SOLVE;
+    m_iparm[IPARM_END_TASK]            = API_TASK_REFINE;
+  
+    internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0,
+                           m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
   }
+  
   // Check the returned error
-  if(m_iparm(IPARM_ERROR_NUMBER)) {
-    m_info = NumericalIssue; 
-    return false;
-  }
-  else {
-    return true;
-  }
+  m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
+  
+  return m_iparm(IPARM_ERROR_NUMBER)==0;
 }
 
 /** \ingroup PaStiXSupport_Module
@@ -516,14 +469,18 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
   public:
     typedef _MatrixType MatrixType;
     typedef PastixBase<PastixLU<MatrixType> > Base;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef SparseMatrix<Scalar, ColMajor> PaStiXType; 
+    typedef typename Base::ColSpMatrix ColSpMatrix;
+    typedef typename MatrixType::Index Index;
     
   public:
-    PastixLU():Base() {}
+    PastixLU() : Base()
+    {
+      init();
+    }
     
     PastixLU(const MatrixType& matrix):Base()
     {
+      init();
       compute(matrix);
     }
     /** Compute the LU supernodal factorization of \p matrix. 
@@ -533,18 +490,9 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
       */
     void compute (const MatrixType& matrix)
     {
-      // Pastix supports only column-major matrices with a symmetric pattern
-      Base::PastixInit(); 
-      PaStiXType temp(matrix.rows(), matrix.cols());
-      // Symmetrize the graph of the matrix
-      if (IsStrSym)   
-        temp = matrix;
-      else 
-      { 
-        internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans, m_hasTranspose);
-      }
-      m_iparm[IPARM_SYM] = API_SYM_NO;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
+      m_structureIsUptodate = false;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::compute(temp);
     }
     /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern. 
@@ -554,20 +502,9 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
       */
     void analyzePattern(const MatrixType& matrix)
     {
-      
-      Base::PastixInit(); 
-      /* Pastix supports only column-major matrices with symmetrized patterns */
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      // Symmetrize the graph of the matrix
-      if (IsStrSym)   
-        temp = matrix;
-      else 
-      { 
-        internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
-      }
-      
-      m_iparm(IPARM_SYM) = API_SYM_NO;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
+      m_structureIsUptodate = false;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::analyzePattern(temp);
     }
 
@@ -578,27 +515,48 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> >
       */ 
     void factorize(const MatrixType& matrix)
     {
-      /* Pastix supports only column-major matrices with symmetrized patterns */
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      // Symmetrize the graph of the matrix
-      if (IsStrSym)   
-        temp = matrix;
-      else 
-      { 
-        internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
-      }
-      m_iparm(IPARM_SYM) = API_SYM_NO;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::factorize(temp);
     }
   protected:
+    
+    void init()
+    {
+      m_structureIsUptodate = false;
+      m_iparm(IPARM_SYM) = API_SYM_NO;
+      m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
+    }
+    
+    void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
+    {
+      if(IsStrSym)
+        out = matrix;
+      else
+      {
+        if(!m_structureIsUptodate)
+        {
+          // update the transposed structure
+          m_transposedStructure = matrix.transpose();
+          
+          // Set the elements of the matrix to zero 
+          for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 
+            for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
+              it.valueRef() = 0.0;
+
+          m_structureIsUptodate = true;
+        }
+        
+        out = m_transposedStructure + matrix;
+      }
+      internal::c_to_fortran_numbering(out);
+    }
+    
     using Base::m_iparm;
     using Base::m_dparm;
-    using Base::m_StrMatTrans;
-    using Base::m_hasTranspose;
     
-  private:
-    PastixLU(PastixLU& ) {}
+    ColSpMatrix m_transposedStructure;
+    bool m_structureIsUptodate;
 };
 
 /** \ingroup PaStiXSupport_Module
@@ -621,15 +579,18 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
   public:
     typedef _MatrixType MatrixType;
     typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::Index Index;
+    typedef typename Base::ColSpMatrix ColSpMatrix;
     
   public:
     enum { UpLo = _UpLo };
-    PastixLLT():Base() {}
+    PastixLLT() : Base()
+    {
+      init();
+    }
     
     PastixLLT(const MatrixType& matrix):Base()
     {
+      init();
       compute(matrix);
     }
 
@@ -638,13 +599,8 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
       */
     void compute (const MatrixType& matrix)
     {
-      // Pastix supports only lower, column-major matrices
-      Base::PastixInit(); // This is necessary to let PaStiX initialize its data structure between successive calls to compute
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      PermutationMatrix<Dynamic,Dynamic,Index> pnull;
-      temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
-      m_iparm(IPARM_SYM) = API_SYM_YES;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::compute(temp);
     }
 
@@ -654,13 +610,8 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
       */
     void analyzePattern(const MatrixType& matrix)
     {
-      Base::PastixInit(); 
-      // Pastix supports only lower, column-major matrices
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      PermutationMatrix<Dynamic,Dynamic,Index> pnull;
-      temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
-      m_iparm(IPARM_SYM) = API_SYM_YES;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::analyzePattern(temp);
     }
       /** Compute the LL^T supernodal numerical factorization of \p matrix 
@@ -668,19 +619,25 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
         */
     void factorize(const MatrixType& matrix)
     {
-      // Pastix supports only lower, column-major matrices 
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      PermutationMatrix<Dynamic,Dynamic,Index> pnull;
-      temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
-      m_iparm(IPARM_SYM) = API_SYM_YES;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::factorize(temp);
     }
   protected:
     using Base::m_iparm;
     
-  private:
-    PastixLLT(PastixLLT& ) {}
+    void init()
+    {
+      m_iparm(IPARM_SYM) = API_SYM_YES;
+      m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
+    }
+    
+    void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
+    {
+      // Pastix supports only lower, column-major matrices 
+      out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
+      internal::c_to_fortran_numbering(out);
+    }
 };
 
 /** \ingroup PaStiXSupport_Module
@@ -700,18 +657,21 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
 template<typename _MatrixType, int _UpLo>
 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
 {
-public:
+  public:
     typedef _MatrixType MatrixType;
     typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::Index Index;
+    typedef typename Base::ColSpMatrix ColSpMatrix;
     
   public:
     enum { UpLo = _UpLo };
-    PastixLDLT():Base() {}
+    PastixLDLT():Base()
+    {
+      init();
+    }
     
     PastixLDLT(const MatrixType& matrix):Base()
     {
+      init();
       compute(matrix);
     }
 
@@ -720,13 +680,8 @@ public:
       */
     void compute (const MatrixType& matrix)
     {
-      Base::PastixInit();
-      // Pastix supports only lower, column-major matrices
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      PermutationMatrix<Dynamic,Dynamic,Index> pnull;
-      temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
-      m_iparm(IPARM_SYM) = API_SYM_YES;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::compute(temp);
     }
 
@@ -736,14 +691,8 @@ public:
       */
     void analyzePattern(const MatrixType& matrix)
     { 
-      Base::PastixInit();
-      // Pastix supports only lower, column-major matrices
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      PermutationMatrix<Dynamic,Dynamic,Index> pnull;
-      temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
-    
-      m_iparm(IPARM_SYM) = API_SYM_YES;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::analyzePattern(temp);
     }
     /** Compute the LDL^T supernodal numerical factorization of \p matrix 
@@ -751,21 +700,26 @@ public:
       */
     void factorize(const MatrixType& matrix)
     {
-      // Pastix supports only lower, column-major matrices
-      SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
-      PermutationMatrix<Dynamic,Dynamic,Index> pnull;
-      temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
-
-      m_iparm(IPARM_SYM) = API_SYM_YES;
-      m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
+      ColSpMatrix temp;
+      grabMatrix(matrix, temp);
       Base::factorize(temp);
     }
 
   protected:
     using Base::m_iparm;
     
-  private:
-    PastixLDLT(PastixLDLT& ) {}
+    void init()
+    {
+      m_iparm(IPARM_SYM) = API_SYM_YES;
+      m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
+    }
+    
+    void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
+    {
+      // Pastix supports only lower, column-major matrices 
+      out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
+      internal::c_to_fortran_numbering(out);
+    }
 };
 
 namespace internal {

Eigen/src/SVD/JacobiSVD.h

@@ -732,8 +732,8 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, u
                           : 0);
   m_workMatrix.resize(m_diagSize, m_diagSize);
   
-  m_qr_precond_morecols.allocate(*this);
-  m_qr_precond_morerows.allocate(*this);
+  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
+  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
 }
 
 template<typename MatrixType, int QRPreconditioner>

Eigen/src/SparseCholesky/SimplicialCholesky.h

@@ -76,6 +76,9 @@ enum SimplicialCholeskyMode {
   * These classes provide LL^T and LDL^T Cholesky factorizations of sparse matrices that are
   * selfadjoint and positive definite. The factorization allows for solving A.X = B where
   * X and B can be either dense or sparse.
+  * 
+  * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
+  * such that the factorized matrix is P A P^-1.
   *
   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
   * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
@@ -208,7 +211,7 @@ class SimplicialCholeskyBase : internal::noncopyable
         return;
 
       if(m_P.size()>0)
-        dest = m_Pinv * b;
+        dest = m_P * b;
       else
         dest = b;
 
@@ -218,11 +221,11 @@ class SimplicialCholeskyBase : internal::noncopyable
       if(m_diag.size()>0)
         dest = m_diag.asDiagonal().inverse() * dest;
 
-      if (m_matrix.nonZeros()>0) // otherwise I==I
+      if (m_matrix.nonZeros()>0) // otherwise U==I
         derived().matrixU().solveInPlace(dest);
 
       if(m_P.size()>0)
-        dest = m_P * dest;
+        dest = m_Pinv * dest;
     }
 
     /** \internal */
@@ -268,7 +271,7 @@ class SimplicialCholeskyBase : internal::noncopyable
       eigen_assert(a.rows()==a.cols());
       int size = a.cols();
       CholMatrixType ap(size,size);
-      ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
+      ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
       factorize_preordered<DoLDLT>(ap);
     }
 
@@ -358,6 +361,9 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_Matr
   * This class provides a LL^T Cholesky factorizations of sparse matrices that are
   * selfadjoint and positive definite. The factorization allows for solving A.X = B where
   * X and B can be either dense or sparse.
+  * 
+  * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
+  * such that the factorized matrix is P A P^-1.
   *
   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
   * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
@@ -443,6 +449,9 @@ public:
   * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
   * selfadjoint and positive definite. The factorization allows for solving A.X = B where
   * X and B can be either dense or sparse.
+  * 
+  * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
+  * such that the factorized matrix is P A P^-1.
   *
   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
   * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
@@ -628,7 +637,7 @@ public:
         return;
 
       if(Base::m_P.size()>0)
-        dest = Base::m_Pinv * b;
+        dest = Base::m_P * b;
       else
         dest = b;
 
@@ -652,7 +661,7 @@ public:
       }
 
       if(Base::m_P.size()>0)
-        dest = Base::m_P * dest;
+        dest = Base::m_Pinv * dest;
     }
     
     Scalar determinant() const
@@ -678,22 +687,23 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixTy
   eigen_assert(a.rows()==a.cols());
   const Index size = a.rows();
   // TODO allows to configure the permutation
+  // Note that amd compute the inverse permutation
   {
     CholMatrixType C;
     C = a.template selfadjointView<UpLo>();
     // remove diagonal entries:
     // seems not to be needed
     // C.prune(keep_diag());
-    internal::minimum_degree_ordering(C, m_P);
+    internal::minimum_degree_ordering(C, m_Pinv);
   }
 
-  if(m_P.size()>0)
-    m_Pinv  = m_P.inverse();
+  if(m_Pinv.size()>0)
+    m_P = m_Pinv.inverse();
   else
-    m_Pinv.resize(0);
+    m_P.resize(0);
 
   ap.resize(size,size);
-  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
+  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
 }
 
 template<typename Derived>

Eigen/src/SparseCore/SparseMatrixBase.h

@@ -372,7 +372,7 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
     const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
     operator*(const MatrixBase<OtherDerived> &other) const;
     
-     /** \returns an expression of P^-1 H P */
+     /** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */
     SparseSymmetricPermutationProduct<Derived,Upper|Lower> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
     {
       return SparseSymmetricPermutationProduct<Derived,Upper|Lower>(derived(), perm);

Eigen/src/SparseCore/SparseSelfAdjointView.h

@@ -125,7 +125,7 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
       _dest = tmp;
     }
     
-    /** \returns an expression of P^-1 H P */
+    /** \returns an expression of P H P^-1 */
     SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
     {
       return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm);

bench/BenchTimer.h

@@ -36,7 +36,7 @@
 #   define EIGEN_BT_UNDEF_WIN32_LEAN_AND_MEAN
 # endif
 # include <windows.h>
-#elif __APPLE__
+#elif defined(__APPLE__)
 #include <CoreServices/CoreServices.h>
 #include <mach/mach_time.h>
 #else

bench/btl/data/action_settings.txt

@@ -5,6 +5,7 @@ axpby ; "{/*1.5 Y = alpha X + beta Y}" ; "vector size" ; 5:1000000
 axpy ; "{/*1.5 Y += alpha X}" ; "vector size" ; 5:1000000
 matrix_matrix ; "{/*1.5 matrix matrix product}" ; "matrix size" ; 4:3000
 matrix_vector ; "{/*1.5 matrix vector product}" ; "matrix size" ; 4:3000
+trmm ; "{/*1.5 triangular matrix matrix product}" ; "matrix size" ; 4:3000
 trisolve_vector ; "{/*1.5 triangular solver - vector (X = inv(L) X)}" ; "size" ; 4:3000
 trisolve_matrix ; "{/*1.5 triangular solver - matrix (M = inv(L) M)}" ; "size" ; 4:3000
 cholesky ; "{/*1.5 Cholesky decomposition}" ; "matrix size" ; 4:3000

bench/btl/data/go_mean

@@ -38,6 +38,7 @@ source mk_mean_script.sh atv $1 11 50 300 1000 $mode $prefix
 source mk_mean_script.sh matrix_matrix $1 11 100 300 1000 $mode $prefix
 source mk_mean_script.sh aat $1 11 100 300 1000 $mode $prefix
 # source mk_mean_script.sh ata $1 11 100 300 1000 $mode $prefix
+source mk_mean_script.sh trmm $1 11 100 300 1000 $mode $prefix
 source mk_mean_script.sh trisolve_vector $1 11 100 300 1000 $mode $prefix
 source mk_mean_script.sh trisolve_matrix $1 11 100 300 1000 $mode $prefix
 source mk_mean_script.sh cholesky $1 11 100 300 1000 $mode $prefix

bench/btl/libs/eigen3/eigen3_interface.hh

@@ -195,16 +195,16 @@ public :
   }
 
   static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){
-    X = L.template triangularView<Lower>().solve(B);
+    X = L.template triangularView<Upper>().solve(B);
   }
 
   static inline void trmm(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){
-    X = L.template triangularView<Lower>() * B;
+    X.noalias() = L.template triangularView<Lower>() * B;
   }
 
   static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){
     C = X;
-    internal::llt_inplace<Lower>::blocked(C);
+    internal::llt_inplace<real,Lower>::blocked(C);
     //C = X.llt().matrixL();
 //     C = X;
 //     Cholesky<gene_matrix>::computeInPlace(C);

bench/spbench/spbenchsolver.cpp

@@ -16,8 +16,8 @@ void bench_printhelp()
     cout<< " -h or --help \n    print this help and return\n\n";
     cout<< " -d matrixdir \n    Use matrixdir as the matrix folder instead of the one specified in the environment variable EIGEN_MATRIXDIR\n\n"; 
     cout<< " -o outputfile.html \n    Output the statistics to a html file \n\n";
-    cout<< " --eps <RelErr> Sets the relative tolerance for iterative solvers (default 1e-08)
-    cout<< " --maxits <MaxIts> Sets the maximum number of iterations (default 1000) 
+    cout<< " --eps <RelErr> Sets the relative tolerance for iterative solvers (default 1e-08) \n\n";
+    cout<< " --maxits <MaxIts> Sets the maximum number of iterations (default 1000) \n\n";
     
 }
 int main(int argc, char ** args)

blas/level3_impl.h

@@ -81,7 +81,7 @@ int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScal
 int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha,  RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
 {
 //   std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n";
-  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex);
+  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, internal::level3_blocking<Scalar,Scalar>&);
   static functype func[32];
 
   static bool init = false;
@@ -143,11 +143,17 @@ int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m,
     return xerbla_(SCALAR_SUFFIX_UP"TRSM ",&info,6);
 
   int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
-
+  
   if(SIDE(*side)==LEFT)
-    func[code](*m, *n, a, *lda, b, *ldb);
+  {
+    internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m);
+    func[code](*m, *n, a, *lda, b, *ldb, blocking);
+  }
   else
-    func[code](*n, *m, a, *lda, b, *ldb);
+  {
+    internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n);
+    func[code](*n, *m, a, *lda, b, *ldb, blocking);
+  }
 
   if(alpha!=Scalar(1))
     matrix(b,*m,*n,*ldb) *= alpha;

cmake/language_support.cmake

@@ -23,10 +23,9 @@ function(workaround_9220 language language_works)
   #message("DEBUG: language = ${language}")
   set(text
     "project(test NONE)
-cmake_minimum_required(VERSION 2.6.0)
-enable_language(${language} OPTIONAL)
-"
-    )
+    cmake_minimum_required(VERSION 2.6.0)
+    enable_language(${language} OPTIONAL)
+  ")
   file(REMOVE_RECURSE ${CMAKE_BINARY_DIR}/language_tests/${language})
   file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/language_tests/${language})
   file(WRITE ${CMAKE_BINARY_DIR}/language_tests/${language}/CMakeLists.txt

doc/C08_TutorialGeometry.dox

@@ -52,10 +52,10 @@ AngleAxis<float> aa(angle_in_radian, Vector3f(ax,ay,az));\endcode
 Quaternion<float> q;  q = AngleAxis<float>(angle_in_radian, axis);\endcode</td></tr>
 <tr class="alt"><td>
 N-D Scaling</td><td>\code
-Scaling<float,2>(sx, sy)
-Scaling<float,3>(sx, sy, sz)
-Scaling<float,N>(s)
-Scaling<float,N>(vecN)\endcode</td></tr>
+Scaling(sx, sy)
+Scaling(sx, sy, sz)
+Scaling(s)
+Scaling(vecN)\endcode</td></tr>
 <tr><td>
 N-D Translation</td><td>\code
 Translation<float,2>(tx, ty)
@@ -65,13 +65,13 @@ Translation<float,N>(vecN)\endcode</td></tr>
 <tr class="alt"><td>
 N-D \ref TutorialGeoTransform "Affine transformation"</td><td>\code
 Transform<float,N,Affine> t = concatenation_of_any_transformations;
-Transform<float,3,Affine> t = Translation3f(p) * AngleAxisf(a,axis) * Scaling3f(s);\endcode</td></tr>
+Transform<float,3,Affine> t = Translation3f(p) * AngleAxisf(a,axis) * Scaling(s);\endcode</td></tr>
 <tr><td>
 N-D Linear transformations \n
 <em class=note>(pure rotations, \n scaling, etc.)</em></td><td>\code
 Matrix<float,N> t = concatenation_of_rotations_and_scalings;
-Matrix<float,2> t = Rotation2Df(a) * Scaling2f(s);
-Matrix<float,3> t = AngleAxisf(a,axis) * Scaling3f(s);\endcode</td></tr>
+Matrix<float,2> t = Rotation2Df(a) * Scaling(s);
+Matrix<float,3> t = AngleAxisf(a,axis) * Scaling(s);\endcode</td></tr>
 </table>
 
 <strong>Notes on rotations</strong>\n To transform more than a single vector the preferred
@@ -93,8 +93,8 @@ Rotation2Df r;  r  = Matrix2f(..);       // assumes a pure rotation matrix
 AngleAxisf aa;  aa = Quaternionf(..);
 AngleAxisf aa;  aa = Matrix3f(..);       // assumes a pure rotation matrix
 Matrix2f m;     m  = Rotation2Df(..);
-Matrix3f m;     m  = Quaternionf(..);       Matrix3f m;   m = Scaling3f(..);
-Affine3f m;     m  = AngleAxis3f(..);       Affine3f m;   m = Scaling3f(..);
+Matrix3f m;     m  = Quaternionf(..);       Matrix3f m;   m = Scaling(..);
+Affine3f m;     m  = AngleAxis3f(..);       Affine3f m;   m = Scaling(..);
 Affine3f m;     m  = Translation3f(..);     Affine3f m;   m = Matrix3f(..);
 \endcode</td></tr>
 </table>
@@ -208,10 +208,10 @@ t.scale(s);
 t.prescale(Vector_(sx,sy,..));
 t.prescale(s);
 \endcode</td><td>\code
-t *= Scaling_(sx,sy,..);
-t *= Scaling_(s);
-t = Scaling_(sx,sy,..) * t;
-t = Scaling_(s) * t;
+t *= Scaling(sx,sy,..);
+t *= Scaling(s);
+t = Scaling(sx,sy,..) * t;
+t = Scaling(s) * t;
 \endcode</td></tr>
 <tr class="alt"><td>Shear transformation \n ( \b 2D \b only ! )</td><td>\code
 t.shear(sx,sy);
@@ -225,7 +225,7 @@ Note that in both API, any many transformations can be concatenated in a single
 t.pretranslate(..).rotate(..).translate(..).scale(..);
 \endcode</td></tr>
 <tr><td>\code
-t = Translation_(..) * t * RotationType(..) * Translation_(..) * Scaling_(..);
+t = Translation_(..) * t * RotationType(..) * Translation_(..) * Scaling(..);
 \endcode</td></tr>
 </table>
 

doc/C09_TutorialSparse.dox

@@ -8,9 +8,17 @@ namespace Eigen {
 
 \b Table \b of \b contents \n
   - \ref TutorialSparseIntro
+  - \ref TutorialSparseExample "Example"
+  - \ref TutorialSparseSparseMatrix
   - \ref TutorialSparseFilling
-  - \ref TutorialSparseFeatureSet
   - \ref TutorialSparseDirectSolvers
+  - \ref TutorialSparseFeatureSet
+    - \ref TutorialSparse_BasicOps
+    - \ref TutorialSparse_Products
+    - \ref TutorialSparse_TriangularSelfadjoint
+    - \ref TutorialSparse_Submat
+
+
 <hr>
 
 Manipulating and solving sparse problems involves various modules which are summarized below:
@@ -27,13 +35,14 @@ Manipulating and solving sparse problems involves various modules which are summ
 
 In many applications (e.g., finite element methods) it is common to deal with very large matrices where only a few coefficients are different from zero.  In such cases, memory consumption can be reduced and performance increased by using a specialized representation storing only the nonzero coefficients. Such a matrix is called a sparse matrix.
 
-\b The \b SparseMatrix \b class
+\b The \b %SparseMatrix \b class
+
 The class SparseMatrix is the main sparse matrix representation of Eigen's sparse module; it offers high performance and low memory usage.
 It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme.
 It consists of four compact arrays:
  - \c Values: stores the coefficient values of the non-zeros.
  - \c InnerIndices: stores the row (resp. column) indices of the non-zeros.
- - \c OuterIndexPtrs: stores for each colmun (resp. row) the index of the first non zero in the previous arrays.
+ - \c OuterStarts: stores for each column (resp. row) the index of the first non-zero in the previous two arrays.
  - \c InnerNNZs: stores the number of non-zeros of each column (resp. row).
 The word \c inner refers to an \em inner \em vector that is a column for a column-major matrix, or a row for a row-major matrix.
 The word \c outer refers to the other direction.
@@ -53,7 +62,7 @@ and one of its possible sparse, \b column \b major representation:
 <tr><td>InnerIndices:</td>  <td> 1</td><td>2</td><td>_</td><td>0</td><td>2</td><td> 4</td><td>_</td><td>_</td><td>2</td><td>_</td><td> 1</td><td>4</td></tr>
 </table>
 <table class="manual">
-<tr><td>OuterIndexPtrs:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td>\em 12 </td></tr>
+<tr><td>OuterStarts:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td>\em 12 </td></tr>
 <tr><td>InnerNNZs:</td>    <td>2</td><td>2</td><td>1</td><td>1</td><td> 2</td><td></td></tr>
 </table>
 
@@ -65,12 +74,12 @@ On the other hand, inserting elements with increasing inner indices in a given i
 The case where no empty space is available is a special case, and is refered as the \em compressed mode.
 It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS).
 Any SparseMatrix can be turned to this form by calling the SparseMatrix::makeCompressed() function.
-In this case, one can remark that the \c InnerNNZs array is redundant with \c OuterIndexPtrs because we the equality: \c InnerNNZs[j] = \c OuterIndexPtrs[j+1]-\c OuterIndexPtrs[j].
+In this case, one can remark that the \c InnerNNZs array is redundant with \c OuterStarts because we the equality: \c InnerNNZs[j] = \c OuterStarts[j+1]-\c OuterStarts[j].
 Therefore, in practice a call to SparseMatrix::makeCompressed() frees this buffer.
 
 It is worth noting that most of our wrappers to external libraries requires compressed matrices as inputs.
 
-The results of Eigen's operations always produces \b compressed sparse matrices.
+The results of %Eigen's operations always produces \b compressed sparse matrices.
 On the other hand, the insertion of a new element into a SparseMatrix converts this later to the \b uncompressed mode.
 
 Here is the previous matrix represented in compressed mode:
@@ -79,25 +88,64 @@ Here is the previous matrix represented in compressed mode:
 <tr><td>InnerIndices:</td>  <td> 1</td><td>2</td><td>0</td><td>2</td><td> 4</td><td>2</td><td> 1</td><td>4</td></tr>
 </table>
 <table class="manual">
-<tr><td>OuterIndexPtrs:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td>\em 8 </td></tr>
+<tr><td>OuterStarts:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td>\em 8 </td></tr>
 </table>
 
 A SparseVector is a special case of a SparseMatrix where only the \c Values and \c InnerIndices arrays are stored.
 There is no notion of compressed/uncompressed mode for a SparseVector.
 
 
-\b Matrix \b and \b vector \b properties \n
+\section TutorialSparseExample First example
 
-Here mat and vec represent any sparse-matrix and sparse-vector type, respectively.
+Before describing each individual class, let's start with the following typical example: solving the Lapace equation \f$ \nabla u = 0 \f$ on a regular 2D grid using a finite difference scheme and Dirichlet boundary conditions.
+Such problem can be mathematically expressed as a linear problem of the form \f$ Ax=b \f$ where \f$ x \f$ is the vector of \c m unknowns (in our case, the values of the pixels), \f$ b \f$ is the right hand side vector resulting from the boundary conditions, and \f$ A \f$ is an \f$ m \times m \f$ matrix containing only a few non-zero elements resulting from the discretization of the Laplacian operator.
+
+<table class="manual">
+<tr><td>
+\include Tutorial_sparse_example.cpp
+</td>
+<td>
+\image html Tutorial_sparse_example.jpeg
+</td></tr></table>
+
+In this example, we start by defining a column-major sparse matrix type of double \c SparseMatrix<double>, and a triplet list of the same scalar type \c  Triplet<double>. A triplet is a simple object representing a non-zero entry as the triplet: \c row index, \c column index, \c value.
+
+In the main function, we declare a list \c coefficients of triplets (as a std vector) and the right hand side vector \f$ b \f$ which are filled by the \a buildProblem function.
+The raw and flat list of non-zero entries is then converted to a true SparseMatrix object \c A.
+Note that the elements of the list do not have to be sorted, and possible duplicate entries will be summed up.
+
+The last step consists of effectively solving the assembled problem.
+Since the resulting matrix \c A is symmetric by construction, we can perform a direct Cholesky factorization via the SimplicialLDLT class which behaves like its LDLT counterpart for dense objects.
+
+The resulting vector \c x contains the pixel values as a 1D array which is saved to a jpeg file shown on the right of the code above.
+
+Describing the \a buildProblem and \a save functions is out of the scope of this tutorial. They are given \ref TutorialSparse_example_details "here" for the curious and reproducibility purpose.
+
+
+
+
+\section TutorialSparseSparseMatrix The SparseMatrix class
+
+\b %Matrix \b and \b vector \b properties \n
+
+The SparseMatrix and SparseVector classes take three template arguments:
+ * the scalar type (e.g., double)
+ * the storage order (ColMajor or RowMajor, the default is RowMajor)
+ * the inner index type (default is \c int).
+
+As for dense Matrix objects, constructors takes the size of the object.
+Here are some examples:
 
-Declarations:
 \code
-SparseMatrix<std::complex<float> > mat(1000,2000);         // declares a 1000x2000 col-major compressed sparse matrix of complex<float>
+SparseMatrix<std::complex<float> > mat(1000,2000);         // declares a 1000x2000 column-major compressed sparse matrix of complex<float>
 SparseMatrix<double,RowMajor> mat(1000,2000);              // declares a 1000x2000 row-major compressed sparse matrix of double
 SparseVector<std::complex<float> > vec(1000);              // declares a column sparse vector of complex<float> of size 1000
 SparseVector<double,RowMajor> vec(1000);                   // declares a row sparse vector of double of size 1000
 \endcode
 
+In the rest of the tutorial, \c mat and \c vec represent any sparse-matrix and sparse-vector objects, respectively.
+
+The dimensions of a matrix can be queried using the following functions:
 <table class="manual">
 <tr><td>Standard \n dimensions</td><td>\code
 mat.rows()
@@ -119,13 +167,16 @@ vec.nonZeros() \endcode</td></tr>
 
 \b Iterating \b over \b the \b nonzero \b coefficients \n
 
-Iterating over the coefficients of a sparse matrix can be done only in the same order as the storage order. Here is an example:
+Random access to the elements of a sparse object can be done through the \c coeffRef(i,j) function.
+However, this function involves a quite expensive binary search.
+In most cases, one only wants to iterate over the non-zeros elements. This is achieved by a standard loop over the outer dimension, and then by iterating over the non-zeros of the current inner vector via an InnerIterator. Thus, the non-zero entries have to be visited in the same order than the storage order.
+Here is an example:
 <table class="manual">
 <tr><td>
 \code
-SparseMatrixType mat(rows,cols);
+SparseMatrix<double> mat(rows,cols);
 for (int k=0; k<mat.outerSize(); ++k)
-  for (SparseMatrixType::InnerIterator it(mat,k); it; ++it)
+  for (SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
   {
     it.value();
     it.row();   // row index
@@ -144,22 +195,21 @@ for (SparseVector<double>::InnerIterator it(vec); it; ++it)
 \endcode
 </td></tr>
 </table>
-
+For a writable expression, the referenced value can be modified using the valueRef() function.
 If the type of the sparse matrix or vector depends on a template parameter, then the \c typename keyword is
 required to indicate that \c InnerIterator denotes a type; see \ref TopicTemplateKeyword for details.
 
 
 \section TutorialSparseFilling Filling a sparse matrix
 
-
 Because of the special storage scheme of a SparseMatrix, special care has to be taken when adding new nonzero entries.
-For instance, the cost of inserting nnz non zeros in a a single purely random insertion into a SparseMatrix is O(nnz), where nnz is the current number of nonzero coefficients.
+For instance, the cost of a single purely random insertion into a SparseMatrix is \c O(nnz), where \c nnz is the current number of non-zero coefficients.
 
-The simplest way to create a sparse matrix while guarantying good performance is to first build a list of so called \em triplets, and then convert it to a SparseMatrix.
+The simplest way to create a sparse matrix while guaranteeing good performance is thus to first build a list of so-called \em triplets, and then convert it to a SparseMatrix.
 
 Here is a typical usage example:
 \code
-typedef Triplet<double> T;
+typedef Eigen::Triplet<double> T;
 std::vector<T> tripletList;
 triplets.reserve(estimation_of_entries);
 for(...)
@@ -167,15 +217,15 @@ for(...)
   // ...
   tripletList.push_back(T(i,j,v_ij));
 }
-SparseMatrixType m(rows,cols);
-m.setFromTriplets(tripletList.begin(), tripletList.end());
-// m is ready to go!
+SparseMatrixType mat(rows,cols);
+mat.setFromTriplets(tripletList.begin(), tripletList.end());
+// mat is ready to go!
 \endcode
-The std::vector triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets().
+The \c std::vector of triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets().
 See the SparseMatrix::setFromTriplets() function and class Triplet for more details.
 
 
-In some cases, however, slightly higher performance, and lower memory consumption can be reached by directly inserting the non zeros into the destination matrix.
+In some cases, however, slightly higher performance, and lower memory consumption can be reached by directly inserting the non-zeros into the destination matrix.
 A typical scenario of this approach is illustrated bellow:
 \code
 1: SparseMatrix<double> mat(rows,cols);         // default is column major
@@ -185,61 +235,24 @@ A typical scenario of this approach is illustrated bellow:
 5: mat.makeCompressed();                        // optional
 \endcode
 
-- The key ingredient here is the line 2 where we reserve room for 6 non zeros per column. In many cases, the number of non zero per column or row can be easily known in advance. If it varies significantly for each inner vector, then it is possible to specify a reserve size for each inner vector by providing a vector object with an operator[](int j) returning the reserve size of the \c j-th inner vector (e.g., via a VectorXi or std::vector<int>). If only a rought estimate of the number of nonzeros per inner-vector can be obtained, it is highly recommended to overestimate it rather than the opposite. If this line is omitted, then the first insertion of a new element will reserve room for 2 elements per inner vector.
+- The key ingredient here is the line 2 where we reserve room for 6 non-zeros per column. In many cases, the number of non-zeros per column or row can easily be known in advance. If it varies significantly for each inner vector, then it is possible to specify a reserve size for each inner vector by providing a vector object with an operator[](int j) returning the reserve size of the \c j-th inner vector (e.g., via a VectorXi or std::vector<int>). If only a rought estimate of the number of nonzeros per inner-vector can be obtained, it is highly recommended to overestimate it rather than the opposite. If this line is omitted, then the first insertion of a new element will reserve room for 2 elements per inner vector.
 - The line 4 performs a sorted insertion. In this example, the ideal case is when the \c j-th column is not full and contains non-zeros whose inner-indices are smaller than \c i. In this case, this operation boils down to trivial O(1) operation.
 - When calling insert(i,j) the element \c i \c ,j must not already exists, otherwise use the coeffRef(i,j) method that will allow to, e.g., accumulate values. This method first performs a binary search and finally calls insert(i,j) if the element does not already exist. It is more flexible than insert() but also more costly.
 - The line 5 suppresses the remaining empty space and transforms the matrix into a compressed column storage.
 
 
-\section TutorialSparseFeatureSet Supported operators and functions
-
-In the following \em sm denotes a sparse matrix, \em sv a sparse vector, \em dm a dense matrix, and \em dv a dense vector.
-In Eigen's sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented. Moreover, not every combination is allowed; for instance, it is not possible to add two sparse matrices having two different storage orders. On the other hand, it is perfectly fine to evaluate a sparse matrix or expression to a matrix having a different storage order:
-\code
-SparseMatrixType sm1, sm2, sm3;
-sm3 = sm1.transpose() + sm2;                    // invalid, because transpose() changes the storage order
-sm3 = SparseMatrixType(sm1.transpose()) + sm2;  // correct, because evaluation reformats as column-major
-\endcode
-
-Here are some examples of supported operations:
-\code
-sm1 *= 0.5;
-sm1 = sm2 * 0.5;
-sm1 = sm2.transpose();
-sm1 = sm2.adjoint();
-sm4 = sm1 + sm2 + sm3;                          // only if sm1, sm2 and sm3 have the same storage order
-sm3 = sm1 * sm2;                                // conservative sparse * sparse product preserving numerical zeros
-sm3 = (sm1 * sm2).pruned();                     // sparse * sparse product that removes numerical zeros (triggers a different algorithm)
-sm3 = (sm1 * sm2).pruned(ref);                  // sparse * sparse product that removes elements much smaller than ref
-sm3 = (sm1 * sm2).pruned(ref,epsilon);          // sparse * sparse product that removes elements smaller than ref*epsilon
-dv3 = sm1 * dv2;
-dm3 = sm1 * dm2;
-dm3 = dm2 * sm1;
-sm3 = sm1.cwiseProduct(sm2);                    // only if sm1 and sm2 have the same storage order
-dv2 = sm1.triangularView<Upper>().solve(dv2);
-\endcode
-
-The product of a sparse \em symmetric matrix A with a dense matrix (or vector) d can be optimized by specifying the symmetry of A using selfadjointView:
-\code
-res = A.selfadjointView<>() * d;        // if all coefficients of A are stored
-res = A.selfadjointView<Upper>() * d;   // if only the upper part of A is stored
-res = A.selfadjointView<Lower>() * d;   // if only the lower part of A is stored
-\endcode
-
-
-
 \section TutorialSparseDirectSolvers Solving linear problems
 
-Eigen currently provides a limited set of built-in solvers as well as wrappers to external solver libraries.
+%Eigen currently provides a limited set of built-in solvers, as well as wrappers to external solver libraries.
 They are summarized in the following table:
 
 <table class="manual">
-<tr><td>Class</td><td>Module</td><td>Solver kind</td><td>Matrix kind</td><td>Features related to performance</td>
-    <td>Dependencies,License</td><td class="width20em"><p>Notes</p></td></tr>
-<tr><td>SimplicialLLt    </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
+<tr><th>Class</th><th>Module</th><th>Solver kind</th><th>Matrix kind</th><th>Features related to performance</th>
+    <th>Dependencies,License</th><th class="width20em"><p>Notes</p></th></tr>
+<tr><td>SimplicialLLT    </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
     <td>built-in, LGPL</td>
-    <td>SimplicialLDLt is often preferable</td></tr>
-<tr><td>SimplicialLDLt   </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LDLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
+    <td>SimplicialLDLT is often preferable</td></tr>
+<tr><td>SimplicialLDLT   </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LDLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
     <td>built-in, LGPL</td>
     <td>Recommended for very sparse and not too large problems (e.g., 2D Poisson eq.)</td></tr>
 <tr><td>ConjugateGradient</td><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>Classic iterative CG</td><td>SPD</td><td>Preconditionning</td>
@@ -250,8 +263,10 @@ They are summarized in the following table:
     <td>Might not always converge</td></tr>
 
 
-
-<tr><td>CholmodDecomposition</td><td>\link CholmodSupport_Module CholmodSupport \endlink</td><td>Direct LLT factorization</td><td>SPD</td><td>Fill-in reducing, Leverage fast dense algebra</td>
+<tr><td>PastixLLT \n PastixLDLT \n PastixLU</td><td>\link PaStiXSupport_Module PaStiXSupport \endlink</td><td>Direct LLt, LDLt, LU factorizations</td><td>SPD \n SPD \n Square</td><td>Fill-in reducing, Leverage fast dense algebra, Multithreading</td>
+    <td>Requires the <a href="http://pastix.gforge.inria.fr">PaStiX</a> package, \b CeCILL-C </td>
+    <td>optimized for tough problems and symmetric patterns</td></tr>
+<tr><td>CholmodSupernodalLLT</td><td>\link CholmodSupport_Module CholmodSupport \endlink</td><td>Direct LLt factorization</td><td>SPD</td><td>Fill-in reducing, Leverage fast dense algebra</td>
     <td>Requires the <a href="http://www.cise.ufl.edu/research/sparse/SuiteSparse/">SuiteSparse</a> package, \b GPL </td>
     <td></td></tr>
 <tr><td>UmfPackLU</td><td>\link UmfPackSupport_Module UmfPackSupport \endlink</td><td>Direct LU factorization</td><td>Square</td><td>Fill-in reducing, Leverage fast dense algebra</td>
@@ -318,6 +333,121 @@ The compute() method is equivalent to calling both analyzePattern() and factoriz
 Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on.
 More details are availble in the documentations of the respective classes.
 
+
+\section TutorialSparseFeatureSet Supported operators and functions
+
+Because of their special storage format, sparse matrices cannot offer the same level of flexbility than dense matrices.
+In Eigen's sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented.
+In the following \em sm denotes a sparse matrix, \em sv a sparse vector, \em dm a dense matrix, and \em dv a dense vector.
+
+\subsection TutorialSparse_BasicOps Basic operations
+
+%Sparse expressions support most of the unary and binary coefficient wise operations:
+\code
+sm1.real()   sm1.imag()   -sm1                    0.5*sm1
+sm1+sm2      sm1-sm2      sm1.cwiseProduct(sm2)
+\endcode
+However, a strong restriction is that the storage orders must match. For instance, in the following example:
+\code
+sm4 = sm1 + sm2 + sm3;
+\endcode
+sm1, sm2, and sm3 must all be row-major or all column major.
+On the other hand, there is no restriction on the target matrix sm4.
+For instance, this means that for computing \f$ A^T + A \f$, the matrix \f$ A^T \f$ must be evaluated into a temporary matrix of compatible storage order:
+\code
+SparseMatrix<double> A, B;
+B = SparseMatrix<double>(A.transpose()) + A;
+\endcode
+
+Binary coefficient wise operators can also mix sparse and dense expressions:
+\code
+sm2 = sm1.cwiseProduct(dm1);
+dm2 = sm1 + dm1;
+\endcode
+
+
+%Sparse expressions also support transposition:
+\code
+sm1 = sm2.transpose();
+sm1 = sm2.adjoint();
+\endcode
+However, there is no transposeInPlace() method.
+
+
+\subsection TutorialSparse_Products Matrix products
+
+%Eigen supports various kind of sparse matrix products which are summarize below:
+  - \b sparse-dense:
+    \code
+dv2 = sm1 * dv1;
+dm2 = dm1 * sm1.adjoint();
+dm2 = 2. * sm1 * dm1;
+    \endcode
+  - \b symmetric \b sparse-dense. The product of a sparse symmetric matrix with a dense matrix (or vector) can also be optimized by specifying the symmetry with selfadjointView():
+    \code
+dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
+dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
+dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored
+    \endcode
+  - \b sparse-sparse. For sparse-sparse products, two different algorithms are available. The default one is conservative and preserve the explicit zeros that might appear:
+    \code
+sm3 = sm1 * sm2;
+sm3 = 4 * sm1.adjoint() * sm2;
+    \endcode
+    The second algorithm prunes on the fly the explicit zeros, or the values smaller than a given threshold. It is enabled and controlled through the prune() functions:
+    \code
+sm3 = (sm1 * sm2).prune();                  // removes numerical zeros
+sm3 = (sm1 * sm2).prune(ref);               // removes elements much smaller than ref
+sm3 = (sm1 * sm2).prune(ref,epsilon);       // removes elements smaller than ref*epsilon
+    \endcode
+
+  - \b permutations. Finally, permutations can be applied to sparse matrices too:
+    \code
+PermutationMatrix<Dynamic,Dynamic> P = ...;
+sm2 = P * sm1;
+sm2 = sm1 * P.inverse();
+sm2 = sm1.transpose() * P;
+    \endcode
+
+
+\subsection TutorialSparse_TriangularSelfadjoint Triangular and selfadjoint views
+
+Just as with dense matrices, the triangularView() function can be used to address a triangular part of the matrix, and perform triangular solves with a dense right hand side:
+\code
+dm2 = sm1.triangularView<Lower>(dm1);
+dv2 = sm1.transpose().triangularView<Upper>(dv1);
+\endcode
+
+The selfadjointView() function permits various operations:
+ - optimized sparse-dense matrix products:
+    \code
+dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
+dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
+dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored
+    \endcode
+ - copy of triangular parts:
+    \code
+sm2 = sm1.selfadjointView<Upper>();                               // makes a full selfadjoint matrix from the upper triangular part
+sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>();      // copies the upper triangular part to the lower triangular part
+    \endcode
+ - application of symmetric permutations:
+ \code
+PermutationMatrix<Dynamic,Dynamic> P = ...;
+sm2 = A.selfadjointView<Upper>().twistedBy(P);                                // compute P S P' from the upper triangular part of A, and make it a full matrix
+sm2.selfadjointView<Lower>() = A.selfadjointView<Lower>().twistedBy(P);       // compute P S P' from the lower triangular part of A, and then only compute the lower part
+ \endcode
+
+\subsection TutorialSparse_Submat Sub-matrices
+
+%Sparse matrices does not support yet the addressing of arbitrary sub matrices. Currently, one can only reference a set of contiguous \em inner vectors, i.e., a set of contiguous rows for a row-major matrix, or a set of contiguous columns for a column major matrix:
+\code
+  sm1.innerVector(j);       // returns an expression of the j-th column (resp. row) of the matrix if sm1 is col-major (resp. row-major)
+  sm1.innerVectors(j, nb);  // returns an expression of the nb columns (resp. row) starting from the j-th column (resp. row)
+                            // of the matrix if sm1 is col-major (resp. row-major)
+  sm1.middleRows(j, nb);    // for row major matrices only, get a range of nb rows
+  sm1.middleCols(j, nb);    // for column major matrices only, get a range of nb columns
+\endcode
+
 \li \b Next: \ref TutorialMapClass
 
 */

doc/CMakeLists.txt

@@ -36,6 +36,7 @@ set(snippets_targets "")
 add_definitions("-DEIGEN_MAKING_DOCS")
 
 add_subdirectory(examples)
+add_subdirectory(special_examples)
 add_subdirectory(snippets)
 
 add_custom_target(

doc/Doxyfile.in

@@ -592,6 +592,7 @@ RECURSIVE              = YES
 EXCLUDE                = "${Eigen_SOURCE_DIR}/Eigen/Eigen2Support" \
                          "${Eigen_SOURCE_DIR}/Eigen/src/Eigen2Support" \
                          "${Eigen_SOURCE_DIR}/doc/examples" \
+                         "${Eigen_SOURCE_DIR}/doc/special_examples" \
                          "${Eigen_SOURCE_DIR}/doc/snippets" 
 
 # The EXCLUDE_SYMLINKS tag can be used select whether or not files or
@@ -638,7 +639,9 @@ EXCLUDE_SYMBOLS        = internal::* Flagged* *InnerIterator* DenseStorage<*
 EXAMPLE_PATH           = "${Eigen_SOURCE_DIR}/doc/snippets" \
                          "${Eigen_BINARY_DIR}/doc/snippets" \
                          "${Eigen_SOURCE_DIR}/doc/examples" \
-                         "${Eigen_BINARY_DIR}/doc/examples"
+                         "${Eigen_BINARY_DIR}/doc/examples" \
+                         "${Eigen_SOURCE_DIR}/doc/special_examples" \
+                         "${Eigen_BINARY_DIR}/doc/special_examples"
 
 # If the value of the EXAMPLE_PATH tag contains directories, you can use the
 # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp

doc/I00_CustomizingEigen.dox

@@ -120,19 +120,22 @@ Eigen::MatrixBase<Eigen::Matrix<std::complex<float>, 10000, 1, 2, 10000, 1>
 
 \anchor user_defined_scalars \section CustomScalarType Using custom scalar types
 
-By default, Eigen currently supports the following scalar types: \c int, \c float, \c double, \c std::complex<float>, \c std::complex<double>, \c long \c double, \c long \c long \c int (64 bits integers), and \c bool. The \c long \c double is especially useful on x86-64 systems or when the SSE2 instruction set is enabled because it enforces the use of x87 registers with extended accuracy.
+By default, Eigen currently supports standard floating-point types (\c float, \c double, \c std::complex<float>, \c std::complex<double>, \c long \c double), as well as all integrale types (e.g., \c int, \c unsigned \c int, \c short, etc.), and \c bool.
+On x86-64 systems, \c long \c double permits to locally enforces the use of x87 registers with extended accuracy (in comparison to SSE).
 
 In order to add support for a custom type \c T you need:
- 1 - make sure the common operator (+,-,*,/,etc.) are supported by the type \c T
- 2 - add a specialization of struct Eigen::NumTraits<T> (see \ref NumTraits)
- 3 - define a couple of math functions for your type such as: internal::sqrt, internal::abs, etc...
+-# make sure the common operator (+,-,*,/,etc.) are supported by the type \c T
+-# add a specialization of struct Eigen::NumTraits<T> (see \ref NumTraits)
+-# define the math functions that makes sense for your type. This includes standard ones like sqrt, pow, sin, tan, conj, real, imag, etc, as well as abs2 which is Eigen specific.
      (see the file Eigen/src/Core/MathFunctions.h)
 
+The math function should be defined in the same namespace than \c T, or in the \c std namespace though that second appraoch is not recommended.
+
 Here is a concrete example adding support for the Adolc's \c adouble type. <a href="https://projects.coin-or.org/ADOL-C">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives.
 
 \code
-#ifndef ADLOCSUPPORT_H
-#define ADLOCSUPPORT_H
+#ifndef ADOLCSUPPORT_H
+#define ADOLCSUPPORT_H
 
 #define ADOLC_TAPELESS
 #include <adolc/adouble.h>
@@ -141,6 +144,7 @@ Here is a concrete example adding support for the Adolc's \c adouble type. <a hr
 namespace Eigen {
 
 template<> struct NumTraits<adtl::adouble>
+ : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
 {
   typedef adtl::adouble Real;
   typedef adtl::adouble NonInteger;
@@ -149,35 +153,27 @@ template<> struct NumTraits<adtl::adouble>
   enum {
     IsComplex = 0,
     IsInteger = 0,
-    IsSigned,
+    IsSigned = 1,
+    RequireInitialization = 1,
     ReadCost = 1,
-    AddCost = 1,
-    MulCost = 1
+    AddCost = 3,
+    MulCost = 3
   };
 };
 
 }
 
-// the Adolc's type adouble is defined in the adtl namespace
-// therefore, the following internal::* functions *must* be defined
-// in the same namespace
 namespace adtl {
 
-  inline const adouble& internal::conj(const adouble& x)  { return x; }
-  inline const adouble& internal::real(const adouble& x)  { return x; }
-  inline adouble internal::imag(const adouble&)    { return 0.; }
-  inline adouble internal::abs(const adouble&  x)  { return fabs(x); }
-  inline adouble internal::abs2(const adouble& x)  { return x*x; }
-  inline adouble internal::sqrt(const adouble& x)  { return sqrt(x); }
-  inline adouble internal::exp(const adouble&  x)  { return exp(x); }
-  inline adouble internal::log(const adouble&  x)  { return log(x); }
-  inline adouble internal::sin(const adouble&  x)  { return sin(x); }
-  inline adouble internal::cos(const adouble&  x)  { return cos(x); }
-  inline adouble internal::pow(const adouble& x, adouble y)  { return pow(x, y); }
+inline const adouble& conj(const adouble& x)  { return x; }
+inline const adouble& real(const adouble& x)  { return x; }
+inline adouble imag(const adouble&)    { return 0.; }
+inline adouble abs(const adouble&  x)  { return fabs(x); }
+inline adouble abs2(const adouble& x)  { return x*x; }
 
 }
 
-#endif // ADLOCSUPPORT_H
+#endif // ADOLCSUPPORT_H
 \endcode
 
 

doc/I09_Vectorization.dox

@@ -1,6 +1,6 @@
 namespace Eigen {
 
-/** \page TopicVectorization Vectorizaion
+/** \page TopicVectorization Vectorization
 
 
 TODO: write this dox page!

doc/Overview.dox

@@ -34,6 +34,7 @@ For a first contact with Eigen, the best place is to have a look at the \ref Get
     - \ref TopicLazyEvaluation
     - \ref TopicLinearAlgebraDecompositions
     - \ref TopicCustomizingEigen
+    - \ref TopicMultiThreading
     - \ref TopicPreprocessorDirectives
     - \ref TopicStorageOrders
     - \ref TopicInsideEigenExample

doc/SparseQuickReference.dox

@@ -18,7 +18,7 @@ In this page, we give a quick summary of the main operations available for spars
 i.e either row major or column major. The default is column major. Most arithmetic operations on sparse matrices will assert that they have the same storage order. Moreover, when interacting with external libraries that are not yet supported by Eigen, it is important to know how to send the required matrix pointers. 
 
 \section Constructors Constructors and assignments
-SparseMatrix is the core class to build and manipulate sparse matrices in Eigen. It takes as template parameters the Scalar type and the storage order, either RowMajor or ColumnMajor. The default is ColumnMajor. ??? It is possible to modify the default storage order at compile-time with the cmake variable \b EIGEN_DEFAULT_ROW_MAJOR ???
+SparseMatrix is the core class to build and manipulate sparse matrices in Eigen. It takes as template parameters the Scalar type and the storage order, either RowMajor or ColumnMajor. The default is ColumnMajor.
 
 \code
   SparseMatrix<double> sm1(1000,1000);              // 1000x1000 compressed sparse matrix of double. 
@@ -195,4 +195,4 @@ The following functions are useful to extract a block of rows (resp. columns) fr
 
 
 */
-}
+}

doc/TopicMultithreading.dox

@@ -0,0 +1,46 @@
+namespace Eigen {
+
+/** \page TopicMultiThreading Eigen and multi-threading
+
+\section TopicMultiThreading_MakingEigenMT Make Eigen run in parallel
+
+Some Eigen's algorithms can exploit the multiple cores present in your hardware. To this end, it is enough to enable OpenMP on your compiler, for instance:
+ * GCC: \c -fopenmp
+ * ICC: \c -openmp
+ * MSVC: check the respective option in the build properties.
+You can control the number of thread that will be used using either the OpenMP API or Eiegn's API using the following priority:
+\code
+ OMP_NUM_THREADS=n ./my_program
+ omp_set_num_threads(n);
+ Eigen::setNbThreads(n);
+\endcode
+Unless setNbThreads has been called, Eigen uses the number of threads specified by OpenMP. You can restore this bahavior by calling \code setNbThreads(0); \endcode
+You can query the number of threads that will be used with:
+\code
+n = Eigen::nbThreads(n);
+\endcode
+You can disable Eigen's multi threading at compile time by defining the EIGEN_DONT_PARALLELIZE preprocessor token.
+
+Currently, the following algorithms can make use of multi-threading:
+ * general matrix - matrix products
+ * PartialPivLU
+
+\section TopicMultiThreading_UsingEigenWithMT Using Eigen in a multi-threaded application
+
+In the case your own application is multithreaded, and multiple threads make calls to Eigen, then you have to initialize Eigen by calling the following routine \b before creating the threads:
+\code
+#include <Eigen/Core>
+
+int main(int argc, char** argv)
+{
+  Eigen::initParallel();
+  
+  ...
+}
+\endcode
+
+In the case your application is parallelized with OpenMP, you might want to disable Eigen's own parallization as detailed in the previous section.
+
+*/
+
+}

doc/TutorialSparse_example_details.dox

@@ -0,0 +1,4 @@
+/**
+\page TutorialSparse_example_details
+\include Tutorial_sparse_example_details.cpp
+*/

doc/UsingIntelMKL.dox

@@ -45,6 +45,7 @@ Intel MKL is available on Linux, Mac and Windows for both Intel64 and IA32 archi
 Using Intel MKL through Eigen is easy:
 -# define the \c EIGEN_USE_MKL_ALL macro before including any Eigen's header
 -# link your program to MKL libraries (see the <a href="http://software.intel.com/en-us/articles/intel-mkl-link-line-advisor/">MKL linking advisor</a>)
+-# on a 64bits system, you must use the LP64 interface (not the ILP64 one)
 
 When doing so, a number of Eigen's algorithms are silently substituted with calls to Intel MKL routines.
 These substitutions apply only for \b Dynamic \b or \b large enough objects with one of the following four standard scalar types: \c float, \c double, \c complex<float>, and \c complex<double>.

doc/special_examples/CMakeLists.txt

@@ -0,0 +1,20 @@
+
+if(NOT EIGEN_TEST_NOQT)
+  find_package(Qt4)
+  if(QT4_FOUND)
+    include(${QT_USE_FILE})
+  endif()
+endif(NOT EIGEN_TEST_NOQT)
+
+
+if(QT4_FOUND)
+  add_executable(Tutorial_sparse_example Tutorial_sparse_example.cpp Tutorial_sparse_example_details.cpp)
+  target_link_libraries(Tutorial_sparse_example ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO} ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY})
+
+  add_custom_command(
+      TARGET Tutorial_sparse_example
+      POST_BUILD
+      COMMAND Tutorial_sparse_example
+      ARGS ${CMAKE_CURRENT_BINARY_DIR}/../html/Tutorial_sparse_example.jpeg
+  )
+endif(QT4_FOUND)

doc/special_examples/Tutorial_sparse_example.cpp

@@ -0,0 +1,32 @@
+#include <Eigen/Sparse>
+#include <vector>
+
+typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
+typedef Eigen::Triplet<double> T;
+
+void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n);
+void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename);
+
+int main(int argc, char** argv)
+{
+  int n = 300;  // size of the image
+  int m = n*n;  // number of unknows (=number of pixels)
+
+  // Assembly:
+  std::vector<T> coefficients;            // list of non-zeros coefficients
+  Eigen::VectorXd b(m);                   // the right hand side-vector resulting from the constraints
+  buildProblem(coefficients, b, n);
+
+  SpMat A(m,m);
+  A.setFromTriplets(coefficients.begin(), coefficients.end());
+
+  // Solving:
+  Eigen::SimplicialCholesky<SpMat> chol(A);  // performs a Cholesky factorization of A
+  Eigen::VectorXd x = chol.solve(b);         // use the factorization to solve for the given right hand side
+
+  // Export the result to a file:
+  saveAsBitmap(x, n, argv[1]);
+
+  return 0;
+}
+

doc/special_examples/Tutorial_sparse_example_details.cpp

@@ -0,0 +1,44 @@
+#include <Eigen/Sparse>
+#include <vector>
+#include <QImage>
+
+typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
+typedef Eigen::Triplet<double> T;
+
+void insertCoefficient(int id, int i, int j, double w, std::vector<T>& coeffs,
+                       Eigen::VectorXd& b, const Eigen::VectorXd& boundary)
+{
+  int n = boundary.size();
+  int id1 = i+j*n;
+
+        if(i==-1 || i==n) b(id) -= w * boundary(j); // constrained coeffcieint
+  else  if(j==-1 || j==n) b(id) -= w * boundary(i); // constrained coeffcieint
+  else  coeffs.push_back(T(id,id1,w));              // unknown coefficient
+}
+
+void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n)
+{
+  b.setZero();
+  Eigen::ArrayXd boundary = Eigen::ArrayXd::LinSpaced(n, 0,M_PI).sin().pow(2);
+  for(int j=0; j<n; ++j)
+  {
+    for(int i=0; i<n; ++i)
+    {
+      int id = i+j*n;
+      insertCoefficient(id, i-1,j, -1, coefficients, b, boundary);
+      insertCoefficient(id, i+1,j, -1, coefficients, b, boundary);
+      insertCoefficient(id, i,j-1, -1, coefficients, b, boundary);
+      insertCoefficient(id, i,j+1, -1, coefficients, b, boundary);
+      insertCoefficient(id, i,j,    4, coefficients, b, boundary);
+    }
+  }
+}
+
+void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename)
+{
+  Eigen::Array<unsigned char,Eigen::Dynamic,Eigen::Dynamic> bits = (x*255).cast<unsigned char>();
+  QImage img(bits.data(), n,n,QImage::Format_Indexed8);
+  img.setColorCount(256);
+  for(int i=0;i<256;i++) img.setColor(i,qRgb(i,i,i));
+  img.save(filename);
+}

test/CMakeLists.txt

@@ -18,13 +18,15 @@ set(LAPACK_FOUND TRUE)
 set(BLAS_LIBRARIES eigen_blas)
 set(LAPACK_LIBRARIES eigen_lapack)
 
-if(TEST_REAL_CASES)
+set(EIGEN_TEST_MATRIX_DIR "" CACHE STRING "Enable testing of realword sparse matrices contained in the specified path")
+if(EIGEN_TEST_MATRIX_DIR)
   if(NOT WIN32)
-    add_definitions( -DTEST_REAL_CASES="${TEST_REAL_CASES}" )
+    message(STATUS "Test realworld sparse matrices: ${EIGEN_TEST_MATRIX_DIR}")
+    add_definitions( -DTEST_REAL_CASES="${EIGEN_TEST_MATRIX_DIR}" )
   else(NOT WIN32)
-    message(STATUS, "REAL CASES CAN NOT BE CURRENTLY TESTED ON WIN32")
+    message(STATUS "REAL CASES CAN NOT BE CURRENTLY TESTED ON WIN32")
   endif(NOT WIN32)
-endif(TEST_REAL_CASES)
+endif(EIGEN_TEST_MATRIX_DIR)
 
 set(SPARSE_LIBS " ")
 

test/nomalloc.cpp

@@ -129,13 +129,20 @@ void ctms_decompositions()
                         0,
                         maxSize, maxSize> ComplexMatrix;
 
-  const Matrix A(Matrix::Random(size, size));
+  const Matrix A(Matrix::Random(size, size)), B(Matrix::Random(size, size));
+  Matrix X(size,size);
   const ComplexMatrix complexA(ComplexMatrix::Random(size, size));
   const Matrix saA = A.adjoint() * A;
+  const Vector b(Vector::Random(size));
+  Vector x(size);
 
   // Cholesky module
   Eigen::LLT<Matrix>  LLT;  LLT.compute(A);
+  X = LLT.solve(B);
+  x = LLT.solve(b);
   Eigen::LDLT<Matrix> LDLT; LDLT.compute(A);
+  X = LDLT.solve(B);
+  x = LDLT.solve(b);
 
   // Eigenvalues module
   Eigen::HessenbergDecomposition<ComplexMatrix> hessDecomp;        hessDecomp.compute(complexA);
@@ -147,12 +154,22 @@ void ctms_decompositions()
 
   // LU module
   Eigen::PartialPivLU<Matrix> ppLU; ppLU.compute(A);
+  X = ppLU.solve(B);
+  x = ppLU.solve(b);
   Eigen::FullPivLU<Matrix>    fpLU; fpLU.compute(A);
+  X = fpLU.solve(B);
+  x = fpLU.solve(b);
 
   // QR module
   Eigen::HouseholderQR<Matrix>        hQR;  hQR.compute(A);
+  X = hQR.solve(B);
+  x = hQR.solve(b);
   Eigen::ColPivHouseholderQR<Matrix>  cpQR; cpQR.compute(A);
+  // FIXME X = cpQR.solve(B);
+  x = cpQR.solve(b);
   Eigen::FullPivHouseholderQR<Matrix> fpQR; fpQR.compute(A);
+  // FIXME X = fpQR.solve(B);
+  x = fpQR.solve(b);
 
   // SVD module
   Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV);

test/sparse.h

@@ -79,7 +79,7 @@ initSparse(double density,
   enum { IsRowMajor = SparseMatrix<Scalar,Opt2,Index>::IsRowMajor };
   sparseMat.setZero();
   //sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
-  sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), (1.5*density)*(IsRowMajor?refMat.cols():refMat.rows())));
+  sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), int((1.5*density)*(IsRowMajor?refMat.cols():refMat.rows()))));
   
   for(int j=0; j<sparseMat.outerSize(); j++)
   {

unsupported/Eigen/AdolcForward

@@ -79,33 +79,22 @@ namespace Eigen {
 
 } // namespace Eigen
 
-// the Adolc's type adouble is defined in the adtl namespace
-// therefore, the following internal::* functions *must* be defined
-// in the same namespace
-namespace Eigen {
+// Eigen's require a few additional functions which must be defined in the same namespace
+// than the custom scalar type own namespace
+namespace adtl {
 
-  namespace internal {
-
-    inline const adtl::adouble& conj(const adtl::adouble& x)  { return x; }
-    inline const adtl::adouble& real(const adtl::adouble& x)  { return x; }
-    inline adtl::adouble imag(const adtl::adouble&)    { return 0.; }
-    inline adtl::adouble abs(const adtl::adouble&  x)  { return adtl::fabs(x); }
-    inline adtl::adouble abs2(const adtl::adouble& x)  { return x*x; }
-    
-    using adtl::sqrt;
-    using adtl::exp;
-    using adtl::log;
-    using adtl::sin;
-    using adtl::cos;
-    using adtl::pow;
-
-  }
+inline const adouble& conj(const adouble& x)  { return x; }
+inline const adouble& real(const adouble& x)  { return x; }
+inline adouble imag(const adouble&)    { return 0.; }
+inline adouble abs(const adouble&  x)  { return fabs(x); }
+inline adouble abs2(const adouble& x)  { return x*x; }
 
 }
 
 namespace Eigen {
 
 template<> struct NumTraits<adtl::adouble>
+    : NumTraits<double>
 {
   typedef adtl::adouble Real;
   typedef adtl::adouble NonInteger;

unsupported/Eigen/MPRealSupport

@@ -107,7 +107,7 @@ int main()
     }
   };
 
-  namespace internal {
+namespace internal {
 
   template<> mpfr::mpreal random<mpfr::mpreal>()
   {
@@ -133,18 +133,6 @@ int main()
     return a + (b-a) * random<mpfr::mpreal>();
   }
 
-  template<> struct conj_impl<mpfr::mpreal> { inline static const mpfr::mpreal& run(const mpfr::mpreal& x) { return x; } };
-  template<> struct real_impl<mpfr::mpreal> { inline static const mpfr::mpreal& run(const mpfr::mpreal& x) { return x; } };
-  template<> struct imag_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal&)   { return mpfr::mpreal(0); } };
-  template<> struct abs_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x)  { return mpfr::fabs(x); } };
-  template<> struct abs2_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x) { return x*x; } };
-  template<> struct sqrt_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x) { return mpfr::sqrt(x); } };
-  template<> struct exp_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x)  { return mpfr::exp(x); } };
-  template<> struct log_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x)  { return mpfr::log(x); } };
-  template<> struct sin_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x)  { return mpfr::sin(x); } };
-  template<> struct cos_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x)  { return mpfr::cos(x); } };
-  template<> struct pow_impl<mpfr::mpreal> { inline static const mpfr::mpreal run(const mpfr::mpreal& x, const mpfr::mpreal& y)  { return mpfr::pow(x, y); } };
-
   bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec)
   {
     return mpfr::abs(a) <= mpfr::abs(b) * prec;
@@ -159,8 +147,17 @@ int main()
   {
     return a <= b || isApprox(a, b, prec);
   }
+  
+  template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x)
+  { return x.toLDouble(); }
+  template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x)
+  { return x.toDouble(); }
+  template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x)
+  { return x.toLong(); }
+  template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x)
+  { return int(x.toLong()); }
 
-  } // end namespace internal
+} // end namespace internal
 }
 
 #endif // EIGEN_MPREALSUPPORT_MODULE_H

unsupported/Eigen/src/NumericalDiff/NumericalDiff.h

@@ -66,7 +66,7 @@ public:
     template<typename T0, typename T1>
         NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {}
     template<typename T0, typename T1, typename T2>
-        NumericalDiff(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2), epsfcn(0) {}
+        NumericalDiff(const T0& a0, const T1& a1, const T2& a2) : Functor(a0, a1, a2), epsfcn(0) {}
 
     enum {
         InputsAtCompileTime = Functor::InputsAtCompileTime,

unsupported/Eigen/src/SparseExtra/MarketIO.h

@@ -156,6 +156,9 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
   char buffer[maxBuffersize];
   
   bool readsizes = false;
+
+  typedef Triplet<Scalar,int> T;
+  std::vector<T> elements;
   
   int M(-1), N(-1), NNZ(-1);
   int count = 0;
@@ -186,13 +189,13 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
       if( internal::GetMarketLine(line, M, N, i, j, value) ) 
       {
         ++ count;
-        mat.insert(i,j) = value;
+        elements.push_back(T(i,j,value));
       }
       else 
         std::cerr << "Invalid read: " << i << "," << j << "\n";        
     }
   }
-  mat.makeCompressed();
+  mat.setFromTriplets(elements.begin(), elements.end());
   if(count!=NNZ)
     std::cerr << count << "!=" << NNZ << "\n";
   

unsupported/Eigen/src/SparseExtra/MatrixMarketIterator.h

@@ -85,7 +85,7 @@ class MatrixMarketIterator
       Getnextvalidmatrix();
       return *this;
     }
-    inline operator bool() { return m_isvalid;}
+    inline operator bool() const { return m_isvalid;}
     
     /** Return the sparse matrix corresponding to the current file */
     inline MatrixType& matrix() 
@@ -186,6 +186,7 @@ class MatrixMarketIterator
     
     void Getnextvalidmatrix( )
     {
+      m_isvalid = false;
       // Here, we return with the next valid matrix in the folder
       while ( (m_curs_id = readdir(m_folder_id)) != NULL) {
         m_isvalid = false;

unsupported/test/forward_adolc.cpp

@@ -23,11 +23,20 @@
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
 #include "main.h"
+#include <Eigen/Dense>
+
 #define NUMBER_DIRECTIONS 16
 #include <unsupported/Eigen/AdolcForward>
 
 int adtl::ADOLC_numDir;
 
+template<typename Vector>
+EIGEN_DONT_INLINE typename Vector::Scalar foo(const Vector& p)
+{
+  typedef typename Vector::Scalar Scalar;
+  return (p-Vector(Scalar(-1),Scalar(1.))).norm() + (p.array().sqrt().abs() * p.array().sin()).sum() + p.dot(p);
+}
+
 template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
 struct TestFunc1
 {
@@ -138,4 +147,12 @@ void test_forward_adolc()
     CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) ));
     CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double>(3,3)) ));
   }
+
+  {
+    // simple instanciation tests
+    Matrix<adtl::adouble,2,1> x;
+    foo(x);
+    Matrix<adtl::adouble,Dynamic,Dynamic> A(4,4);;
+    A.selfadjointView<Lower>().eigenvalues();
+  }
 }

unsupported/test/mpreal/dlmalloc.h

@@ -556,7 +556,7 @@ int mspace_mallopt(int, int);
 #endif  /* MSPACES */
 
 #ifdef __cplusplus
-}  /* end of extern "C" */
+};  /* end of extern "C" */
 #endif
 
 #endif /* MALLOC_280_H */

unsupported/test/mpreal/mpreal.cpp

@@ -3,14 +3,15 @@
 	Project homepage: http://www.holoborodko.com/pavel/
 	Contact e-mail:   pavel@holoborodko.com
 
-	Copyright (c) 2008-2010 Pavel Holoborodko
+	Copyright (c) 2008-2011 Pavel Holoborodko
 
 	Core Developers: 
 	Pavel Holoborodko, Dmitriy Gubanov, Konstantin Holoborodko. 
 
 	Contributors:
 	Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, 
-	Heinz van Saanen, Pere Constans, Peter van Hoof.
+	Heinz van Saanen, Pere Constans, Peter van Hoof, Gael Guennebaud, 
+	Tsai Chia Cheng, Alexei Zubanov.
 
 	****************************************************************************
 	This library is free software; you can redistribute it and/or
@@ -27,31 +28,21 @@
 	License along with this library; if not, write to the Free Software
 	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 
+	****************************************************************************
 	****************************************************************************
 	Redistribution and use in source and binary forms, with or without
 	modification, are permitted provided that the following conditions
 	are met:
-	
+
 	1. Redistributions of source code must retain the above copyright
 	notice, this list of conditions and the following disclaimer.
-	
+
 	2. Redistributions in binary form must reproduce the above copyright
 	notice, this list of conditions and the following disclaimer in the
 	documentation and/or other materials provided with the distribution.
-	
-	3. Redistributions of any form whatsoever must retain the following
-	acknowledgment:
-	"
-         This product includes software developed by Pavel Holoborodko
-         Web: http://www.holoborodko.com/pavel/
-         e-mail: pavel@holoborodko.com
-	"
 
-	4. This software cannot be, by any means, used for any commercial 
-	purpose without the prior permission of the copyright holder.
-	
-	Any of the above conditions can be waived if you get permission from 
-	the copyright holder. 
+	3. The name of the author may be used to endorse or promote products
+	derived from this software without specific prior written permission.
 
 	THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 	ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
@@ -66,9 +57,11 @@
 	SUCH DAMAGE.
 */
 #include <cstring>
-#include <cstdlib>
 #include "mpreal.h"
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 #include "dlmalloc.h"
+#endif
 
 using std::ws;
 using std::cerr;
@@ -79,62 +72,107 @@ using std::istream;
 
 namespace mpfr{
 
-mp_rnd_t   mpreal::default_rnd  = mpfr_get_default_rounding_mode();	
-mp_prec_t  mpreal::default_prec = mpfr_get_default_prec();	
+mp_rnd_t   mpreal::default_rnd  = MPFR_RNDN;	//(mpfr_get_default_rounding_mode)();	
+mp_prec_t  mpreal::default_prec = 64;			//(mpfr_get_default_prec)();	
 int		   mpreal::default_base = 10;
 int        mpreal::double_bits = -1;
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 bool       mpreal::is_custom_malloc = false;
+#endif
 
 // Default constructor: creates mp number and initializes it to 0.
 mpreal::mpreal() 
 { 
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,default_prec); 
 	mpfr_set_ui(mp,0,default_rnd);
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const mpreal& u) 
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,mpfr_get_prec(u.mp));
 	mpfr_set(mp,u.mp,default_rnd);
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const mpfr_t u)
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,mpfr_get_prec(u));
 	mpfr_set(mp,u,default_rnd);
+	
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const mpf_t u)
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
-	mpfr_init2(mp,mpf_get_prec(u));
+#endif
+
+	mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
 	mpfr_set_f(mp,u,default_rnd);
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_z(mp,u,mode);
+	
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_q(mp,u,mode);
+	
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
     if(double_bits == -1 || fits_in_bits(u, double_bits))
     {
     	mpfr_init2(mp,prec);
 	    mpfr_set_d(mp,u,mode);
+		
+		MPREAL_MSVC_DEBUGVIEW_CODE;
     }
     else
         throw conversion_overflow();
@@ -142,51 +180,121 @@ mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
 
 mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
 { 
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
     mpfr_init2(mp,prec);
 	mpfr_set_ld(mp,u,mode);
+	
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
 { 
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_ui(mp,u,mode);
+	
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
 { 
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_ui(mp,u,mode);
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
 { 
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_si(mp,u,mode);
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
 { 
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_si(mp,u,mode);
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
+#if defined (MPREAL_HAVE_INT64_SUPPORT)
+mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode)
+{
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
+	set_custom_malloc();
+#endif
+
+	mpfr_init2(mp,prec);
+	mpfr_set_uj(mp, u, mode); 
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)
+{
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
+	set_custom_malloc();
+#endif
+
+	mpfr_init2(mp,prec);
+	mpfr_set_sj(mp, u, mode); 
+	
+	MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+#endif
+
 mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_str(mp, s, base, mode); 
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)
 {
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
+
 	mpfr_init2(mp,prec);
 	mpfr_set_str(mp, s.c_str(), base, mode); 
+
+	MPREAL_MSVC_DEBUGVIEW_CODE;
 }
 
 mpreal::~mpreal() 
@@ -198,18 +306,22 @@ mpreal::~mpreal()
 mpreal& mpreal::operator=(const char* s)
 {
 	mpfr_t t;
-	
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
+#endif
 
 	if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
 	{
-		// We will rewrite mp anyway, so use flash it and resize
-		mpfr_set_prec(mp,mpfr_get_prec(t)); //<- added 01.04.2011
+		// We will rewrite mp anyway, so flash it and resize
+		mpfr_set_prec(mp,mpfr_get_prec(t)); 
 		mpfr_set(mp,t,mpreal::default_rnd);
 		mpfr_clear(t);
+
+		MPREAL_MSVC_DEBUGVIEW_CODE;
+
 	}else{
 		mpfr_clear(t);
-		// cerr<<"fail to convert string"<<endl;
 	}
 
 	return *this;
@@ -260,21 +372,6 @@ const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
 	return a;
 }
 
-const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
-{
-	mpreal a;
-	mp_prec_t yp, xp;
-
-	yp = y.get_prec(); 
-	xp = x.get_prec(); 
-
-	a.set_prec(yp>xp?yp:xp);
-
-	mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);
-
-	return a;
-}
-
 const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
 {
 	mpreal x;
@@ -288,21 +385,6 @@ const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
 	return x;
 }
 
-const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
-{	
-	mpreal a;
-	mp_prec_t yp, xp;
-
-	yp = y.get_prec(); 
-	xp = x.get_prec(); 
-
-	a.set_prec(yp>xp?yp:xp);
-
-	mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);
-
-	return a;
-}
-
 const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
 {
 	mpreal a;
@@ -319,36 +401,71 @@ const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mod
 }
 
 template <class T>
-std::string to_string(T t, std::ios_base & (*f)(std::ios_base&))
+std::string toString(T t, std::ios_base & (*f)(std::ios_base&))
 {
 	std::ostringstream oss;
 	oss << f << t;
 	return oss.str();
 }
 
-mpreal::operator std::string() const
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+
+std::string mpreal::toString(const std::string& format) const
 {
-	return to_string();
+	char *s = NULL;
+	string out;
+
+	if( !format.empty() )
+	{
+		if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0))
+		{
+			out = std::string(s);
+
+			mpfr_free_str(s);
+		}
+	}
+
+	return out;
 }
 
-std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
+#endif
+
+std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
 {
-	char *s, *ns = NULL;	
+  (void)b;
+  (void)mode;
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+
+	// Use MPFR native function for output
+	char format[128];
+	int digits;
+
+	digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp));
+
+	sprintf(format,"%%.%dRNg",digits);		// Default format
+
+	return toString(std::string(format));
+
+#else
+
+	char *s, *ns = NULL; 
 	size_t slen, nslen;
 	mp_exp_t exp;
 	string out;
-	
+
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
 	set_custom_malloc();
-	
+#endif
+
 	if(mpfr_inf_p(mp))
 	{ 
-		if(mpfr_sgn(mp)>0) return "+@Inf@";
-		else			   return "-@Inf@";
+		if(mpfr_sgn(mp)>0) return "+Inf";
+		else			   return "-Inf";
 	}
 
 	if(mpfr_zero_p(mp)) return "0";
-	if(mpfr_nan_p(mp))  return "@NaN@";
-		
+	if(mpfr_nan_p(mp))  return "NaN";
+
 	s  = mpfr_get_str(NULL,&exp,b,0,mp,mode);
 	ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
 
@@ -419,8 +536,8 @@ std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
 			// Make final string
 			if(--exp)
 			{
-				if(exp>0) out += "e+"+mpfr::to_string<mp_exp_t>(exp,std::dec);
-				else 	  out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec);
+				if(exp>0) out += "e+"+mpfr::toString<mp_exp_t>(exp,std::dec);
+				else 	  out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec);
 			}
 		}
 
@@ -429,79 +546,52 @@ std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
 	}else{
 		return "conversion error!";
 	}
+#endif
 }
 
+
 //////////////////////////////////////////////////////////////////////////
 // I/O
 ostream& operator<<(ostream& os, const mpreal& v)
 {
-	return os<<v.to_string(static_cast<size_t>(os.precision()));
+	return os<<v.toString(static_cast<int>(os.precision()));
 }
 
 istream& operator>>(istream &is, mpreal& v)
 {
-	char c;	
-	string s = "";
-	mpfr_t t;
-	
-	mpreal::set_custom_malloc();
-	
-	if(is.good())
-	{
-		is>>ws;
-		while ((c = is.get())!=EOF)
-		{
-			if(c ==' ' || c == '\t' || c == '\n' || c == '\r')
-			{
-				is.putback(c);
-				break;
-			}
-			s += c;
-		}
-
-		if(s.size() != 0)
-		{
-			// Protect current value from alternation in case of input error
-			// so some error handling(roll back) procedure can be used 			
-
-			if(0==mpfr_init_set_str(t,s.c_str(),mpreal::default_base,mpreal::default_rnd))
-			{
-				mpfr_set(v.mp,t,mpreal::default_rnd);				
-				mpfr_clear(t);
-
-			}else{
-				mpfr_clear(t);
-				cerr<<"error reading from istream"<<endl;
-				// throw an exception
-			}
-		}
-	}
+	string tmp;
+	is >> tmp;
+	mpfr_set_str(v.mp, tmp.c_str(),mpreal::default_base,mpreal::default_rnd);
 	return is;
 }
 
-// Optimized dynamic memory allocation/(re-)deallocation.
-void * mpreal::mpreal_allocate(size_t alloc_size)
-{
-	return(dlmalloc(alloc_size));
-}
 
-void * mpreal::mpreal_reallocate(void *ptr, size_t /*old_size*/, size_t new_size)
-{
-	return(dlrealloc(ptr,new_size));
-}
-
-void mpreal::mpreal_free(void *ptr, size_t /*size*/)
-{
-	dlfree(ptr);
-}
-
-inline void mpreal::set_custom_malloc(void)
-{
-	if(!is_custom_malloc)
+#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
+	// Optimized dynamic memory allocation/(re-)deallocation.
+	void * mpreal::mpreal_allocate(size_t alloc_size)
 	{
-		mp_set_memory_functions(mpreal_allocate,mpreal_reallocate,mpreal_free);
-		is_custom_malloc = true;
+		return(dlmalloc(alloc_size));
 	}
-}
+
+	void * mpreal::mpreal_reallocate(void *ptr, size_t old_size, size_t new_size)
+	{
+		return(dlrealloc(ptr,new_size));
+	}
+
+	void mpreal::mpreal_free(void *ptr, size_t size)
+	{
+		dlfree(ptr);
+	}
+
+	inline void mpreal::set_custom_malloc(void)
+	{
+		if(!is_custom_malloc)
+		{
+			mp_set_memory_functions(mpreal_allocate,mpreal_reallocate,mpreal_free);
+			is_custom_malloc = true;
+		}
+	}
+#endif
+
 }
 

unsupported/test/mpreal_support.cpp

@@ -2,6 +2,7 @@
 #include <Eigen/MPRealSupport>
 #include <Eigen/LU>
 #include <Eigen/Eigenvalues>
+#include <sstream>
 
 using namespace mpfr;
 using namespace std;
@@ -24,6 +25,15 @@ void test_mpreal_support()
     MatrixXmp B = MatrixXmp::Random(s,s);
     MatrixXmp S = A.adjoint() * A;
     MatrixXmp X;
+    
+    // Basic stuffs
+    VERIFY_IS_APPROX(A.real(), A);
+    VERIFY(Eigen::internal::isApprox(A.array().abs2().sum(), A.squaredNorm()));
+    VERIFY_IS_APPROX(A.array().exp(),         exp(A.array()));
+    VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs());
+    VERIFY_IS_APPROX(A.array().sin(),         sin(A.array()));
+    VERIFY_IS_APPROX(A.array().cos(),         cos(A.array()));
+    
 
     // Cholesky
     X = S.selfadjointView<Lower>().llt().solve(B);
@@ -39,6 +49,13 @@ void test_mpreal_support()
     VERIFY_IS_APPROX((S.selfadjointView<Lower>() * eig.eigenvectors()),
                       eig.eigenvectors() * eig.eigenvalues().asDiagonal());
   }
+  
+  {
+    MatrixXmp A(8,3); A.setRandom();
+    // test output (interesting things happen in this code)
+    std::stringstream stream;
+    stream << A;
+  }
 }
 
 extern "C" {