diff --git a/lib/type/matrix/CcsMatrix.js b/lib/type/matrix/CcsMatrix.js index ccf20bf90..8943f129e 100644 --- a/lib/type/matrix/CcsMatrix.js +++ b/lib/type/matrix/CcsMatrix.js @@ -18,7 +18,6 @@ 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')); @@ -1231,18 +1230,19 @@ function factory (type, config, load, typed) { }); }; + /** + * LU decomposition with partial pivoting. + */ CcsMatrix.prototype.lup = function () { // rows & columns var rows = this._size[0]; var columns = this._size[1]; - // matrix arrays + // matrix arrays (will not be modified) var values = this._values; var index = this._index; var ptr = this._ptr; - // diagonal length - var n = rows > columns ? columns : rows; // p (eye [n, n]) - var p = CcsMatrix.diagonal([n, n], 1, 0); + var p = CcsMatrix.diagonal([rows, rows], 1, 0); // l matrix arrays var lvalues = []; var lindex = []; @@ -1251,21 +1251,22 @@ function factory (type, config, load, typed) { var uvalues = []; var uindex = []; var uptr = []; - - debugger; - // vars var i, j, k; // loop columns for (j = 0; j < columns; j++) { // sparse accumulator - var spa = new Spa(rows); - // update ptrs - lptr.push(lvalues.length); + var spa = new Spa(); + // check lower diagonal matrix has a value @ column j + if (j < rows) { + // update ptr + lptr.push(lvalues.length); + // first value in j column for lower diagonal matrix + lvalues.push(1); + lindex.push(j); + } + // update ptr uptr.push(uvalues.length); - // first value in j column for lower diagonal matrix - lvalues.push(1); - lindex.push(j); // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1] var k0 = ptr[j]; var k1 = ptr[j + 1]; @@ -1286,24 +1287,24 @@ function factory (type, config, load, typed) { if (i < j) { // value in spa @ i var vij = spa.get(i); - // value indexes from column i - var ki0 = ptr[i]; - var ki1 = ptr[i + 1]; + // value indexes from column i in matrix L + var kl0 = lptr[i]; + var kl1 = lptr[i + 1]; // loop values in column i - for (var ki = ki0; ki < ki1; ki++) { + for (var kl = kl0; kl < kl1; kl++) { // row - var r = index[ki]; + var r = lindex[kl]; // check row is below i if (r > i) { - // update spa @ r - spa.accumulate(r, unaryMinus(multiply(values[ki], vij))); + // update spa @ r (v[r,j] -= L[r,i] * v[i,j]) + spa.accumulate(r, unaryMinus(multiply(lvalues[kl], vij))); } } } } } - // row with larger value in matrix @ column j, row >= j - var pi = 0; + // row with larger value in spa, row >= j + var pi = j; var pabsv = 0; var vjj = 0; // loop values in spa (order by row, below diagonal) @@ -1323,7 +1324,7 @@ function factory (type, config, load, typed) { // swap rows (j <-> pi) if (j !== pi) { } - // loop values in spa (order by row, above diagonal) + // loop values in spa (order by row, above diagonal, include j) spa.forEach(function (x, v) { // check value is non zero if (!equal(v, 0)) { @@ -1352,14 +1353,14 @@ function factory (type, config, load, typed) { values: lvalues, index: lindex, ptr: lptr, - size: [rows, n] + size: [rows, Math.min(rows, columns)] }); // u matrix var u = new CcsMatrix({ values: uvalues, index: uindex, ptr: uptr, - size: [n, columns] + size: [Math.min(rows, columns), columns] }); // return matrices return { diff --git a/lib/type/matrix/Spa.js b/lib/type/matrix/Spa.js index 153591d18..43e553efa 100644 --- a/lib/type/matrix/Spa.js +++ b/lib/type/matrix/Spa.js @@ -4,7 +4,7 @@ function factory (type, config, load, typed) { var add = load(require('../../function/arithmetic/add')); - function Spa(length) { + function Spa() { if (!(this instanceof Spa)) throw new SyntaxError('Constructor must be called with the new operator'); @@ -27,8 +27,18 @@ function factory (type, config, load, typed) { }; Spa.prototype.forEach = function (callback, from, to) { + // check arguments + if (arguments.length === 1) { + // from & to + from = 0; + to = this._values.length - 1; + } + else if (arguments.length === 2) { + // to + to = this._values.length - 1; + } // loop indexes [from, to] - for (var i = from || 0, values = this._values, to = to || values.length - 1; i <= to; i++) { + for (var i = from, values = this._values; i <= to; i++) { // value @ i var v = values[i]; // invoke callback if needed diff --git a/test/type/matrix/CcsMatrix.test.js b/test/type/matrix/CcsMatrix.test.js index ad5f127b3..5ebd708e9 100644 --- a/test/type/matrix/CcsMatrix.test.js +++ b/test/type/matrix/CcsMatrix.test.js @@ -1966,7 +1966,7 @@ describe('CcsMatrix', function() { new CcsMatrix( [ [1, 0], - [1.5, 1] + [0.5, 1] ] )); assert.deepEqual( @@ -1986,5 +1986,80 @@ describe('CcsMatrix', function() { ] )); }); + + it('should decompose matrix, m x n, m < n, no pivoting', function () { + var m = new CcsMatrix( + [ + [2, 1, 1], + [1, 4, 5] + ] + ); + + var r = m.lup(); + + assert.deepEqual( + r.L, + new CcsMatrix( + [ + [1, 0], + [0.5, 1] + ] + )); + assert.deepEqual( + r.U, + new CcsMatrix( + [ + [2, 1, 1], + [0, 3.5, 4.5] + ] + )); + assert.deepEqual( + r.P, + new CcsMatrix( + [ + [1, 0], + [0, 1] + ] + )); + }); + + it('should decompose matrix, m x n, m > n, no pivoting', function () { + var m = new CcsMatrix( + [ + [8, 2], + [6, 4], + [4, 1] + ] + ); + + var r = m.lup(); + + assert.deepEqual( + r.L, + new CcsMatrix( + [ + [1, 0], + [0.75, 1], + [0.5, 0] + ] + )); + assert.deepEqual( + r.U, + new CcsMatrix( + [ + [8, 2], + [0, 2.5] + ] + )); + assert.deepEqual( + r.P, + new CcsMatrix( + [ + [1, 0, 0], + [0, 1, 0], + [0, 0, 1] + ] + )); + }); }); }); \ No newline at end of file