Add ol.math.solveLinearSystem

src/ol/math.js

@@ -78,6 +78,69 @@ ol.math.squaredDistance = function(x1, y1, x2, y2) {
 };
 
 
+/**
+ * Solves system of linear equations using Gaussian elimination method.
+ *
+ * @param {Array.<Array.<number>>} A Augmented matrix (n x n + 1 column)
+ *                                   in row-major order.
+ * @return {Array.<number>} The resulting vector.
+ */
+ol.math.solveLinearSystem = function(A) {
+  var n = A.length;
+
+  if (goog.asserts.ENABLE_ASSERTS) {
+    for (var row = 0; row < n; row++) {
+      goog.asserts.assert(A[row].length == n + 1,
+                          'every row should have correct number of columns');
+    }
+  }
+
+  for (var i = 0; i < n; i++) {
+    // Find max in the i-th column (ignoring i - 1 first rows)
+    var maxRow = i;
+    var maxEl = Math.abs(A[i][i]);
+    for (var r = i + 1; r < n; r++) {
+      var absValue = Math.abs(A[r][i]);
+      if (absValue > maxEl) {
+        maxEl = absValue;
+        maxRow = r;
+      }
+    }
+
+    if (maxEl === 0) {
+      return null; // matrix is singular
+    }
+
+    // Swap max row with i-th (current) row
+    var tmp = A[maxRow];
+    A[maxRow] = A[i];
+    A[i] = tmp;
+
+    // Subtract the i-th row to make all the remaining rows 0 in the i-th column
+    for (var j = i + 1; j < n; j++) {
+      var coef = -A[j][i] / A[i][i];
+      for (var k = i; k < n + 1; k++) {
+        if (i == k) {
+          A[j][k] = 0;
+        } else {
+          A[j][k] += coef * A[i][k];
+        }
+      }
+    }
+  }
+
+  // Solve Ax=b for upper triangular matrix A
+  var x = new Array(n);
+  for (var l = n - 1; l >= 0; l--) {
+    x[l] = A[l][n] / A[l][l];
+    for (var m = l - 1; m >= 0; m--) {
+      A[m][n] -= A[m][l] * x[l];
+    }
+  }
+  return x;
+};
+
+
 /**
  * Converts radians to to degrees.
  *

test/spec/ol/math.test.js

@@ -91,6 +91,45 @@ describe('ol.math.roundUpToPowerOfTwo', function() {
 });
 
 
+describe('ol.math.solveLinearSystem', function() {
+  it('calculates correctly', function() {
+    var result = ol.math.solveLinearSystem([
+      [2, 1, 3, 1],
+      [2, 6, 8, 3],
+      [6, 8, 18, 5]
+    ]);
+    expect(result[0]).to.roughlyEqual(0.3, 1e-9);
+    expect(result[1]).to.roughlyEqual(0.4, 1e-9);
+    expect(result[2]).to.roughlyEqual(0, 1e-9);
+  });
+  it('can handle singular matrix', function() {
+    var result = ol.math.solveLinearSystem([
+      [2, 1, 3, 1],
+      [2, 6, 8, 3],
+      [2, 1, 3, 1]
+    ]);
+    expect(result).to.be(null);
+  });
+  it('raises an exception when the matrix is malformed', function() {
+    expect(function() {
+      ol.math.solveLinearSystem([
+        [2, 1, 3, 1],
+        [2, 6, 8, 3],
+        [6, 8, 18]
+      ]);
+    }).to.throwException();
+
+    expect(function() {
+      ol.math.solveLinearSystem([
+        [2, 1, 3, 1],
+        [2, 6, 8, 3],
+        [6, 8, 18, 5, 0]
+      ]);
+    }).to.throwException();
+  });
+});
+
+
 describe('ol.math.toDegrees', function() {
   it('returns the correct value at -π', function() {
     expect(ol.math.toDegrees(-Math.PI)).to.be(-180);