Google Git
Sign in
eigen/mirror/a8f943127cfd79e4c059e33307d2c60369242f3e^2..a8f943127cfd79e4c059e33307d2c60369242f3e/.
commit
a8f943127cfd79e4c059e33307d2c60369242f3e
[log]
author
Gael Guennebaud <g.gael@free.fr>
Mon Aug 03 16:11:30 2009 +0200
committer
Gael Guennebaud <g.gael@free.fr>
Mon Aug 03 16:11:30 2009 +0200
tree
0865f13e0084b3253df10fa80f02fcfd29254655
parent
3cf5bb31f6b6e3b6b8f229ed1658cb867fe6e8f5 [diff]
parent
d10c710b15676935ccc93faf59fe0f37e0f1b7cd [diff]
merge

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 4527067..ef04b8f 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h

@@ -116,80 +116,81 @@
     bool m_isInitialized;
 };
 
-template<typename MatrixType>
-bool ei_inplace_llt_lo(MatrixType& mat)
+// forward declaration (defined at the end of this file)
+template<int UpLo> struct ei_llt_inplace;
+
+template<> struct ei_llt_inplace<LowerTriangular>
 {
-  typedef typename MatrixType::Scalar Scalar;
-  typedef typename MatrixType::RealScalar RealScalar;
-  assert(mat.rows()==mat.cols());
-  const int size = mat.rows();
-
-  // The biggest overall is the point of reference to which further diagonals
-  // are compared; if any diagonal is negligible compared
-  // to the largest overall, the algorithm bails.  This cutoff is suggested
-  // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
-  // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
-  // Algorithms" page 217, also by Higham.
-  const RealScalar cutoff = machine_epsilon<Scalar>() * size * mat.diagonal().cwise().abs().maxCoeff();
-  RealScalar x;
-  x = ei_real(mat.coeff(0,0));
-  mat.coeffRef(0,0) = ei_sqrt(x);
-  if(size==1)
+  template<typename MatrixType>
+  static bool unblocked(MatrixType& mat)
   {
-    return true;
-  }
-  mat.col(0).end(size-1) = mat.col(0).end(size-1) / ei_real(mat.coeff(0,0));
-  for (int j = 1; j < size; ++j)
-  {
-    x = ei_real(mat.coeff(j,j)) - mat.row(j).start(j).squaredNorm();
-    if (ei_abs(x) < cutoff) continue;
-
-    mat.coeffRef(j,j) = x = ei_sqrt(x);
-
-    int endSize = size-j-1;
-    if (endSize>0)
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    ei_assert(mat.rows()==mat.cols());
+    const int size = mat.rows();
+    for(int k = 0; k < size; ++k)
     {
-      mat.col(j).end(endSize) -= (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()).lazy();
-      mat.col(j).end(endSize) *= RealScalar(1)/x;
+      int rs = size-k-1; // remaining size
+
+      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
+      Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
+      Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
+
+      RealScalar x = ei_real(mat.coeff(k,k));
+      if (k>0) x -= mat.row(k).start(k).squaredNorm();
+      if (x<=RealScalar(0))
+        return false;
+      mat.coeffRef(k,k) = x = ei_sqrt(x);
+      if (k>0 && rs>0) A21 -= (A20 * A10.adjoint()).lazy();
+      if (rs>0) A21 *= RealScalar(1)/x;
     }
-  }
-
-  return true;
-}
-
-template<typename MatrixType>
-bool ei_inplace_llt_up(MatrixType& mat)
-{
-  typedef typename MatrixType::Scalar Scalar;
-  typedef typename MatrixType::RealScalar RealScalar;
-  assert(mat.rows()==mat.cols());
-  const int size = mat.rows();
-
-  const RealScalar cutoff = machine_epsilon<Scalar>() * size * mat.diagonal().cwise().abs().maxCoeff();
-  RealScalar x;
-  x = ei_real(mat.coeff(0,0));
-  mat.coeffRef(0,0) = ei_sqrt(x);
-  if(size==1)
-  {
     return true;
   }
-  mat.row(0).end(size-1) = mat.row(0).end(size-1) / ei_real(mat.coeff(0,0));
-  for (int j = 1; j < size; ++j)
+
+  template<typename MatrixType>
+  static bool blocked(MatrixType& m)
   {
-    x = ei_real(mat.coeff(j,j)) - mat.col(j).start(j).squaredNorm();
-    if (ei_abs(x) < cutoff) continue;
+    ei_assert(m.rows()==m.cols());
+    int size = m.rows();
+    if(size<32)
+      return unblocked(m);
 
-    mat.coeffRef(j,j) = x = ei_sqrt(x);
+    int blockSize = size/8;
+    blockSize = (blockSize/16)*16;
+    blockSize = std::min(std::max(blockSize,8), 128);
 
-    int endSize = size-j-1;
-    if (endSize>0) {
-      mat.row(j).end(endSize) -= (mat.col(j).start(j).adjoint() * mat.block(0, j+1, j, endSize)).lazy();
-      mat.row(j).end(endSize) *= RealScalar(1)/x;
+    for (int k=0; k<size; k+=blockSize)
+    {
+      int bs = std::min(blockSize, size-k);
+      int rs = size - k - bs;
+
+      Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
+      Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
+      Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
+
+      if(!unblocked(A11)) return false;
+      if(rs>0) A11.adjoint().template triangularView<UpperTriangular>().template solveInPlace<OnTheRight>(A21);
+      if(rs>0) A22.template selfadjointView<LowerTriangular>().rankUpdate(A21,-1); // bottleneck
     }
+    return true;
   }
+};
 
-  return true;
-}
+template<> struct ei_llt_inplace<UpperTriangular>
+{
+  template<typename MatrixType>
+  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat)
+  {
+    Transpose<MatrixType> matt(mat);
+    return ei_llt_inplace<LowerTriangular>::unblocked(matt);
+  }
+  template<typename MatrixType>
+  static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat)
+  {
+    Transpose<MatrixType> matt(mat);
+    return ei_llt_inplace<LowerTriangular>::blocked(matt);
+  }
+};
 
 template<typename MatrixType> struct LLT_Traits<MatrixType,LowerTriangular>
 {
@@ -198,7 +199,7 @@
   inline static MatrixL getL(const MatrixType& m) { return m; }
   inline static MatrixU getU(const MatrixType& m) { return m.adjoint().nestByValue(); }
   static bool inplace_decomposition(MatrixType& m)
-  { return ei_inplace_llt_lo(m); }
+  { return ei_llt_inplace<LowerTriangular>::blocked(m); }
 };
 
 template<typename MatrixType> struct LLT_Traits<MatrixType,UpperTriangular>
