M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #861 from jturner314/fix-cholesky-det-complex 

Fix Cholesky::determinant for complex elements
This commit is contained in:
Sébastien Crozet 2021-04-11 13:49:35 +02:00 committed by GitHub
parent b9a524427b cea3bdc8e5
commit 2ebba30117
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 7 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/linalg/cholesky.rs
View File

@ -140,13 +140,13 @@ where
} }
/// Computes the determinant of the decomposed matrix. /// Computes the determinant of the decomposed matrix.
pub fn determinant(&self) -> N { pub fn determinant(&self) -> N::SimdRealField {
let dim = self.chol.nrows(); let dim = self.chol.nrows();
let mut prod_diag = N::one(); let mut prod_diag = N::one();
for i in 0..dim { for i in 0..dim {
prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) }; prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) };
} }
prod_diag * prod_diag prod_diag.simd_modulus_squared()
} }
} }

7
tests/linalg/cholesky.rs
View File

@ -7,6 +7,7 @@ macro_rules! gen_tests(
use na::dimension::{U4, Dynamic}; use na::dimension::{U4, Dynamic};
use na::{DMatrix, DVector, Matrix4x3, Vector4}; use na::{DMatrix, DVector, Matrix4x3, Vector4};
use rand::random; use rand::random;
use simba::scalar::ComplexField;
#[allow(unused_imports)] #[allow(unused_imports)]
use crate::core::helper::{RandScalar, RandComplex}; use crate::core::helper::{RandScalar, RandComplex};
@ -83,18 +84,20 @@ macro_rules! gen_tests(
fn cholesky_determinant(n in PROPTEST_MATRIX_DIM) { fn cholesky_determinant(n in PROPTEST_MATRIX_DIM) {
let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
let lu_det = m.clone().lu().determinant(); let lu_det = m.clone().lu().determinant();
assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7);
let chol_det = m.cholesky().unwrap().determinant(); let chol_det = m.cholesky().unwrap().determinant();
prop_assert!(relative_eq!(lu_det, chol_det, epsilon = 1.0e-7)); prop_assert!(relative_eq!(lu_det.real(), chol_det, epsilon = 1.0e-7));
} }
#[test] #[test]
fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) { fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) {
let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
let lu_det = m.clone().lu().determinant(); let lu_det = m.clone().lu().determinant();
assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7);
let chol_det = m.cholesky().unwrap().determinant(); let chol_det = m.cholesky().unwrap().determinant();
prop_assert!(relative_eq!(lu_det, chol_det, epsilon = 1.0e-7)); prop_assert!(relative_eq!(lu_det.real(), chol_det, epsilon = 1.0e-7));
} }
#[test] #[test]