From ccaaaedbb9b68d9be417dceedc17428cabace853 Mon Sep 17 00:00:00 2001 From: rjbaucells Date: Wed, 8 Apr 2015 23:48:19 -0400 Subject: [PATCH] lusolve() - ccs --- lib/function/algebra/solver/lusolve.js | 3 + lib/type/matrix/CcsMatrix.js | 143 +++++++++++++++++++ lib/type/matrix/CrsMatrix.js | 12 +- lib/type/matrix/DenseMatrix.js | 18 ++- test/function/algebra/solver/lusolve.test.js | 56 +++++++- 5 files changed, 223 insertions(+), 9 deletions(-) diff --git a/lib/function/algebra/solver/lusolve.js b/lib/function/algebra/solver/lusolve.js index 3b81a24b7..7a30bd336 100644 --- a/lib/function/algebra/solver/lusolve.js +++ b/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 diff --git a/lib/type/matrix/CcsMatrix.js b/lib/type/matrix/CcsMatrix.js index e69fe00b6..63d1874f1 100644 --- a/lib/type/matrix/CcsMatrix.js +++ b/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; } diff --git a/lib/type/matrix/CrsMatrix.js b/lib/type/matrix/CrsMatrix.js index 23e9c26b4..a1731c5f3 100644 --- a/lib/type/matrix/CrsMatrix.js +++ b/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; diff --git a/lib/type/matrix/DenseMatrix.js b/lib/type/matrix/DenseMatrix.js index 3a1474ed9..d3a194e07 100644 --- a/lib/type/matrix/DenseMatrix.js +++ b/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({ diff --git a/test/function/algebra/solver/lusolve.test.js b/test/function/algebra/solver/lusolve.test.js index 4e588806b..a5adaf6d2 100644 --- a/test/function/algebra/solver/lusolve.test.js +++ b/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')); + }); }); \ No newline at end of file