@@ -208,7 +209,7 @@
   inline static MatrixL getL(const MatrixType& m) { return m.adjoint().nestByValue(); }
   inline static MatrixU getU(const MatrixType& m) { return m; }
   static bool inplace_decomposition(MatrixType& m)
-  { return ei_inplace_llt_up(m); }
+  { return ei_llt_inplace<UpperTriangular>::blocked(m); }
 };
 
 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix

diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h
index 2da463a..c226969 100644
--- a/Eigen/src/Core/BandMatrix.h
+++ b/Eigen/src/Core/BandMatrix.h

@@ -57,7 +57,7 @@
 };
 
 template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options>
-class BandMatrix : public MultiplierBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
+class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
 {
   public:
 

diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 5fc80c9..ebbed15 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h

@@ -27,7 +27,7 @@
 #define EIGEN_DIAGONALMATRIX_H
 
 template<typename Derived>
-class DiagonalBase : public MultiplierBase<Derived>
+class DiagonalBase : public AnyMatrixBase<Derived>
 {
   public:
     typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType;
@@ -52,6 +52,12 @@
     DenseMatrixType toDenseMatrix() const { return derived(); }
     template<typename DenseDerived>
     void evalToDense(MatrixBase<DenseDerived> &other) const;
+    template<typename DenseDerived>
+    void addToDense(MatrixBase<DenseDerived> &other) const
+    { other.diagonal() += diagonal(); }
+    template<typename DenseDerived>
+    void subToDense(MatrixBase<DenseDerived> &other) const
+    { other.diagonal() -= diagonal(); }
 
     inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
     inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
@@ -84,6 +90,7 @@
   */
 template<typename _Scalar, int _Size>
 struct ei_traits<DiagonalMatrix<_Scalar,_Size> >
+ : ei_traits<Matrix<_Scalar,_Size,_Size> >
 {
   typedef Matrix<_Scalar,_Size,1> DiagonalVectorType;
 };
@@ -170,6 +177,14 @@
 struct ei_traits<DiagonalWrapper<_DiagonalVectorType> >
 {
   typedef _DiagonalVectorType DiagonalVectorType;
+  typedef typename DiagonalVectorType::Scalar Scalar;
+  enum {
+    RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    MaxRowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    MaxColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    Flags = 0
+  };
 };
 
 template<typename _DiagonalVectorType>

diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index 8937596..c31acab 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h

@@ -462,6 +462,7 @@
       : m_storage(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols())
     {
       _check_template_params();
+      resize(other.rows(), other.cols());
       *this = other;
     }
 

diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index eb29c0e..f23188e 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h

@@ -37,20 +37,31 @@
   */
 template<typename Derived> struct AnyMatrixBase
 {
+  typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
+  
   Derived& derived() { return *static_cast<Derived*>(this); }
   const Derived& derived() const { return *static_cast<const Derived*>(this); }
-};
+  /** \returns the number of rows. \sa cols(), RowsAtCompileTime */
+  inline int rows() const { return derived().rows(); }
+  /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
+  inline int cols() const { return derived().cols(); }
+  
+  template<typename Dest> inline void evalTo(Dest& dst) const
+  { derived().evalTo(dst); }
 
-/** Common base class for all classes T such that there are overloaded operator* allowing to
-  * multiply a MatrixBase by a T on both sides.
-  *
-  * In other words, an AnyMatrixBase object is an object that can be multiplied a MatrixBase, the result being again a MatrixBase.
-  *
-  * Besides MatrixBase-derived classes, this also includes certain special matrix classes, such as diagonal matrices.
-  */
-template<typename Derived> struct MultiplierBase : public AnyMatrixBase<Derived>
-{
-  using AnyMatrixBase<Derived>::derived;
+  template<typename Dest> inline void addToDense(Dest& dst) const
+  {
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalToDense(res);
+    dst += res;
+  }
+  
+  template<typename Dest> inline void subToDense(Dest& dst) const
+  {
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalToDense(res);
+    dst -= res;
+  }
 };
 
 /** \class MatrixBase
@@ -80,7 +91,7 @@
   */
 template<typename Derived> class MatrixBase
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-  : public MultiplierBase<Derived>,
+  : public AnyMatrixBase<Derived>,
     public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
                 typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>
 #endif // not EIGEN_PARSED_BY_DOXYGEN
@@ -299,12 +310,16 @@
     Derived& operator=(const AnyMatrixBase<OtherDerived> &other)
     { other.derived().evalToDense(derived()); return derived(); }
 
+    template<typename OtherDerived>
+    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other)
+    { other.derived().addToDense(derived()); return derived(); }
+
+    template<typename OtherDerived>
+    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other)
+    { other.derived().subToDense(derived()); return derived(); }
+
     template<typename OtherDerived,typename OtherEvalType>
     Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
