* bugfixes in Product, and test/product_selfadjoint

* speed up in the extraction of the matrix Q in Tridiagonalization

5 files changed, 33 insertions, 18 deletions

diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h
index 8136724dd..7a37b3624 100644
--- a/Eigen/src/Core/BandMatrix.h
+++ b/Eigen/src/Core/BandMatrix.h
@@ -198,7 +198,7 @@ class BandMatrix : public MultiplierBase<BandMatrix<_Scalar,Rows,Cols,Supers,Sub
   * \sa class BandMatrix
   */
 template<typename Scalar, int Size, int Options>
-class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,1,Options&SelfAdjoint?0:1,Options|RowMajor>
+class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor>
 {
     typedef BandMatrix<Scalar,Size,Size,1,Options&SelfAdjoint?0:1,Options|RowMajor> Base;
   public:
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index 6feede458..93e318035 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -880,7 +880,7 @@ inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien
   }
   else
   {
-    lazyAssign(static_cast<const MatrixBase<Product<Lhs,Rhs,CacheFriendlyProduct> > &>(product));
+    lazyAssign(Product<Lhs,Rhs,NormalProduct>(product.lhs(),product.rhs()));
   }
   return derived();
 }
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index ff6b98aa1..b003fb4d7 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -50,6 +50,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
     typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
+//     typedef typename TridiagonalizationType::TridiagonalMatrixType TridiagonalMatrixType;
 
     SelfAdjointEigenSolver()
         : m_eivec(int(Size), int(Size)),
diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h
index 8051f7a55..cf737bd30 100644
--- a/Eigen/src/QR/Tridiagonalization.h
+++ b/Eigen/src/QR/Tridiagonalization.h
@@ -121,8 +121,9 @@ template<typename _MatrixType> class Tridiagonalization
       */
     inline const MatrixType& packedMatrix(void) const { return m_matrix; }
 
-    MatrixType matrixQ(void) const;
-    MatrixType matrixT(void) const;
+    MatrixType matrixQ() const;
+    template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const;
+    MatrixType matrixT() const;
     const DiagonalReturnType diagonal(void) const;
     const SubDiagonalReturnType subDiagonal(void) const;
 
@@ -270,21 +271,32 @@ template<typename MatrixType>
 typename Tridiagonalization<MatrixType>::MatrixType
 Tridiagonalization<MatrixType>::matrixQ(void) const
 {
+  MatrixType matQ;
+  matrixQInPlace(&matQ);
+  return matQ;
+}
+
+template<typename MatrixType>
+template<typename QDerived>
+void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) const
+{
+  QDerived& matQ = q->derived();
   int n = m_matrix.rows();
-  MatrixType matQ = MatrixType::Identity(n,n);
+  matQ = MatrixType::Identity(n,n);
+  Matrix<Scalar,1,Dynamic> aux(n);
   for (int i = n-2; i>=0; i--)
   {
     Scalar tmp = m_matrix.coeff(i+1,i);
     m_matrix.const_cast_derived().coeffRef(i+1,i) = 1;
 
-    // TODO this product could be optimized by processing the submatrix per panel of at least 4 columns
-    matQ.corner(BottomRight,n-i-1,n-i-1) -=
-      ((m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1)) *
-      (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy();
+    aux.end(n-i-1) = (m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy();
+    // rank one update, TODO ! make it works efficiently as expected
+    for (int j=i+1;j<n;++j)
+      matQ.col(j).end(n-i-1) -= ( aux.coeff(j)) * m_matrix.col(i).end(n-i-1);
+//     matQ.corner(BottomRight,n-i-1,n-i-1) -= (m_matrix.col(i).end(n-i-1) * aux.end(n-i-1)).lazy();
 
     m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp;
   }
-  return matQ;
 }
 
 /** Performs a full decomposition in place */
@@ -303,7 +315,7 @@ void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalT
     diag = tridiag.diagonal();
     subdiag = tridiag.subDiagonal();
     if (extractQ)
-      mat = tridiag.matrixQ();
+      tridiag.matrixQInPlace(&mat);
   }
 }
 
diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp
index 297bab1a9..29fbf11bf 100644
--- a/test/product_selfadjoint.cpp
+++ b/test/product_selfadjoint.cpp
@@ -77,12 +77,14 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
   m2.template selfadjointView<UpperTriangular>().rank2update(-r1.adjoint(),r2.adjoint()*s3,s1);
   VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<UpperTriangular>().toDense());
 
-  m2 = m1.template triangularView<LowerTriangular>();
-  m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rank2update(v1.end(rows-1),v2.start(cols-1));
-  m3 = m1;
-  m3.block(1,1,rows-1,cols-1) += v1.end(rows-1) * v2.start(cols-1).adjoint()+ v2.start(cols-1) * v1.end(rows-1).adjoint();
-  VERIFY_IS_APPROX(m2, m3.template triangularView<LowerTriangular>().toDense());
-
+  if (rows>1)
+  {
+    m2 = m1.template triangularView<LowerTriangular>();
+    m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rank2update(v1.end(rows-1),v2.start(cols-1));
+    m3 = m1;
+    m3.block(1,1,rows-1,cols-1) += v1.end(rows-1) * v2.start(cols-1).adjoint()+ v2.start(cols-1) * v1.end(rows-1).adjoint();
+    VERIFY_IS_APPROX(m2, m3.template triangularView<LowerTriangular>().toDense());
+  }
 }
 
 void test_product_selfadjoint()
@@ -93,7 +95,7 @@ void test_product_selfadjoint()
     CALL_SUBTEST( product_selfadjoint(Matrix3d()) );
     CALL_SUBTEST( product_selfadjoint(MatrixXcf(4, 4)) );
     CALL_SUBTEST( product_selfadjoint(MatrixXcd(21,21)) );
-    CALL_SUBTEST( product_selfadjoint(MatrixXd(4,4)) );
+    CALL_SUBTEST( product_selfadjoint(MatrixXd(14,14)) );
     CALL_SUBTEST( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(17,17)) );
     CALL_SUBTEST( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) );
   }