srenblad/nalgebra
1
0
Fork 0
forked from M-Labs/nalgebra
Code Pull Requests Activity

Relax invertibility test in try_inverse()

Browse Source 
The previous implementation of try_inverse() used an approximate
check of the determinant against 0 for small matrices to
determine if the matrix was invertible. This is not a reliable test,
and may fail for perfectly invertible matrices. This change
simply makes the test criterion an exact comparison instead.
This commit is contained in:
Andreas Longva 2017-04-27 22:19:32 +02:00 committed by Sébastien Crozet
parent 9489e8f97e
commit a52b079578
2 changed files with 63 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
src/core/inverse.rs
View File

@ -45,18 +45,19 @@ impl<N, D: Dim, S> SquareMatrix<N, D, S>
match dim { match dim {
0 => true, 0 => true,
1 => { 1 => {
if relative_eq!(self.get_unchecked(0, 0), &N::zero()) { let determinant = self.get_unchecked(0, 0).clone();
if determinant == N::zero() {
false false
} }
else { else {
*self.get_unchecked_mut(0, 0) = N::one() / self.determinant(); *self.get_unchecked_mut(0, 0) = N::one() / determinant;
true true
} }
}, },
2 => { 2 => {
let determinant = self.determinant(); let determinant = self.determinant();
if relative_eq!(&determinant, &N::zero()) { if determinant == N::zero() {
false false
} }
else { else {
@ -94,7 +95,7 @@ impl<N, D: Dim, S> SquareMatrix<N, D, S>
m12 * minor_m11_m23 + m12 * minor_m11_m23 +
m13 * minor_m11_m22; m13 * minor_m11_m22;
if relative_eq!(&determinant, &N::zero()) { if determinant == N::zero() {
false false
} }
else { else {

58
tests/matrix_inverse.rs
View File

@ -61,3 +61,61 @@ fn matrix5_try_inverse() {
assert_relative_eq!(a_inv, expected_inverse, max_relative=1e-4); 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);
}