-    template<typename OtherDerived,typename OtherEvalType>
-    Derived& operator+=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
-    template<typename OtherDerived,typename OtherEvalType>
-    Derived& operator-=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** Copies \a other into *this without evaluating other. \returns a reference to *this. */
@@ -411,7 +426,7 @@
     operator*(const MatrixBase<OtherDerived> &other) const;
 
     template<typename OtherDerived>
-    Derived& operator*=(const MultiplierBase<OtherDerived>& other);
+    Derived& operator*=(const AnyMatrixBase<OtherDerived>& other);
 
     template<typename DiagonalDerived>
     const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
@@ -646,7 +661,7 @@
     void visit(Visitor& func) const;
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-    using MultiplierBase<Derived>::derived;
+    using AnyMatrixBase<Derived>::derived;
     inline Derived& const_cast_derived() const
     { return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
 #endif // not EIGEN_PARSED_BY_DOXYGEN

diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index ff45cba..1a32eb5 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h

@@ -294,7 +294,7 @@
 template<typename Derived>
 template<typename OtherDerived>
 inline Derived &
-MatrixBase<Derived>::operator*=(const MultiplierBase<OtherDerived> &other)
+MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other)
 {
   return derived() = derived() * other.derived();
 }

diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h
index 58b205e..3f2b478 100644
--- a/Eigen/src/Core/ReturnByValue.h
+++ b/Eigen/src/Core/ReturnByValue.h

@@ -48,14 +48,6 @@
   public:
     template<typename Dest> inline void evalTo(Dest& dst) const
     { static_cast<const Functor*>(this)->evalTo(dst); }
-    template<typename Dest> inline void addTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_addTo(dst); }
-    template<typename Dest> inline void subTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_subTo(dst); }
-    template<typename Dest> inline void _addTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst += res; }
-    template<typename Dest> inline void _subTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst -= res; }
 };
 
 template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int _MaxRows,int _MaxCols>
@@ -68,14 +60,6 @@
     template<typename Dest>
     inline void evalTo(Dest& dst) const
     { static_cast<const Functor* const>(this)->evalTo(dst); }
-    template<typename Dest> inline void addTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_addTo(dst); }
-    template<typename Dest> inline void subTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_subTo(dst); }
-    template<typename Dest> inline void _addTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst += res; }
-    template<typename Dest> inline void _subTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst -= res; }
     inline int rows() const { return static_cast<const Functor* const>(this)->rows(); }
     inline int cols() const { return static_cast<const Functor* const>(this)->cols(); }
 };
@@ -88,20 +72,4 @@
   return derived();
 }
 
-template<typename Derived>
-template<typename OtherDerived,typename OtherEvalType>
-Derived& MatrixBase<Derived>::operator+=(const ReturnByValue<OtherDerived,OtherEvalType>& other)
-{
-  other.addTo(derived());
-  return derived();
-}
-
-template<typename Derived>
-template<typename OtherDerived,typename OtherEvalType>
-Derived& MatrixBase<Derived>::operator-=(const ReturnByValue<OtherDerived,OtherEvalType>& other)
-{
-  other.subTo(derived());
-  return derived();
-}
-
 #endif // EIGEN_RETURNBYVALUE_H

diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index c21f3a3..0a4ba17 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h

@@ -55,7 +55,7 @@
 
 template <typename Lhs, int LhsMode, bool LhsIsVector,
           typename Rhs, int RhsMode, bool RhsIsVector>
-struct ei_selfadjoint_product_returntype;
+struct SelfadjointProductMatrix;
 
 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
