From e6b77bcc6bc915ec38640ecf414726fa2ba56fba Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Wed, 2 Sep 2009 06:36:55 -0400
Subject: [PATCH] JacobiSVD: implement general R-SVD using full-pivoting QR, so
 we now support any rectangular matrix size by reducing to the smaller of the
 two dimensions (which is also an optimization)

---
 Eigen/SVD                                |  2 +-
 Eigen/src/QR/FullPivotingHouseholderQR.h |  2 +-
 Eigen/src/SVD/JacobiSVD.h                | 49 ++++++++++++++++++------
 test/jacobisvd.cpp                       |  8 ++--
 4 files changed, 45 insertions(+), 16 deletions(-)

diff --git a/Eigen/SVD b/Eigen/SVD
index 01582310c..8d66d0736 100644
--- a/Eigen/SVD
+++ b/Eigen/SVD
@@ -1,7 +1,7 @@
 #ifndef EIGEN_SVD_MODULE_H
 #define EIGEN_SVD_MODULE_H
 
-#include "Core"
+#include "QR"
 #include "Householder"
 #include "Jacobi"
 
diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h
index cee41b152..52405170b 100644
--- a/Eigen/src/QR/FullPivotingHouseholderQR.h
+++ b/Eigen/src/QR/FullPivotingHouseholderQR.h
@@ -286,7 +286,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co
   m_cols_permutation.resize(matrix.cols());
   int number_of_transpositions = 0;
   
-  RealScalar biggest;
+  RealScalar biggest(0);
   
   for (int k = 0; k < size; ++k)
   {
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index b3d289706..654f8981a 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -188,15 +188,42 @@ void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q,
 template<typename MatrixType, unsigned int Options>
 JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix)
 {
-  MatrixType work_matrix(matrix);
-  int size = matrix.rows();
-  if(ComputeU) m_matrixU = MatrixUType::Identity(size,size);
-  if(ComputeV) m_matrixV = MatrixUType::Identity(size,size);
-  m_singularValues.resize(size);
+  MatrixType work_matrix;
+  int rows = matrix.rows();
+  int cols = matrix.cols();
+  int diagSize = std::min(rows, cols);
+  if(ComputeU) m_matrixU = MatrixUType::Zero(rows,rows);
+  if(ComputeV) m_matrixV = MatrixVType::Zero(cols,cols);
+  m_singularValues.resize(diagSize);
   const RealScalar precision = 2 * epsilon<Scalar>();
 
+  if(rows > cols)
+  {
+    FullPivotingHouseholderQR<MatrixType> qr(matrix);
+    work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>();
+    if(ComputeU) m_matrixU = qr.matrixQ();
+    if(ComputeV)
+      for(int i = 0; i < cols; i++)
+        m_matrixV.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1);
+  }
+  else if(rows < cols)
+  {
+    FullPivotingHouseholderQR<MatrixType> qr(MatrixType(matrix.adjoint()));
+    work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>().adjoint();
+    if(ComputeV) m_matrixV = qr.matrixQ();
+    if(ComputeU)
+      for(int i = 0; i < rows; i++)
+        m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1);
+    
+  }
+  else
+  {
+    work_matrix = matrix;
+    if(ComputeU) m_matrixU.diagonal().setOnes();
+    if(ComputeV) m_matrixV.diagonal().setOnes();
+  }
 sweep_again:
-  for(int p = 1; p < size; ++p)
+  for(int p = 1; p < diagSize; ++p)
   {
     for(int q = 0; q < p; ++q)
     {
@@ -219,27 +246,27 @@ sweep_again:
 
   RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff();
   RealScalar maxAllowedOffDiag = biggestOnDiag * precision;
-  for(int p = 0; p < size; ++p)
+  for(int p = 0; p < diagSize; ++p)
   {
     for(int q = 0; q < p; ++q)
       if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
         goto sweep_again;
-    for(int q = p+1; q < size; ++q)
+    for(int q = p+1; q < diagSize; ++q)
       if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
         goto sweep_again;
   }
 
-  for(int i = 0; i < size; ++i)
+  for(int i = 0; i < diagSize; ++i)
   {
     RealScalar a = ei_abs(work_matrix.coeff(i,i));
     m_singularValues.coeffRef(i) = a;
     if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a;
   }
 
-  for(int i = 0; i < size; i++)
+  for(int i = 0; i < diagSize; i++)
   {
     int pos;
-    m_singularValues.end(size-i).maxCoeff(&pos);
+    m_singularValues.end(diagSize-i).maxCoeff(&pos);
     if(pos)
     {
       pos += i;
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index 9a4d79e45..8b4c7584e 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -91,12 +91,14 @@ void test_jacobisvd()
     CALL_SUBTEST( svd(Matrix3f()) );
     CALL_SUBTEST( svd(Matrix4d()) );
     CALL_SUBTEST( svd(MatrixXf(50,50)) );
-//    CALL_SUBTEST( svd(MatrixXd(14,7)) );
+    CALL_SUBTEST( svd(MatrixXcd(14,7)) );
+    CALL_SUBTEST( svd(MatrixXd(10,50)) );
+    
     CALL_SUBTEST( svd(MatrixXcf(3,3)) );
     CALL_SUBTEST( svd(MatrixXd(30,30)) );
   }
-  CALL_SUBTEST( svd(MatrixXf(200,200)) );
-  CALL_SUBTEST( svd(MatrixXcd(100,100)) );
+  CALL_SUBTEST( svd(MatrixXf(300,200)) );
+  CALL_SUBTEST( svd(MatrixXcd(100,150)) );
   
   CALL_SUBTEST( svd_verify_assert<Matrix3f>() );
   CALL_SUBTEST( svd_verify_assert<Matrix3d>() );