add polar decomposition on both sides, in SVD, with test

Eigen/src/SVD/SVD.h

diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index 4705855..9d26916 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -79,6 +79,9 @@
     void compute(const MatrixType& matrix);
     SVD& sort();
 
+    void computeUnitaryPositive(MatrixUType *unitary, MatrixType *positive) const;
+    void computePositiveUnitary(MatrixType *positive, MatrixVType *unitary) const;
+
   protected:
     /** \internal */
     MatrixUType m_matU;
@@ -534,6 +537,36 @@
   return true;
 }
 
+/** Computes the polar decomposition of the matrix, as a product unitary x positive.
+  *
+  * If either pointer is zero, the corresponding computation is skipped.
+  *
+  * Only for square matrices.
+  */
+template<typename MatrixType>
+void SVD<MatrixType>::computeUnitaryPositive(typename SVD<MatrixType>::MatrixUType *unitary,
+                                             MatrixType *positive) const
+{
+  ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices");
+  if(unitary) *unitary = m_matU * m_matV.adjoint();
+  if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint();
+}
+
+/** Computes the polar decomposition of the matrix, as a product positive x unitary.
+  *
+  * If either pointer is zero, the corresponding computation is skipped.
+  *
+  * Only for square matrices.
+  */
+template<typename MatrixType>
+void SVD<MatrixType>::computePositiveUnitary(MatrixType *positive,
+                                             typename SVD<MatrixType>::MatrixVType *unitary) const
+{
+  ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
+  if(unitary) *unitary = m_matU * m_matV.adjoint();
+  if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint();
+}
+
 /** \svd_module
   * \returns the SVD decomposition of \c *this
   */

test/svd.cpp

diff --git a/test/svd.cpp b/test/svd.cpp
index 4776ef7..688c3f4 100644
--- a/test/svd.cpp
+++ b/test/svd.cpp
@@ -44,13 +44,15 @@
   if (ei_is_same_type<RealScalar,float>::ret)
     largerEps = 1e-3f;
 
-  SVD<MatrixType> svd(a);
-  MatrixType sigma = MatrixType::Zero(rows,cols);
-  MatrixType matU  = MatrixType::Zero(rows,rows);
-  sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal();
-  matU.block(0,0,rows,cols) = svd.matrixU();
+  {
+    SVD<MatrixType> svd(a);
+    MatrixType sigma = MatrixType::Zero(rows,cols);
+    MatrixType matU  = MatrixType::Zero(rows,rows);
+    sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal();
+    matU.block(0,0,rows,cols) = svd.matrixU();
+    VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose());
+  }
 
-  VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose());
 
   if (rows==cols)
   {
@@ -63,6 +65,24 @@
     svd.solve(b, &x);
     VERIFY_IS_APPROX(a * x,b);
   }
+
+
+  if(rows==cols)
+  {
+    SVD<MatrixType> svd(a);
+    MatrixType unitary, positive;
+    svd.computeUnitaryPositive(&unitary, &positive);
+    VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows()));
+    VERIFY_IS_APPROX(positive, positive.adjoint());
+    for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity
+    VERIFY_IS_APPROX(unitary*positive, a);
+
+    svd.computePositiveUnitary(&positive, &unitary);
+    VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows()));
+    VERIFY_IS_APPROX(positive, positive.adjoint());
+    for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity
+    VERIFY_IS_APPROX(positive*unitary, a);
+  }
 }
 
 void test_svd()