From aeb9f7ea3994499540e3a7eeb1b4eed46768a275 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Crozet=20S=C3=A9bastien?= <developer@crozet.re>
Date: Thu, 25 Feb 2021 13:20:20 +0100
Subject: [PATCH] Add a matrix.udu() method to compute the UDU decomposition.

---
 src/linalg/decomposition.rs | 17 +++++++++++++++--
 tests/linalg/udu.rs         | 13 +++++++------
 2 files changed, 22 insertions(+), 8 deletions(-)

diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs
index 7890748c..99f2bcfe 100644
--- a/src/linalg/decomposition.rs
+++ b/src/linalg/decomposition.rs
@@ -1,8 +1,8 @@
 use crate::storage::Storage;
 use crate::{
     Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff,
-    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, Schur, SymmetricEigen,
-    SymmetricTridiagonal, LU, QR, SVD, U1,
+    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, RealField, Schur, SymmetricEigen,
+    SymmetricTridiagonal, LU, QR, SVD, U1, UDU
 };
 
 /// # Rectangular matrix decomposition
@@ -134,6 +134,7 @@ impl<N: ComplexField, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
 /// | -------------------------|---------------------------|--------------|
 /// | Hessenberg               | `Q * H * Qᵀ`             | `Q` is a unitary matrix and `H` an upper-Hessenberg matrix. |
 /// | Cholesky                 | `L * Lᵀ`                 | `L` is a lower-triangular matrix. |
+/// | UDU                      | `U * D * Uᵀ`             | `U` is a upper-triangular matrix, and `D` a diagonal matrix. |
 /// | Schur decomposition      | `Q * T * Qᵀ`             | `Q` is an unitary matrix and `T` a quasi-upper-triangular matrix. |
 /// | Symmetric eigendecomposition | `Q ~ Λ ~ Qᵀ`   | `Q` is an unitary matrix, and `Λ` is a real diagonal matrix. |
 /// | Symmetric tridiagonalization | `Q ~ T ~ Qᵀ`   | `Q` is an unitary matrix, and `T` is a tridiagonal matrix. |
@@ -149,6 +150,18 @@ impl<N: ComplexField, D: Dim, S: Storage<N, D, D>> Matrix<N, D, D, S> {
         Cholesky::new(self.into_owned())
     }
 
+    /// Attempts to compute the UDU decomposition of this matrix.
+    ///
+    /// The input matrix `self` is assumed to be symmetric and this decomposition will only read
+    /// the upper-triangular part of `self`.
+    pub fn udu(self) -> UDU<N, D>
+    where
+        N: RealField,
+        DefaultAllocator: Allocator<N, D> + Allocator<N, D, D>,
+    {
+        UDU::new(self.into_owned())
+    }
+
     /// Computes the Hessenberg decomposition of this matrix using householder reflections.
     pub fn hessenberg(self) -> Hessenberg<N, D>
     where
diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs
index 95cdc9c2..006cc95f 100644
--- a/tests/linalg/udu.rs
+++ b/tests/linalg/udu.rs
@@ -1,4 +1,4 @@
-use na::{Matrix3, UDU};
+use na::Matrix3;
 
 #[test]
 #[rustfmt::skip]
@@ -8,7 +8,8 @@ fn udu_simple() {
        -1.0,  2.0, -1.0,
         0.0, -1.0,  2.0);
 
-    let udu = UDU::new(m);
+    let udu = m.udu();
+    
     // Rebuild
     let p = udu.u * udu.d_matrix() * udu.u.transpose();
 
@@ -24,7 +25,7 @@ fn udu_non_sym_panic() {
         1.0, -2.0,  3.0,
        -2.0,  1.0,  0.0);
 
-    let udu = UDU::new(m);
+    let udu = m.udu();
     // Rebuild
     let p = udu.u * udu.d_matrix() * udu.u.transpose();
 
@@ -39,7 +40,7 @@ mod quickcheck_tests {
     macro_rules! gen_tests(
         ($module: ident, $scalar: ty) => {
             mod $module {
-                use na::{UDU, DMatrix, Matrix4};
+                use na::{DMatrix, Matrix4};
                 #[allow(unused_imports)]
                 use crate::core::helper::{RandScalar, RandComplex};
 
@@ -48,7 +49,7 @@ mod quickcheck_tests {
                         let n = std::cmp::max(1, std::cmp::min(n, 10));
                         let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part();
 
-                        let udu = UDU::new(m.clone());
+                        let udu = m.clone().udu();
                         let p = &udu.u * &udu.d_matrix() * &udu.u.transpose();
 
                         relative_eq!(m, p, epsilon = 1.0e-7)
@@ -57,7 +58,7 @@ mod quickcheck_tests {
                     fn udu_static(m: Matrix4<$scalar>) -> bool {
                         let m = m.map(|e| e.0).hermitian_part();
 
-                        let udu = UDU::new(m.clone());
+                        let udu = m.udu();
                         let p = udu.u * udu.d_matrix() * udu.u.transpose();
 
                         relative_eq!(m, p, epsilon = 1.0e-7)