@@ -100,20 +100,20 @@
 
     /** Efficient self-adjoint matrix times vector/matrix product */
     template<typename OtherDerived>
-    ei_selfadjoint_product_returntype<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
+    SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
     operator*(const MatrixBase<OtherDerived>& rhs) const
     {
-      return ei_selfadjoint_product_returntype
+      return SelfadjointProductMatrix
               <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
               (m_matrix, rhs.derived());
     }
 
     /** Efficient vector/matrix times self-adjoint matrix product */
     template<typename OtherDerived> friend
-    ei_selfadjoint_product_returntype<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
+    SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
     operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
     {
-      return ei_selfadjoint_product_returntype
+      return SelfadjointProductMatrix
               <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
               (lhs.derived(),rhs.m_matrix);
     }
@@ -201,10 +201,13 @@
 ***************************************************************************/
 
 template<typename Lhs, int LhsMode, typename Rhs>
-struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>
-  : public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<typename Lhs, int LhsMode, typename Rhs>
+struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
+  : public AnyMatrixBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> >
 {
   typedef typename Lhs::Scalar Scalar;
 
@@ -224,19 +227,19 @@
     LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit)
   };
 
-  ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
   inline int rows() const { return m_lhs.rows(); }
   inline int cols() const { return m_rhs.cols(); }
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.setZero();
     evalTo(dst,1);
@@ -272,12 +275,15 @@
 ***************************************************************************/
 
 template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
-struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>
-  : public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
+struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>
+  : public AnyMatrixBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> >
 {
-  ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
@@ -305,12 +311,12 @@
     RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit
   };
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.setZero();
     evalTo(dst,1);

diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h
index 9b67dd5..810b082 100644
--- a/Eigen/src/Core/SolveTriangular.h
+++ b/Eigen/src/Core/SolveTriangular.h

@@ -158,7 +158,7 @@
     const ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
     ei_triangular_solve_matrix<Scalar,Side,Mode,LhsProductTraits::NeedToConjugate,StorageOrder,
                                Rhs::Flags&RowMajorBit>
-      ::run(lhs.rows(), rhs.cols(), &actualLhs.coeff(0,0), actualLhs.stride(), &rhs.coeffRef(0,0), rhs.stride());
+      ::run(lhs.rows(), Side==OnTheLeft? rhs.cols() : rhs.rows(), &actualLhs.coeff(0,0), actualLhs.stride(), &rhs.coeffRef(0,0), rhs.stride());
   }
 };
 
@@ -215,25 +215,25 @@
   * See TriangularView:solve() for the details.
   */
 template<typename MatrixType, unsigned int Mode>
-template<int Side, typename RhsDerived>
-void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<RhsDerived>& _rhs) const
+template<int Side, typename OtherDerived>
+void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<OtherDerived>& _other) const
 {
-  RhsDerived& rhs = _rhs.const_cast_derived();
+  OtherDerived& other = _other.const_cast_derived();
   ei_assert(cols() == rows());
-  ei_assert(cols() == rhs.rows());
+  ei_assert( (Side==OnTheLeft && cols() == other.rows()) || (Side==OnTheRight && cols() == other.cols()) );
   ei_assert(!(Mode & ZeroDiagBit));
   ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit));
 
-  enum { copy = ei_traits<RhsDerived>::Flags & RowMajorBit  && RhsDerived::IsVectorAtCompileTime };
+  enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit  && OtherDerived::IsVectorAtCompileTime };
   typedef typename ei_meta_if<copy,
-    typename ei_plain_matrix_type_column_major<RhsDerived>::type, RhsDerived&>::ret RhsCopy;
-  RhsCopy rhsCopy(rhs);
+    typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy;
+  OtherCopy otherCopy(other);
 
-  ei_triangular_solver_selector<MatrixType, typename ei_unref<RhsCopy>::type,
-    Side, Mode>::run(_expression(), rhsCopy);
+  ei_triangular_solver_selector<MatrixType, typename ei_unref<OtherCopy>::type,
+    Side, Mode>::run(_expression(), otherCopy);
 
   if (copy)
-    rhs = rhsCopy;
+    other = otherCopy;
 }
 
 /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular.

diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h
index 77865c4..0787daa 100644
--- a/Eigen/src/Core/Transpose.h
+++ b/Eigen/src/Core/Transpose.h

@@ -70,7 +70,9 @@
     inline int rows() const { return m_matrix.cols(); }
     inline int cols() const { return m_matrix.rows(); }
     inline int nonZeros() const { return m_matrix.nonZeros(); }
