change the make householder algorithm so that the remaining coefficient

is real, and make Tridiagonalization use it

1 files changed, 21 insertions, 10 deletions

diff --git a/test/householder.cpp b/test/householder.cpp
index d8b11c9f1..7d300899f 100644
--- a/test/householder.cpp
+++ b/test/householder.cpp
@@ -27,6 +27,8 @@
 
 template<typename MatrixType> void householder(const MatrixType& m)
 {
+  static bool even = true;
+  even = !even;
   /* this test covers the following files:
      Householder.h
   */
@@ -38,46 +40,55 @@ template<typename MatrixType> void householder(const MatrixType& m)
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
   typedef Matrix<Scalar, ei_decrement_size<MatrixType::RowsAtCompileTime>::ret, 1> EssentialVectorType;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
+
+  Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols));
+  Scalar* tmp = &_tmp.coeffRef(0,0);
   
-  RealScalar beta;
+  Scalar beta;
+  RealScalar alpha;
   EssentialVectorType essential;
 
   VectorType v1 = VectorType::Random(rows), v2;
   v2 = v1;
-  v1.makeHouseholder(&essential, &beta);
-  v1.applyHouseholderOnTheLeft(essential,beta);
-  
+  v1.makeHouseholder(&essential, &beta, &alpha);
+  v1.applyHouseholderOnTheLeft(essential,beta,tmp);
   VERIFY_IS_APPROX(v1.norm(), v2.norm());
   VERIFY_IS_MUCH_SMALLER_THAN(v1.end(rows-1).norm(), v1.norm());
   v1 = VectorType::Random(rows);
   v2 = v1;
-  v1.applyHouseholderOnTheLeft(essential,beta);
+  v1.applyHouseholderOnTheLeft(essential,beta,tmp);
   VERIFY_IS_APPROX(v1.norm(), v2.norm());
   
   MatrixType m1(rows, cols),
              m2(rows, cols);
 
   v1 = VectorType::Random(rows);
+  if(even) v1.end(rows-1).setZero();
   m1.colwise() = v1;
   m2 = m1;
-  m1.col(0).makeHouseholder(&essential, &beta);
-  m1.applyHouseholderOnTheLeft(essential,beta);
+  m1.col(0).makeHouseholder(&essential, &beta, &alpha);
+  m1.applyHouseholderOnTheLeft(essential,beta,tmp);
   VERIFY_IS_APPROX(m1.norm(), m2.norm());
   VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm());
+  VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m1(0,0)), ei_real(m1(0,0)));
+  VERIFY_IS_APPROX(ei_real(m1(0,0)), alpha);
   
   v1 = VectorType::Random(rows);
+  if(even) v1.end(rows-1).setZero();
   SquareMatrixType m3(rows,rows), m4(rows,rows);
   m3.rowwise() = v1.transpose();
   m4 = m3;
-  m3.row(0).makeHouseholder(&essential, &beta);
-  m3.applyHouseholderOnTheRight(essential,beta);
+  m3.row(0).makeHouseholder(&essential, &beta, &alpha);
+  m3.applyHouseholderOnTheRight(essential,beta,tmp);
   VERIFY_IS_APPROX(m3.norm(), m4.norm());
   VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
+  VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0)));
+  VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha);
 }
 
 void test_householder()
 {
-  for(int i = 0; i < g_repeat; i++) {
+  for(int i = 0; i < 2*g_repeat; i++) {
     CALL_SUBTEST( householder(Matrix<double,2,2>()) );
     CALL_SUBTEST( householder(Matrix<float,2,3>()) );
     CALL_SUBTEST( householder(Matrix<double,3,5>()) );