From 3ac01b1400dccefcb2db85ab1caf88fc3bd379db Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Thu, 6 Aug 2009 16:41:54 +0200
Subject: [PATCH] compilation fix in EigenSolver, bugfix in PartialLU

---
 Eigen/src/Core/Product.h   |  7 +++++++
 Eigen/src/LU/PartialLU.h   | 31 +++++++++++++++----------------
 Eigen/src/QR/EigenSolver.h | 14 +++++++-------
 test/CMakeLists.txt        |  2 +-
 4 files changed, 30 insertions(+), 24 deletions(-)

diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index 29a0f9cac..18f14f75a 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -132,6 +132,13 @@ struct ProductReturnType<Lhs,Rhs,UnrolledProduct>
 *  Implementation of Inner Vector Vector Product
 ***********************************************************************/
 
+// FIXME : maybe the "inner product" could return a Scalar
+// instead of a 1x1 matrix ?? 
+// Pro: more natural for the user
+// Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix
+// product ends up to a row-vector times col-vector product... To tackle this use
+// case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x);
+
 template<typename Lhs, typename Rhs>
 struct ei_traits<GeneralProduct<Lhs,Rhs,InnerProduct> >
  : ei_traits<ProductBase<GeneralProduct<Lhs,Rhs,InnerProduct>, Lhs, Rhs> >
diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h
index 1cca1c7db..10d39055a 100644
--- a/Eigen/src/LU/PartialLU.h
+++ b/Eigen/src/LU/PartialLU.h
@@ -245,11 +245,10 @@ struct ei_partial_lu_impl
 
       if(k<rows-1)
       {
-        lu.col(k).end(rows-k-1) /= lu.coeff(k,k);
-
-        // TODO implement a fast rank one update routine
-        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);
+        int rrows = rows-k-1;
+        int rsize = size-k-1;
+        lu.col(k).end(rrows) /= lu.coeff(k,k);
+        lu.corner(BottomRight,rrows,rsize) -= (lu.col(k).end(rrows) * lu.row(k).end(rsize)).lazy();
       }
     }
   }
@@ -271,7 +270,6 @@ struct ei_partial_lu_impl
   {
     MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
     MatrixType lu(lu1,0,0,rows,cols);
-    
 
     const int size = std::min(rows,cols);
 
@@ -294,24 +292,25 @@ struct ei_partial_lu_impl
     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;
+      int bs = std::min(size-k,blockSize); // actual size of the block
+      int trows = rows - k - bs; // trailing rows
+      int tsize = size - k - bs; // trailing size
+      
       // partition the matrix:
       //        A00 | A01 | A02
       // lu  =  A10 | A11 | A12
       //        A20 | A21 | A22
-      BlockType A_0(lu,0,0,size,k);
-      BlockType A_2(lu,0,k+bs,size,rs);
+      BlockType A_0(lu,0,0,rows,k);
+      BlockType A_2(lu,0,k+bs,rows,tsize);
       BlockType A11(lu,k,k,bs,bs);
-      BlockType A12(lu,k,k+bs,bs,rs);
-      BlockType A21(lu,k+bs,k,rs,bs);
-      BlockType A22(lu,k+bs,k+bs,rs,rs);
+      BlockType A12(lu,k,k+bs,bs,tsize);
+      BlockType A21(lu,k+bs,k,trows,bs);
+      BlockType A22(lu,k+bs,k+bs,trows,tsize);
 
       int nb_transpositions_in_panel;
       // recursively calls the blocked LU algorithm with a very small
       // blocking size:
-      blocked_lu(ps, bs, &lu.coeffRef(k,k), luStride,
+      blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
                  row_transpositions+k, nb_transpositions_in_panel, 16);
       nb_transpositions_in_panel += nb_transpositions_in_panel;
 
@@ -322,7 +321,7 @@ struct ei_partial_lu_impl
         A_0.row(i).swap(A_0.row(piv));
       }
 
-      if(rs)
+      if(trows)
       {
         // apply permutations to A_2
         for(int i=k;i<k+bs; ++i)
diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h
index 70f21cebc..dbe1aca89 100644
--- a/Eigen/src/QR/EigenSolver.h
+++ b/Eigen/src/QR/EigenSolver.h
@@ -244,9 +244,9 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
       // H = (I-u*u'/h)*H*(I-u*u')/h)
       int bSize = high-m+1;
       matH.block(m, m, bSize, n-m) -= ((ort.segment(m, bSize)/h)
-        * (ort.segment(m, bSize).transpose() *  matH.block(m, m, bSize, n-m)).lazy()).lazy();
+        * (ort.segment(m, bSize).transpose() *  matH.block(m, m, bSize, n-m))).lazy();
 
-      matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)).lazy()
+      matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize))
         * (ort.segment(m, bSize)/h).transpose()).lazy();
 
       ort.coeffRef(m) = scale*ort.coeff(m);
@@ -264,8 +264,8 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
       ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m);
 
       int bSize = high-m+1;
-      m_eivec.block(m, m, bSize, bSize) += ( (ort.segment(m, bSize) /  (matH.coeff(m,m-1) * ort.coeff(m) ) )
-        * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)).lazy());
+      m_eivec.block(m, m, bSize, bSize) += ( (ort.segment(m, bSize) /  (matH.coeff(m,m-1) * ort.coeff(m))) 
+        * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ).lazy();
     }
   }
 }
@@ -585,7 +585,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
       for (int i = n-1; i >= 0; i--)
       {
         w = matH.coeff(i,i) - p;
-        r = (matH.row(i).segment(l,n-l+1) * matH.col(n).segment(l, n-l+1))(0,0);
+        r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1));
 
         if (m_eivalues.coeff(i).imag() < 0.0)
         {
@@ -644,8 +644,8 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
       for (int i = n-2; i >= 0; i--)
       {
         Scalar ra,sa,vr,vi;
-        ra = (matH.block(i,l, 1, n-l+1) * matH.block(l,n-1, n-l+1, 1)).lazy()(0,0);
-        sa = (matH.block(i,l, 1, n-l+1) * matH.block(l,n, n-l+1, 1)).lazy()(0,0);
+        ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1));
+        sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1));
         w = matH.coeff(i,i) - p;
 
         if (m_eivalues.coeff(i).imag() < 0.0)
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f3c477374..ac4eb7633 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -120,7 +120,7 @@ ei_add_test(bandmatrix)
 ei_add_test(cholesky " " "${GSL_LIBRARIES}")
 ei_add_test(lu ${EI_OFLAG})
 ei_add_test(determinant)
-ei_add_test(inverse)
+ei_add_test(inverse ${EI_OFLAG})
 ei_add_test(qr)
 ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}")
 ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}")