-    inline int stride(void) const { return m_matrix.stride(); }
+    inline int stride() const { return m_matrix.stride(); }
+    inline Scalar* data() { return m_matrix.data(); }
+    inline const Scalar* data() const { return m_matrix.data(); }
 
     inline Scalar& coeffRef(int row, int col)
     {

diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 8b6c9a2..c262ea7 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h

@@ -43,7 +43,7 @@
   *
   * \sa MatrixBase::part()
   */
-template<typename Derived> class TriangularBase : public MultiplierBase<Derived>
+template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
 {
   public:
 
@@ -145,7 +145,7 @@
 template<int Mode, bool LhsIsTriangular,
          typename Lhs, bool LhsIsVector,
          typename Rhs, bool RhsIsVector>
-struct ei_triangular_product_returntype;
+struct TriangularProduct;
 
 template<typename _MatrixType, unsigned int _Mode> class TriangularView
   : public TriangularBase<TriangularView<_MatrixType, _Mode> >
@@ -253,20 +253,20 @@
 
     /** Efficient triangular matrix times vector/matrix product */
     template<typename OtherDerived>
-    ei_triangular_product_returntype<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
+    TriangularProduct<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
     operator*(const MatrixBase<OtherDerived>& rhs) const
     {
-      return ei_triangular_product_returntype
+      return TriangularProduct
               <Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
               (m_matrix, rhs.derived());
     }
 
     /** Efficient vector/matrix times triangular matrix product */
     template<typename OtherDerived> friend
-    ei_triangular_product_returntype<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
+    TriangularProduct<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
     operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
     {
-      return ei_triangular_product_returntype
+      return TriangularProduct
               <Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
               (lhs.derived(),rhs.m_matrix);
     }

diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h
index ce18941..f69c043 100644
--- a/Eigen/src/Core/products/TriangularMatrixMatrix.h
+++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h

@@ -321,12 +321,15 @@
 ***************************************************************************/
 
 template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
-struct ei_triangular_product_returntype<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
-  : public ReturnByValue<ei_triangular_product_returntype<Mode,LhsIsTriangular,Lhs,false,Rhs,false>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
+struct TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
+  : public AnyMatrixBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false> >
 {
-  ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  TriangularProduct(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
@@ -347,12 +350,12 @@
   typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
   typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.resize(m_lhs.rows(), m_rhs.cols());
     dst.setZero();

diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h
index 18d76b9..42239fa 100644
--- a/Eigen/src/Core/products/TriangularMatrixVector.h
+++ b/Eigen/src/Core/products/TriangularMatrixVector.h

@@ -118,10 +118,13 @@
 ***************************************************************************/
 
 template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs>
-struct ei_triangular_product_returntype<Mode,true,Lhs,false,Rhs,true>
-  : public ReturnByValue<ei_triangular_product_returntype<Mode,true,Lhs,false,Rhs,true>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<TriangularProduct<Mode,true,Lhs,false,Rhs,true> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs>
+struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
+  : public AnyMatrixBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true> >
 {
   typedef typename Lhs::Scalar Scalar;
 
@@ -137,19 +140,19 @@
   typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
   typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
 
-  ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  TriangularProduct(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
   inline int rows() const { return m_lhs.rows(); }
   inline int cols() const { return m_rhs.cols(); }
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.setZero();
     evalTo(dst,1);

diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index b4fbae2..310d0fb 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h

@@ -29,7 +29,6 @@
 template<typename T> struct NumTraits;
 
 template<typename Derived> struct AnyMatrixBase;
-template<typename Derived> struct MultiplierBase;
 
 template<typename _Scalar, int _Rows, int _Cols,
          int _Options = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION | AutoAlign,

diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h
index b3a40f0..337438e 100644
--- a/Eigen/src/LU/PartialLU.h
+++ b/Eigen/src/LU/PartialLU.h

@@ -201,6 +201,100 @@
   compute(matrix);
 }
 
+/** \internal performs the LU decomposition in place of the matrix \a lu.
+  * In addition, this function returns the row transpositions in the
+  * vector \a row_transpositions which must have a size equal to the number
+  * of columns of the matrix \a lu, and an integer \a nb_transpositions
+  * which returns the actual number of transpositions.
+  */
+template<typename MatrixType, typename IntVector>
+void ei_lu_unblocked(MatrixType& lu, IntVector& row_transpositions, int& nb_transpositions)
+{
+  const int rows = lu.rows();
+  const int size = std::min(lu.rows(),lu.cols());
+  nb_transpositions = 0;
+  for(int k = 0; k < size; ++k)
+  {
+    int row_of_biggest_in_col;
+    lu.block(k,k,rows-k,1).cwise().abs().maxCoeff(&row_of_biggest_in_col);
+    row_of_biggest_in_col += k;
+
+    row_transpositions.coeffRef(k) = row_of_biggest_in_col;
+
+    if(k != row_of_biggest_in_col)
+    {
+      lu.row(k).swap(lu.row(row_of_biggest_in_col));
+      ++nb_transpositions;
+    }
+
+    if(k<rows-1)
+    {
+      lu.col(k).end(rows-k-1) /= lu.coeff(k,k);
+      for(int col = k + 1; col < size; ++col)
+        lu.col(col).end(rows-k-1) -= lu.col(k).end(rows-k-1) * lu.coeff(k,col);
+    }
+  }
+}
+
+/** This is the blocked version of ei_lu_unblocked() */
+template<typename MatrixType, typename IntVector>
+void ei_lu_blocked(MatrixType& lu, IntVector& row_transpositions, int& nb_transpositions)
+{
+  const int size = lu.rows();
+
+  // automatically adjust the number of subdivisions to the size
+  // of the matrix so that there is enough sub blocks:
+  int blockSize = size/8;
+  blockSize = (blockSize/16)*16;
+  blockSize = std::min(std::max(blockSize,8), 256);
+  // if the matrix is too small, no blocking:
+  if(size<32)
+    blockSize = size;
+
+  nb_transpositions = 0;
+  for(int k = 0; k < size; k+=blockSize)
+  {
+    int bs = std::min(size-k,blockSize);
+    int ps = size - k;
+    int rs = size - k - bs;
+    // partition the matrix:
+    //        A00 | A01 | A02
+    // lu  =  A10 | A11 | A12
+    //        A20 | A21 | A22
+    Block<MatrixType,Dynamic,Dynamic> A_0(lu,0,0,size,k);
+    Block<MatrixType,Dynamic,Dynamic> A11_21(lu,k,k,ps,bs);
+    Block<MatrixType,Dynamic,Dynamic> A_2(lu,0,k+bs,size,rs);
+    Block<MatrixType,Dynamic,Dynamic> A11(lu,k,k,bs,bs);
+    Block<MatrixType,Dynamic,Dynamic> A12(lu,k,k+bs,bs,rs);
+    Block<MatrixType,Dynamic,Dynamic> A21(lu,k+bs,k,rs,bs);
+    Block<MatrixType,Dynamic,Dynamic> A22(lu,k+bs,k+bs,rs,rs);
+    
+    VectorBlock<IntVector,Dynamic> row_transpositions_in_panel(row_transpositions,k,bs);
+    int nb_transpositions_in_panel;
+    ei_lu_unblocked(A11_21, row_transpositions_in_panel, nb_transpositions_in_panel);
+    nb_transpositions_in_panel += nb_transpositions_in_panel;
+
+    // update permutations and apply them to A10
+    for(int i=k;i<k+bs; ++i)
+    {
+      int piv = (row_transpositions.coeffRef(i) += k);
+      A_0.row(i).swap(A_0.row(piv));
+    }
+
+    if(rs)
+    {
+      // apply permutations to A_2
+      for(int i=k;i<k+bs; ++i)
+        A_2.row(i).swap(A_2.row(row_transpositions.coeff(i)));
+
+      // A12 = A11^-1 A12
+      A11.template triangularView<UnitLowerTriangular>().solveInPlace(A12);
+
+      A22 -= A21 * A12;
+    }
+  }
+}
+
 template<typename MatrixType>
 void PartialLU<MatrixType>::compute(const MatrixType& matrix)
 {
@@ -211,40 +305,15 @@
   const int size = matrix.rows();
 
   IntColVectorType rows_transpositions(size);
-  int number_of_transpositions = 0;
 
-  for(int k = 0; k < size; ++k)
-  {
-    int row_of_biggest_in_col;
-    m_lu.block(k,k,size-k,1).cwise().abs().maxCoeff(&row_of_biggest_in_col);
-    row_of_biggest_in_col += k;
-
-    rows_transpositions.coeffRef(k) = row_of_biggest_in_col;
-
-    if(k != row_of_biggest_in_col) {
-      m_lu.row(k).swap(m_lu.row(row_of_biggest_in_col));
-      ++number_of_transpositions;
-    }
-
-    if(k<size-1) {
-      m_lu.col(k).end(size-k-1) /= m_lu.coeff(k,k);
-      /* I know it's tempting to replace this for loop by a single matrix product. But actually there's no reason why it
-       * should be faster because it's just an exterior vector product; and in practice this gives much slower code with
-       * GCC 4.2-4.4 (this is weird, would be interesting to investigate). On the other hand, it would be worth having a variant
-       * for row-major matrices, traversing in the other direction for better performance, with a meta selector to compile only
-       * one path
-       */
-      for(int col = k + 1; col < size; ++col)
-        m_lu.col(col).end(size-k-1) -= m_lu.col(k).end(size-k-1) * m_lu.coeff(k,col);
-    }
-  }
+  int nb_transpositions;
+  ei_lu_blocked(m_lu, rows_transpositions, nb_transpositions);
+  m_det_p = (nb_transpositions%2) ? -1 : 1;
 
   for(int k = 0; k < size; ++k) m_p.coeffRef(k) = k;
   for(int k = size-1; k >= 0; --k)
     std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
 
-  m_det_p = (number_of_transpositions%2) ? -1 : 1;
-
   m_isInitialized = true;
 }
 

diff --git a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
index e08c4ae..c8b6988 100644
--- a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
+++ b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh

@@ -237,6 +237,15 @@
     spotrf_(&uplo, &N, C, &N, &info);
   }
 
+  static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
+    int N2 = N*N;
+    scopy_(&N2, X, &intone, C, &intone);
+    char uplo = 'L';
+    int info = 0;
+    int * ipiv = (int*)alloca(sizeof(int)*N);
+    sgetrf_(&N, &N, C, &N, ipiv, &info);
+  }
+
   #ifdef HAS_LAPACK
 
   static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
@@ -249,14 +258,7 @@
     sgetc2_(&N, C, &N, ipiv, jpiv, &info);
   }
 
-  static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
-    int N2 = N*N;
-    scopy_(&N2, X, &intone, C, &intone);
-    char uplo = 'L';
-    int info = 0;
-    int * ipiv = (int*)alloca(sizeof(int)*N);
-    sgetrf_(&N, &N, C, &N, ipiv, &info);
-  }
+  
 
   static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){
 #ifdef PUREBLAS

diff --git a/bench/btl/libs/C_BLAS/main.cpp b/bench/btl/libs/C_BLAS/main.cpp
index 5c2f8a3..bb3022d 100644
--- a/bench/btl/libs/C_BLAS/main.cpp
+++ b/bench/btl/libs/C_BLAS/main.cpp

@@ -24,7 +24,7 @@
 
 #include "action_cholesky.hh"
 #include "action_lu_decomp.hh"
-// #include "action_partial_lu_decomp.hh"
+#include "action_partial_lu.hh"
 #include "action_trisolve_matrix.hh"
 
 #ifdef HAS_LAPACK
@@ -52,10 +52,10 @@
   bench<Action_trisolve_matrix<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
 
   bench<Action_cholesky<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
+  bench<Action_partial_lu<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
 
   #ifdef HAS_LAPACK
   bench<Action_lu_decomp<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
-//   bench<Action_partial_lu_decomp<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_hessenberg<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_tridiagonalization<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   #endif

diff --git a/bench/btl/libs/eigen2/eigen2_interface.hh b/bench/btl/libs/eigen2/eigen2_interface.hh
index ffc800a..7a01472 100644
--- a/bench/btl/libs/eigen2/eigen2_interface.hh
+++ b/bench/btl/libs/eigen2/eigen2_interface.hh

@@ -190,10 +190,10 @@
   }
 
   static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){
-//     X = L.template triangularView<LowerTriangular>().solve(B);
-    X = B;
-    ei_triangular_solve_matrix<real,ColMajor,ColMajor,LowerTriangular>
-      ::run(L.cols(), X.cols(), L.data(), L.stride(), X.data(), X.stride());
+    X = L.template triangularView<LowerTriangular>().solve(B);
+//     
+//     ei_triangular_solve_matrix<real,ColMajor,ColMajor,LowerTriangular>
+//       ::run(L.cols(), X.cols(), L.data(), L.stride(), X.data(), X.stride());
   }
 
   static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){

diff --git a/bench/btl/libs/eigen2/main_adv.cpp b/bench/btl/libs/eigen2/main_adv.cpp
index c3e5fac..d98c0cd 100644
--- a/bench/btl/libs/eigen2/main_adv.cpp
+++ b/bench/btl/libs/eigen2/main_adv.cpp

@@ -23,7 +23,7 @@
 #include "action_cholesky.hh"
 #include "action_hessenberg.hh"
 #include "action_lu_decomp.hh"
-// #include "action_partial_lu_decomp.hh"
+#include "action_partial_lu.hh"
 
 BTL_MAIN;
 
@@ -33,7 +33,7 @@
   bench<Action_trisolve_matrix<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_cholesky<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_lu_decomp<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
-//   bench<Action_partial_lu_decomp<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
+  bench<Action_partial_lu<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
 
   bench<Action_hessenberg<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_tridiagonalization<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);

diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 4b3952a..df937fd 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp

@@ -49,52 +49,58 @@
   MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
   SquareMatrixType symm =  a0 * a0.adjoint();
   // let's make sure the matrix is not singular or near singular
-  MatrixType a1 = MatrixType::Random(rows,cols);
-  symm += a1 * a1.adjoint();
+  for (int k=0; k<3; ++k)
+  {
+    MatrixType a1 = MatrixType::Random(rows,cols);
+    symm += a1 * a1.adjoint();
+  }
+
+  SquareMatrixType symmUp = symm.template triangularView<UpperTriangular>();
+  SquareMatrixType symmLo = symm.template triangularView<LowerTriangular>();
 
   // to test if really Cholesky only uses the upper triangular part, uncomment the following
   // FIXME: currently that fails !!
   //symm.template part<StrictlyLowerTriangular>().setZero();
 
   #ifdef HAS_GSL
-  if (ei_is_same_type<RealScalar,double>::ret)
-  {
-    typedef GslTraits<Scalar> Gsl;
-    typename Gsl::Matrix gMatA=0, gSymm=0;
-    typename Gsl::Vector gVecB=0, gVecX=0;
-    convert<MatrixType>(symm, gSymm);
-    convert<MatrixType>(symm, gMatA);
-    convert<VectorType>(vecB, gVecB);
-    convert<VectorType>(vecB, gVecX);
-    Gsl::cholesky(gMatA);
-    Gsl::cholesky_solve(gMatA, gVecB, gVecX);
-    VectorType vecX(rows), _vecX, _vecB;
-    convert(gVecX, _vecX);
-    symm.llt().solve(vecB, &vecX);
-    Gsl::prod(gSymm, gVecX, gVecB);
-    convert(gVecB, _vecB);
-    // test gsl itself !
-    VERIFY_IS_APPROX(vecB, _vecB);
-    VERIFY_IS_APPROX(vecX, _vecX);
-
-    Gsl::free(gMatA);
-    Gsl::free(gSymm);
-    Gsl::free(gVecB);
-    Gsl::free(gVecX);
-  }
+//   if (ei_is_same_type<RealScalar,double>::ret)
+//   {
+//     typedef GslTraits<Scalar> Gsl;
+//     typename Gsl::Matrix gMatA=0, gSymm=0;
+//     typename Gsl::Vector gVecB=0, gVecX=0;
+//     convert<MatrixType>(symm, gSymm);
+//     convert<MatrixType>(symm, gMatA);
+//     convert<VectorType>(vecB, gVecB);
+//     convert<VectorType>(vecB, gVecX);
+//     Gsl::cholesky(gMatA);
+//     Gsl::cholesky_solve(gMatA, gVecB, gVecX);
+//     VectorType vecX(rows), _vecX, _vecB;
+//     convert(gVecX, _vecX);
+//     symm.llt().solve(vecB, &vecX);
+//     Gsl::prod(gSymm, gVecX, gVecB);
+//     convert(gVecB, _vecB);
+//     // test gsl itself !
+//     VERIFY_IS_APPROX(vecB, _vecB);
+//     VERIFY_IS_APPROX(vecX, _vecX);
+// 
+//     Gsl::free(gMatA);
+//     Gsl::free(gSymm);
+//     Gsl::free(gVecB);
+//     Gsl::free(gVecX);
+//   }
   #endif
 
   {
-    LLT<SquareMatrixType> chol(symm);
-    VERIFY_IS_APPROX(symm, chol.matrixL().toDense() * chol.matrixL().adjoint().toDense());
-    chol.solve(vecB, &vecX);
+    LLT<SquareMatrixType,LowerTriangular> chollo(symmLo);
+    VERIFY_IS_APPROX(symm, chollo.matrixL().toDense() * chollo.matrixL().adjoint().toDense());
+    chollo.solve(vecB, &vecX);
     VERIFY_IS_APPROX(symm * vecX, vecB);
-    chol.solve(matB, &matX);
+    chollo.solve(matB, &matX);
     VERIFY_IS_APPROX(symm * matX, matB);
 
     // test the upper mode
-    LLT<SquareMatrixType,UpperTriangular> cholup(symm);
-    VERIFY_IS_APPROX(symm, cholup.matrixL().toDense() * chol.matrixL().adjoint().toDense());
+    LLT<SquareMatrixType,UpperTriangular> cholup(symmUp);
+    VERIFY_IS_APPROX(symm, cholup.matrixL().toDense() * cholup.matrixL().adjoint().toDense());
     cholup.solve(vecB, &vecX);
     VERIFY_IS_APPROX(symm * vecX, vecB);
     cholup.solve(matB, &matX);
@@ -143,17 +149,16 @@
 {
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST( cholesky(Matrix<double,1,1>()) );
-    CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
-    CALL_SUBTEST( cholesky(Matrix2d()) );
-    CALL_SUBTEST( cholesky(Matrix3f()) );
-    CALL_SUBTEST( cholesky(Matrix4d()) );
-    CALL_SUBTEST( cholesky(MatrixXcd(7,7)) );
-    CALL_SUBTEST( cholesky(MatrixXd(17,17)) );
-    CALL_SUBTEST( cholesky(MatrixXf(200,200)) );
+//     CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
+//     CALL_SUBTEST( cholesky(Matrix2d()) );
+//     CALL_SUBTEST( cholesky(Matrix3f()) );
+//     CALL_SUBTEST( cholesky(Matrix4d()) );
+    CALL_SUBTEST( cholesky(MatrixXd(200,200)) );
+    CALL_SUBTEST( cholesky(MatrixXcd(100,100)) );
   }
 
-  CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
-  CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
-  CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
-  CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
+//   CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
+//   CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
+//   CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
+//   CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
 }

diff --git a/test/inverse.cpp b/test/inverse.cpp
index 352887d..b4eef73 100644
--- a/test/inverse.cpp
+++ b/test/inverse.cpp

@@ -86,8 +86,8 @@
     CALL_SUBTEST( inverse(Matrix2d()) );
     CALL_SUBTEST( inverse(Matrix3f()) );
     CALL_SUBTEST( inverse(Matrix4f()) );
-    CALL_SUBTEST( inverse(MatrixXf(8,8)) );
-    CALL_SUBTEST( inverse(MatrixXcd(7,7)) );
+    CALL_SUBTEST( inverse(MatrixXf(72,72)) );
+    CALL_SUBTEST( inverse(MatrixXcd(56,56)) );
   }
 
   // test some tricky cases for 4x4 matrices