Fix compilation of HouseholderQR and ColPivotingHouseholderQR for non-square fixed-size matrices.

For Colpiv that was just changing MatrixQType to MatrixType in the instantiation of HouseholderSequence.
For HouseholderQR I also re-ported the solve method from Colpiv as there were multiple issues.

4 files changed, 36 insertions, 27 deletions

diff --git a/Eigen/src/QR/ColPivotingHouseholderQR.h b/Eigen/src/QR/ColPivotingHouseholderQR.h
index 0898b5d1f..b141da0aa 100644
--- a/Eigen/src/QR/ColPivotingHouseholderQR.h
+++ b/Eigen/src/QR/ColPivotingHouseholderQR.h
@@ -62,7 +62,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
     typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
-    typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
+    typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
 
     /**
     * \brief Default Constructor.
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index 39502a30f..4cc553926 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -62,7 +62,7 @@ template<typename MatrixType> class HouseholderQR
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType;
     typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
-    typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
+    typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
 
     /**
     * \brief Default Constructor.
@@ -206,18 +206,22 @@ void HouseholderQR<MatrixType>::solve(
 ) const
 {
   ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
+  result->resize(m_qr.cols(), b.cols());
   const int rows = m_qr.rows();
-  const int cols = b.cols();
+  const int rank = std::min(m_qr.rows(), m_qr.cols());
   ei_assert(b.rows() == rows);
-  result->resize(rows, cols);
 
-  *result = b;
-  result->applyOnTheLeft(matrixQAsHouseholderSequence().inverse());
+  typename OtherDerived::PlainMatrixType c(b);
+
+  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+  c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,rank), m_hCoeffs.start(rank)).transpose());
 
-  const int rank = std::min(result->rows(), result->cols());
   m_qr.corner(TopLeft, rank, rank)
       .template triangularView<UpperTriangular>()
-      .solveInPlace(result->corner(TopLeft, rank, result->cols()));
+      .solveInPlace(c.corner(TopLeft, rank, c.cols()));
+
+  result->corner(TopLeft, rank, c.cols()) = c.corner(TopLeft,rank, c.cols());
+  result->corner(BottomLeft, result->rows()-rank, c.cols()).setZero();
 }
 
 /** \returns the matrix Q */
diff --git a/test/qr.cpp b/test/qr.cpp
index f185ac86e..036a3c9f2 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -41,20 +41,26 @@ template<typename MatrixType> void qr(const MatrixType& m)
   for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
 
   VERIFY_IS_APPROX(a, qrOfA.matrixQ() * r);
+}
 
-  SquareMatrixType b = a.adjoint() * a;
+template<typename MatrixType, int Cols2> void qr_fixedsize()
+{
+  enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
+  typedef typename MatrixType::Scalar Scalar;
+  Matrix<Scalar,Rows,Cols> m1 = Matrix<Scalar,Rows,Cols>::Random();
+  HouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
 
-  // check tridiagonalization
-  Tridiagonalization<SquareMatrixType> tridiag(b);
-  VERIFY_IS_APPROX(b, tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
+  Matrix<Scalar,Rows,Cols> r = qr.matrixQR();
+  // FIXME need better way to construct trapezoid
+  for(int i = 0; i < Rows; i++) for(int j = 0; j < Cols; j++) if(i>j) r(i,j) = Scalar(0);
 
-  // check hessenberg decomposition
-  HessenbergDecomposition<SquareMatrixType> hess(b);
-  VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
-  VERIFY_IS_APPROX(tridiag.matrixT(), hess.matrixH());
-  b = SquareMatrixType::Random(cols,cols);
-  hess.compute(b);
-  VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
+  VERIFY_IS_APPROX(m1, qr.matrixQ() * r);
+
+  Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
+  Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
+  m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
+  qr.solve(m3, &m2);
+  VERIFY_IS_APPROX(m3, m1*m2);
 }
 
 template<typename MatrixType> void qr_invertible()
@@ -105,11 +111,11 @@ template<typename MatrixType> void qr_verify_assert()
 void test_qr()
 {
   for(int i = 0; i < 1; i++) {
-    // FIXME : very weird bug here
-//     CALL_SUBTEST( qr(Matrix2f()) );
-    CALL_SUBTEST( qr(Matrix4d()) );
-    CALL_SUBTEST( qr(MatrixXf(47,40)) );
-    CALL_SUBTEST( qr(MatrixXcd(17,7)) );
+   CALL_SUBTEST( qr(MatrixXf(47,40)) );
+   CALL_SUBTEST( qr(MatrixXcd(17,7)) );
+   CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,4>, 2 >() ));
+   CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 4 >() ));
+   CALL_SUBTEST(( qr_fixedsize<Matrix<double,2,5>, 7 >() ));
   }
 
   for(int i = 0; i < g_repeat; i++) {
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index e0edad842..4b6f7dd6b 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -154,10 +154,9 @@ void test_qr_colpivoting()
     CALL_SUBTEST( qr<MatrixXf>() );
     CALL_SUBTEST( qr<MatrixXd>() );
     CALL_SUBTEST( qr<MatrixXcd>() );
+    CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
+    CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
   }
-  
-  CALL_SUBTEST(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
-  CALL_SUBTEST(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
 
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST( qr_invertible<MatrixXf>() );