diff --git a/src/core/inverse.rs b/src/core/inverse.rs index ed066e49..488c5b5b 100644 --- a/src/core/inverse.rs +++ b/src/core/inverse.rs @@ -45,18 +45,19 @@ impl SquareMatrix match dim { 0 => true, 1 => { - if relative_eq!(self.get_unchecked(0, 0), &N::zero()) { + let determinant = self.get_unchecked(0, 0).clone(); + if determinant == N::zero() { false } else { - *self.get_unchecked_mut(0, 0) = N::one() / self.determinant(); + *self.get_unchecked_mut(0, 0) = N::one() / determinant; true } }, 2 => { let determinant = self.determinant(); - if relative_eq!(&determinant, &N::zero()) { + if determinant == N::zero() { false } else { @@ -94,7 +95,7 @@ impl SquareMatrix m12 * minor_m11_m23 + m13 * minor_m11_m22; - if relative_eq!(&determinant, &N::zero()) { + if determinant == N::zero() { false } else { diff --git a/tests/matrix_inverse.rs b/tests/matrix_inverse.rs index a37b54ff..fe1af112 100644 --- a/tests/matrix_inverse.rs +++ b/tests/matrix_inverse.rs @@ -61,3 +61,61 @@ fn matrix5_try_inverse() { assert_relative_eq!(a_inv, expected_inverse, max_relative=1e-4); } + +#[test] +fn matrix1_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix1::new(1.0e-20); + let expected_inverse = Matrix1::new(1.0e20); + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +} + +#[test] +fn matrix2_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix2::new(1.0e-20, 0.0, + 0.0, 1.0e-20); + let expected_inverse = Matrix2::new(1.0e20, 0.0, + 0.0, 1.0e20); + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +} + +#[test] +fn matrix3_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix3::new(1.0e-20, 0.0, 0.0, + 0.0, 1.0e-20, 0.0, + 0.0, 0.0, 1.0e-20); + let expected_inverse = Matrix3::new(1.0e20, 0.0, 0.0, + 0.0, 1.0e20, 0.0, + 0.0, 0.0, 1.0e20); + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +} + +#[test] +fn matrix5_try_inverse_scaled_identity() { + // A perfectly invertible matrix with + // very small coefficients + let a = Matrix5::new(1.0e-20, 0.0, 0.0, 0.0, 0.0, + 0.0, 1.0e-20, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0e-20, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.0e-20, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0e-20); + let expected_inverse = Matrix5::new(1.0e+20, 0.0, 0.0, 0.0, 0.0, + 0.0, 1.0e+20, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0e+20, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.0e+20, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0e+20);; + let a_inv = a.try_inverse().expect("Matrix is invertible"); + + assert_relative_eq!(a_inv, expected_inverse); +}