From 50fed194f49c13e810add1655f65b1ae97c3ea68 Mon Sep 17 00:00:00 2001
From: Jim Turner <git@turner.link>
Date: Wed, 7 Apr 2021 22:43:27 -0400
Subject: [PATCH] Add determinant method to Cholesky

---
 src/linalg/cholesky.rs   | 10 ++++++++++
 tests/linalg/cholesky.rs | 18 ++++++++++++++++++
 2 files changed, 28 insertions(+)

diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs
index a6757b08..6bfbeb81 100644
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@@ -138,6 +138,16 @@ where
         self.solve_mut(&mut res);
         res
     }
+
+    /// Computes the determinant of the decomposed matrix.
+    pub fn determinant(&self) -> N {
+        let dim = self.chol.nrows();
+        let mut prod_diag = N::one();
+        for i in 0..dim {
+            prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) };
+        }
+        prod_diag * prod_diag
+    }
 }
 
 impl<N: ComplexField, D: Dim> Cholesky<N, D>
diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs
index 5ea0edaf..99893bd4 100644
--- a/tests/linalg/cholesky.rs
+++ b/tests/linalg/cholesky.rs
@@ -79,6 +79,24 @@ macro_rules! gen_tests(
                     prop_assert!(id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7));
                 }
 
+                #[test]
+                fn cholesky_determinant(n in PROPTEST_MATRIX_DIM) {
+                    let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
+                    let lu_det = m.clone().lu().determinant();
+                    let chol_det = m.cholesky().unwrap().determinant();
+
+                    prop_assert!(relative_eq!(lu_det, chol_det, epsilon = 1.0e-7));
+                }
+
+                #[test]
+                fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) {
+                    let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
+                    let lu_det = m.clone().lu().determinant();
+                    let chol_det = m.cholesky().unwrap().determinant();
+
+                    prop_assert!(relative_eq!(lu_det, chol_det, epsilon = 1.0e-7));
+                }
+
                 #[test]
                 fn cholesky_rank_one_update(_n in PROPTEST_MATRIX_DIM) {
                     let mut m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();