Merge pull request #860 from jturner314/cholesky-determinant

Add determinant method to Cholesky

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>

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();