* take advantage of new possibilies in LLT (mat -= product)

* fix Block::operator+= product which was not optimized
* fix some compilation issues

Eigen/src/Cholesky/LLT.h

@@ -26,7 +26,7 @@
 #define EIGEN_LLT_H
 
 template<typename MatrixType, int UpLo> struct LLT_Traits;
-    
+
 /** \ingroup cholesky_Module
   *
   * \class LLT
@@ -70,7 +70,7 @@ template<typename MatrixType, int _UpLo> class LLT
 
   public:
 
-    /** 
+    /**
     * \brief Default Constructor.
     *
     * The default constructor is useful in cases in which the user intends to
@@ -94,7 +94,7 @@ template<typename MatrixType, int _UpLo> class LLT
 
     /** \returns a view of the lower triangular matrix L */
     inline typename Traits::MatrixL matrixL() const
-    { 
+    {
       ei_assert(m_isInitialized && "LLT is not initialized.");
       return Traits::getL(m_matrix);
     }
@@ -123,7 +123,7 @@ bool ei_inplace_llt_lo(MatrixType& mat)
   typedef typename MatrixType::RealScalar RealScalar;
   assert(mat.rows()==mat.cols());
   const int size = mat.rows();
-  Matrix<Scalar,Dynamic,1> aux(size);
+
   // 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
@@ -147,16 +147,10 @@ bool ei_inplace_llt_lo(MatrixType& mat)
     mat.coeffRef(j,j) = x = ei_sqrt(x);
 
     int endSize = size-j-1;
-    if (endSize>0) {
+    if (endSize>0)
-      // Note that when all matrix columns have good alignment, then the following
+    {
-      // product is guaranteed to be optimal with respect to alignment.
+      mat.col(j).end(endSize) -= (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()).lazy();
-      aux.end(endSize) = 
+      mat.col(j).end(endSize) *= RealScalar(1)/x;
-        (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint()).lazy();
-
-      mat.col(j).end(endSize) = (mat.col(j).end(endSize) - aux.end(endSize) ) / x;
-      // TODO improve the products so that the following is efficient:
-//       mat.col(j).end(endSize) -= (mat.block(j+1, 0, endSize, j) * mat.row(j).start(j).adjoint());
-//       mat.col(j).end(endSize) *= Scalar(1)/x;
     }
   }
 
@@ -170,8 +164,7 @@ bool ei_inplace_llt_up(MatrixType& mat)
   typedef typename MatrixType::RealScalar RealScalar;
   assert(mat.rows()==mat.cols());
   const int size = mat.rows();
-  Matrix<Scalar,Dynamic,1> aux(size);
+
-  
   const RealScalar cutoff = machine_epsilon<Scalar>() * size * mat.diagonal().cwise().abs().maxCoeff();
   RealScalar x;
   x = ei_real(mat.coeff(0,0));
@@ -190,10 +183,8 @@ bool ei_inplace_llt_up(MatrixType& mat)
 
     int endSize = size-j-1;
     if (endSize>0) {
-      aux.start(endSize) =
+      mat.row(j).end(endSize) -= (mat.col(j).start(j).adjoint() * mat.block(0, j+1, j, endSize)).lazy();
-        (mat.block(0, j+1, j, endSize).adjoint() * mat.col(j).start(j)).lazy();
+      mat.row(j).end(endSize) *= RealScalar(1)/x;
-
-      mat.row(j).end(endSize) = (mat.row(j).end(endSize) - aux.start(endSize).adjoint() ) / x;
     }
   }
 

Eigen/src/Core/MapBase.h

@@ -172,8 +172,9 @@ template<typename Derived> class MapBase
     }
 
     using Base::operator=;
-
     using Base::operator*=;
+    using Base::operator+=;
+    using Base::operator-=;
 
     template<typename OtherDerived>
     Derived& operator+=(const MatrixBase<OtherDerived>& other)

Eigen/src/Core/SolveTriangular.h

@@ -176,7 +176,8 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
           IsLowerTriangular ? size-endBlock : endBlock+1,
           &(lhs.const_cast_derived().coeffRef(IsLowerTriangular ? endBlock : 0, IsLowerTriangular ? startBlock : endBlock+1)),
           lhs.stride(),
-          btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c)));
+          btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c)),
+          Scalar(1));
 // 				if (IsLowerTriangular)
 //           other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock)
 //                                           * other.col(c).block(startBlock,endBlock-startBlock)).lazy();

Eigen/src/Core/products/GeneralMatrixMatrix.h

@@ -151,7 +151,6 @@ static void ei_cache_friendly_product(
         const Scalar* b3 = &rhs[(j2+3)*rhsStride + k2];
         if (hasAlpha)
         {
-          std::cerr << "* by " << alpha << "\n";
           for(int k=0; k<actual_kc; k++)
           {
             ei_pstore(&blockB[count+0*PacketSize], ei_pset1(alpha*b0[k]));

test/cholesky.cpp

@@ -91,6 +91,14 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     VERIFY_IS_APPROX(symm * vecX, vecB);
     chol.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());
+    cholup.solve(vecB, &vecX);
+    VERIFY_IS_APPROX(symm * vecX, vecB);
+    cholup.solve(matB, &matX);
+    VERIFY_IS_APPROX(symm * matX, matB);
   }
 
   int sign = ei_random<int>()%2 ? 1 : -1;