diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 4d667fbc1..bf7d8e8f8 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -235,10 +235,16 @@ where } let sqrt_denom = |v: T| { - if v.is_zero() { + // positive definiteness requires the diagonal to be real and strictly positive + if !v.clone().imaginary().is_zero() { return None; } - v.try_sqrt() + let re = v.real(); + if re <= T::RealField::zero() { + return None; + } + // fixme: could be improved by sqrt on real, e.g. `T::from_real(re.sqrt())` + Some(T::from_real(re).sqrt()) }; let diag = unsafe { matrix.get_unchecked((j, j)).clone() }; diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index bd135fb6a..d0ee92b9a 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -11,6 +11,20 @@ fn cholesky_with_substitute() { assert!(na::Cholesky::new_with_substitute(m, 1e-8).is_some()); } +// 2x2 matrix filled with -1 is not positive definite and should not return a Cholesky +// decomposition (regardless of element type). +// See https://github.com/dimforge/nalgebra/issues/1536 . +#[test] +fn cholesky_check_positive_definiteness() { + let x = -1.0; + let m = na::Matrix2::new(x, x, x, x); + assert!(nalgebra::Cholesky::new(m).is_none()); + + let x = na::Complex::new(-1.0, 0.0); + let m = na::Matrix2::new(x, x, x, x); + assert!(nalgebra::Cholesky::new(m).is_none()); +} + macro_rules! gen_tests( ($module: ident, $scalar: ty) => { mod $module {