diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs
index 5fd3debc..dabcb491 100644
--- a/src/linalg/svd.rs
+++ b/src/linalg/svd.rs
@@ -643,6 +643,35 @@ where
     }
 }
 
+impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
+where
+    DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
+        + Allocator<T, R, DimMinimum<R, C>>
+        + Allocator<T::RealField, DimMinimum<R, C>>,
+{
+    /// converts SVD results to a polar form
+
+    pub fn to_polar(&self) -> Result<(OMatrix<T, R, R>, OMatrix<T, R, C>), &'static str>
+    where DefaultAllocator: Allocator<T, R, C> //result
+            + Allocator<T, DimMinimum<R, C>, R> // adjoint
+            + Allocator<T, DimMinimum<R, C>> // mapped vals
+            + Allocator<T, R, R> // square matrix & result
+            + Allocator<T, DimMinimum<R, C>, DimMinimum<R, C>> // ?
+        ,
+    {
+        match (&self.u, &self.v_t) {
+            (Some(u), Some(v_t)) => Ok((
+                u * OMatrix::from_diagonal(&self.singular_values.map(|e| T::from_real(e)))
+                    * u.adjoint(),
+                u * v_t,
+            )),
+            (None, None) => Err("SVD solve: U and V^t have not been computed."),
+            (None, _) => Err("SVD solve: U has not been computed."),
+            (_, None) => Err("SVD solve: V^t has not been computed."),
+        }
+    }
+}
+
 impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
 where
     DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs
index d76d1759..65a92ddb 100644
--- a/tests/linalg/svd.rs
+++ b/tests/linalg/svd.rs
@@ -441,3 +441,14 @@ fn svd_sorted() {
         epsilon = 1.0e-5
     );
 }
+
+#[test]
+fn svd_polar_decomposition() {
+
+    let m  = DMatrix::<f64>::new_random(4, 4);
+    let svd = m.clone().svd(true, true);
+    let (p,u) = svd.to_polar().unwrap();
+
+    assert_relative_eq!(m, p*u, epsilon = 1.0e-5);
+
+}