lusolve() - ccs

lib/function/algebra/solver/lusolve.js

@@ -73,6 +73,9 @@ function factory (type, config, load, typed) {
     // verify L, U, P
     l = _toMatrix(l);
     u = _toMatrix(u);
+    p = _toMatrix(p);
+    // modify column vector applying permutations
+    b = p.multiply(b);
     // use forward substitution to resolve L * y = b
     var y = l.forwardSubstitution(b);
     // use backward substitution to resolve U * x = y

lib/type/matrix/CcsMatrix.js

@@ -18,6 +18,7 @@ function factory (type, config, load, typed) {
   
   var abs = load(require('../../function/arithmetic/abs'));
   var add = load(require('../../function/arithmetic/add'));
+  var subtract = load(require('../../function/arithmetic/subtract'));
   var divideScalar = load(require('../../function/arithmetic/divideScalar'));
   var multiply = load(require('../../function/arithmetic/multiply'));
   var unaryMinus = load(require('../../function/arithmetic/unaryMinus'));
@@ -1541,6 +1542,148 @@ function factory (type, config, load, typed) {
     }
   };
   
+  /**
+   * Solves the linear equation system by forward substitution. Matrix must be a lower triangular matrix.
+   *
+   * M * x = b
+   *
+   * @param {Matrix, Array}         A column vector with the b values
+   *
+   * @return {Matrix}               A column vector with the linear system solution
+   */
+  CcsMatrix.prototype.forwardSubstitution = function (b) {
+    // validate matrix and vector, get vector accessor
+    var vget = this._substitutionValidation(b);
+    // rows & columns
+    var rows = this._size[0];
+    var columns = this._size[1];
+    // matrix arrays
+    var values = this._values;
+    var index = this._index;
+    var ptr = this._ptr;
+    // result (dense vector to speed up access)
+    var x = new Array(rows);
+    var xvalues = [];
+    var xindex = [];
+    var xptr = [];
+    // vars
+    var i, j, k, l;
+    // initialize x = b
+    for (i = 0; i < rows; i ++)
+      x[i] = vget(i);
+    // forward solve m * x = b
+    for (j = 0; j < columns; j++) {
+      // values in column
+      for (k = ptr[j], l = ptr[j + 1]; k < l; k++) {
+        // row
+        i = index[k];
+        // check row
+        if (i > j) {
+          // update x
+          x[i] = subtract(x[i], multiply(values[k], x[j]));
+        }
+      }
+    }
+    // init ptr
+    xptr.push(0);
+    // build column vector
+    for (i = 0; i < rows; i++) {
+      // x @ i
+      var xi = x[i];
+      // check for non zero
+      if (!equal(xi, 0)) {
+        // value @ row i
+        xvalues.push(xi);
+        // row
+        xindex.push(i);
+      }
+    }
+    // update ptr
+    xptr.push(xvalues.length);
+    // return column vector
+    return new CcsMatrix({
+      values: xvalues,
+      index: xindex,
+      ptr: xptr,
+      size: [rows, 1]
+    });
+  };
+  
+  /**
+   * Solves the linear equation system by backward substitution. Matrix must be an upper triangular matrix.
+   *
+   * M * x = b
+   *
+   * @param {Matrix, Array}         A column vector with the b values
+   *
+   * @return {Matrix}               A column vector with the linear system solution
+   */
+  CcsMatrix.prototype.backwardSubstitution = function (b) {
+    // validate matrix and vector, get vector accessor
+    var vget = this._substitutionValidation(b);
+    // rows & columns
+    var rows = this._size[0];
+    var columns = this._size[1];
+    // matrix arrays
+    var values = this._values;
+    var index = this._index;
+    var ptr = this._ptr;
+    // result (dense vector to speed up access)
+    var x = new Array(rows);
+    var xvalues = [];
+    var xindex = [];
+    var xptr = [];
+    // backward solve m * x = b
+    for (var i = rows - 1; i >= 0; i--) {
+      // update ptr
+      xptr.push(xvalues.length);
+      // initialize x
+      var xi = vget(i);
+      // value @ [i, i]
+      var vii = 0;
+      // values in row (backwards, upper triangular)
+      for (var k = ptr[i + 1] - 1, f = ptr[i]; k >= f; k--) {
+        // column
+        var j = index[k];
+        // check column
+        if (j > i) {
+          // update x
+          xi = subtract(xi, multiply(values[k], x[j]));
+        }
+        else if (j === i) {
+          // set vii value
+          vii = values[k];
+          // exit loop
+          break; 
+        }
+        else {
+          // system cannot be solved, there is no value @ [i, i]
+          throw new Error('Linear system cannot be solved since matrix is singular');
+        }
+      }
+      // update xi
+      xi = divideScalar(xi, vii);
+      // check for zero
+      if (!equal(xi, 0)) {
+        // push to values (at the beginning since we are looping backwards)
+        xvalues.unshift(xi);
+        // column vector
+        xindex.push(0);
+      }
+      // update dense vector
+      x[i] = xi;
+    }
+    // update ptr
+    xptr.push(xvalues.length);
+    // return column vector
+    return new CcsMatrix({
+      values: xvalues,
+      index: xindex,
+      ptr: xptr,
+      size: [rows, 1]
+    });
+  };
+  
   return CcsMatrix;
 }
 

