GC-81: Redundant condition check in specialized Inverse method.

KeithNel · Dec 29, 2013 · 250b767 · 250b767
1 parent d9c7d0f
commit 250b767
Show file tree

Hide file tree

Showing 2 changed files with 54 additions and 6 deletions.
diff --git a/Sources/Accord.Math/Matrix/Matrix.Linear.cs b/Sources/Accord.Math/Matrix/Matrix.Linear.cs
@@ -215,16 +215,21 @@ public static double[] Solve(this double[,] matrix, double[] rightSide, bool lea
             if (rows != cols)
                 throw new ArgumentException("Matrix must be square", "matrix");
 
-            if (rows == 3 && rows == 3)
+            if (rows == 3)
             {
                 // Special case for 3x3 matrices
                 double a = matrix[0, 0], b = matrix[0, 1], c = matrix[0, 2];
                 double d = matrix[1, 0], e = matrix[1, 1], f = matrix[1, 2];
                 double g = matrix[2, 0], h = matrix[2, 1], i = matrix[2, 2];
 
-                double m = 1.0 / (a * (e * i - f * h) -
-                                  b * (d * i - f * g) +
-                                  c * (d * h - e * g));
+                double den = a * (e * i - f * h) -
+                             b * (d * i - f * g) +
+                             c * (d * h - e * g);
+
+                if (den == 0)
+                    throw new SingularMatrixException();
+
+                double m = 1.0 / den;
 
                 double[,] inv = (inPlace) ? matrix : new double[3, 3];
                 inv[0, 0] = m * (e * i - f * h);
@@ -236,22 +241,29 @@ public static double[] Solve(this double[,] matrix, double[] rightSide, bool lea
                 inv[2, 0] = m * (d * h - e * g);
                 inv[2, 1] = m * (b * g - a * h);
                 inv[2, 2] = m * (a * e - b * d);
+
                 return inv;
             }
 
-            if (rows == 2 && cols == 2)
+            if (rows == 2)
             {
                 // Special case for 2x2 matrices
                 double a = matrix[0, 0], b = matrix[0, 1];
                 double c = matrix[1, 0], d = matrix[1, 1];
 
-                double m = 1.0 / (a * d - b * c);
+                double den = a * d - b * c;
+
+                if (den == 0)
+                    throw new SingularMatrixException();
+
+                double m = 1.0 / den;
 
                 double[,] inv = (inPlace) ? matrix : new double[2, 2];
                 inv[0, 0] = +m * d;
                 inv[0, 1] = -m * b;
                 inv[1, 0] = -m * c;
                 inv[1, 1] = +m * a;
+
                 return inv;
             }
 

diff --git a/Sources/Accord.Tests/Accord.Tests.Math/Matrix/MatrixTest.cs b/Sources/Accord.Tests/Accord.Tests.Math/Matrix/MatrixTest.cs
@@ -29,6 +29,7 @@ namespace Accord.Tests.Math
     using System.Data;
     using AForge;
     using System.Linq;
+    using Accord.Math.Decompositions;
 
     [TestClass()]
     public partial class MatrixTest
@@ -1261,6 +1262,41 @@ public void CartesianProductTest3()
         #endregion
 
         #region Inverse, division and solving
+        [TestMethod()]
+        public void InverseTest2x2()
+        {
+            double[,] value = 
+            { 
+                { 3.0, 1.0 },
+                { 2.0, 2.0 }  
+            };
+
+            double[,] expected = new SingularValueDecomposition(value).Inverse();
+
+            double[,] actual = Matrix.Inverse(value);
+
+            Assert.IsTrue(Matrix.IsEqual(expected, actual, 1e-6));
+        }
+
+        [TestMethod()]
+        public void InverseTest3x3()
+        {
+            double[,] value = 
+            { 
+                { 6.0, 1.0, 2.0 },
+                { 0.0, 8.0, 1.0 },  
+                { 2.0, 4.0, 5.0 }  
+            };
+
+            Assert.IsFalse(value.IsSingular());
+
+            double[,] expected = new SingularValueDecomposition(value).Inverse();
+
+            double[,] actual = Matrix.Inverse(value);
+
+            Assert.IsTrue(Matrix.IsEqual(expected, actual, 1e-6));
+        }
+
         [TestMethod()]
         public void PseudoInverse()
         {