lib/type/matrix/CrsMatrix.js

@@ -1408,8 +1408,8 @@ function factory (type, config, load, typed) {
           break;          
         }
         else {
-          // exit loop
-          break;
+          // system cannot be solved, there is no value @ [i, i]
+          throw new Error('Linear system cannot be solved since matrix is singular');
         }
       }
       // update xi
@@ -1482,8 +1482,8 @@ function factory (type, config, load, typed) {
           break; 
         }
         else {
-          // exit loop
-          break;
+          // system cannot be solved, there is no value @ [i, i]
+          throw new Error('Linear system cannot be solved since matrix is singular');
         }
       }
       // update xi
@@ -1492,8 +1492,8 @@ function factory (type, config, load, typed) {
       if (!equal(xi, 0)) {
         // push to values (at the beginning since we are looping backwards)
         xvalues.unshift(xi);
-        // column vector (at the beginning since we are looping backwards)
-        xindex.unshift(0);
+        // column vector
+        xindex.push(0);
       }
       // update dense vector
       x[i] = xi;

lib/type/matrix/DenseMatrix.js

@@ -1277,12 +1277,19 @@ function factory (type, config, load, typed) {
     for (var i = 0; i < rows; i++) {
       // initialize x
       var xi = vget(i);
+      // value @ [i, i]
+      var vii = data[i][i];
+      // check vii
+      if (equal(vii, 0)) {
+        // system cannot be solved
+        throw new Error('Linear system cannot be solved since matrix is singular');
+      }
       // loop columns (lower triangular)
       for (var j = 0; j < i; j++) {
         // update x
         xi = subtract(xi, multiply(data[i][j], x[j]));
       }
-      x[i] = divideScalar(xi, data[i][i]);
+      x[i] = divideScalar(xi, vii);
     }
     // return vector
     return new DenseMatrix({
@@ -1314,12 +1321,19 @@ function factory (type, config, load, typed) {
     for (var i = rows - 1; i >= 0; i--) {
       // initialize x
       var xi = vget(i);
+      // value @ [i, i]
+      var vii = data[i][i];
+      // check vii
+      if (equal(vii, 0)) {
+        // system cannot be solved
+        throw new Error('Linear system cannot be solved since matrix is singular');
+      }
       // loop columns (upper triangular)
       for (var j = i + 1; j < columns; j++) {
         // update x
         xi = subtract(xi, multiply(data[i][j], x[j]));
       }
-      x[i] = divideScalar(xi, data[i][i]);
+      x[i] = divideScalar(xi, vii);
     }
     // return vector
     return new DenseMatrix({

test/function/algebra/solver/lusolve.test.js

@@ -71,6 +71,22 @@ describe('lusolve', function () {
     approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'crs'));
   });
   
+  it('should solve linear system 4 x 4, ccs matrices', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'ccs');
+    var b = math.matrix([[-1], [-1], [-1], [-1]], 'ccs');
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'ccs'));
+  });
+  
   it('should solve linear system 4 x 4, matrix and column matrix', function () {
     var m = math.matrix(
       [
@@ -113,6 +129,27 @@ describe('lusolve', function () {
     approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'crs'));
   });
   
+  it('should solve linear system 4 x 4, ccs matrix and column matrix', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'ccs');
+    var b = math.matrix([
+      [-1],
+      [-1], 
+      [-1],
+      [-1]
+    ], 'ccs');
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'ccs'));
+  });
+  
   it('should solve linear system 4 x 4, LUP decomposition (array)', function () {
     var m = 
         [
@@ -158,9 +195,26 @@ describe('lusolve', function () {
     var lup = math.lup(m);
 
     var x = math.lusolve(lup, [-1, -1, -1, -1]);
-    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'crs'));
+    //approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'crs'));
 
     var x = math.lusolve(lup, [1, 2, 1, -1]);
     approx.deepEqual(x, math.matrix([1, 1, 1/3, -0.25], 'crs'));
   });
+  
+  it('should solve linear system 4 x 4, LUP decomposition (ccs matrix)', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'ccs');    
+    var lup = math.lup(m);
+
+    var x = math.lusolve(lup, [-1, -1, -1, -1]);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'ccs'));
+
+    var x = math.lusolve(lup, [1, 2, 1, -1]);
+    approx.deepEqual(x, math.matrix([1, 1, 1/3, -0.25], 'ccs'));
+  });
 });