From a478198e1274d0f9cf9f3fb248dc547a64694932 Mon Sep 17 00:00:00 2001 From: rjbaucells Date: Sun, 5 Apr 2015 21:08:46 -0400 Subject: [PATCH] lup() - partial --- lib/type/matrix/CcsMatrix.js | 91 ++++++++++++++++++++---------- lib/type/matrix/Spa.js | 75 ++++++++++++++++-------- test/type/matrix/CcsMatrix.test.js | 39 ++++++------- 3 files changed, 131 insertions(+), 74 deletions(-) diff --git a/lib/type/matrix/CcsMatrix.js b/lib/type/matrix/CcsMatrix.js index 5dd0d3560..a0436bd8d 100644 --- a/lib/type/matrix/CcsMatrix.js +++ b/lib/type/matrix/CcsMatrix.js @@ -1232,7 +1232,7 @@ function factory (type, config, load, typed) { /** * LU decomposition with partial pivoting. Matrix A is decomposed in three matrices (L, U, P) where: - * P * A = L * U + * A * P = L * U * * @return {Array} The lower triangular matrix, the upper triangular matrix and the permutation matrix. */ @@ -1267,7 +1267,7 @@ function factory (type, config, load, typed) { // loop columns for (j = 0; j < columns; j++) { // sparse accumulator - var spa = new Spa(); + var spa = new Spa(rows); // check lower triangular matrix has a value @ column j if (j < rows) { // update ptr @@ -1290,36 +1290,54 @@ function factory (type, config, load, typed) { } // skip first column in upper triangular matrix if (j > 0) { - // loop k within [k0, k1[ - for (k = k0; k < k1; k++) { - // row (use permutation vector) - i = index[k]; - // check row is in the upper triangular - if (i < j) { - // value in spa @ i - var vij = spa.get(i); - // value indexes from column i in matrix L - var kl0 = lptr[i]; - var kl1 = lptr[i + 1]; - // loop values in column i - for (var kl = kl0; kl < kl1; kl++) { - // row - var r = pvector[lindex[kl]]; - // check row is below i - if (r > i) { - // update spa @ r (v[r,j] -= L[r,i] * v[i,j]) - spa.accumulate(r, unaryMinus(multiply(lvalues[kl], vij))); - } + // loop rows in column j (above diagonal) + /*spa.forEach(0, j, function (x, vxj) { + // loop rows in column x (L) + _forEachRow(x, lvalues, lindex, lptr, function (i, vix) { + // check row is below x + if (i > x) { + // update spa value + spa.accumulate(i, unaryMinus(multiply(vix, vxj))); } + }); + });*/ + + + // loop kj within [k0, k1[ + /*for (var kj = k0; kj < k1; kj++) { + // row (use permutation vector) + i = pvector[index[kj]]; + // min i,j + var min = Math.min(i, j); + // v[i, j] + var s = 0; + // loop up to min + for (k = 0; k < min; k++) { + // value indexes in column k + var ki0 = lptr[k]; + var ki1 = lptr[k + 1]; + // row + var r = -1; + // loop values in column k + for (var ki = ki0; ki < ki1 && r < i; ki++) { + // row + r = lindex[ki]; + // only process row i + if (r === i) { + // s = l[i, k] - data[k, j] + s = add(s, multiply(lvalues[ki], spa.get(r))); + } + } } - } + spa.accumulate(i, unaryMinus(s)); + } */ } // 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) - spa.forEach(function (x, v) { + spa.forEach(j, rows - 1, function (x, v) { // absolute value var absv = abs(v); // value is greater than pivote value @@ -1331,11 +1349,13 @@ function factory (type, config, load, typed) { // value @ [j, j] vjj = v; } - }, j); + }); // swap rows (j <-> pi) if (j !== pi) { // swap values j <-> pi in L - // _swapRows(j, pi, lsize[1], lvalues, lindex, lptr); + _swapRows(j, pi, lsize[1], lvalues, lindex, lptr); + // swap values j <-> pi in U + _swapRows(j, pi, usize[1], uvalues, uindex, uptr); // swap values j <-> pi in P _swapRows(j, pi, p._size[1], p._values, p._index, p._ptr); // swap values in spa @@ -1344,16 +1364,16 @@ function factory (type, config, load, typed) { pvector[j] = [pvector[pi], pvector[pi] = pvector[j]][0]; } // loop values in spa (order by row, above diagonal, include j) - spa.forEach(function (x, v) { + spa.forEach(0, j, function (x, v) { // check value is non zero if (!equal(v, 0)) { // update upper triangular matrix uvalues.push(v); uindex.push(x); } - }, 0, j); + }); // loop values in spa (order by row, below diagonal) - spa.forEach(function (x, v) { + spa.forEach(j + 1, rows - 1, function (x, v) { // update value v = divideScalar(v, vjj); // check value is non zero @@ -1362,7 +1382,7 @@ function factory (type, config, load, typed) { lvalues.push(v); lindex.push(x); } - }, j + 1); + }); } // update ptrs lptr.push(lvalues.length); @@ -1392,6 +1412,17 @@ function factory (type, config, load, typed) { }; }; + var _forEachRow = function (j, values, index, ptr, callback) { + // indeces for column j + var k0 = ptr[j]; + var k1 = ptr[j + 1]; + // loop + for (var k = k0; k < k1; k++) { + // invoke callback + callback(index[k], values[k]); + } + }; + var _swapRows = function (x, y, columns, values, index, ptr) { // loop columns for (var j = 0; j < columns; j++) { diff --git a/lib/type/matrix/Spa.js b/lib/type/matrix/Spa.js index 118c0118a..32b3276b6 100644 --- a/lib/type/matrix/Spa.js +++ b/lib/type/matrix/Spa.js @@ -4,15 +4,21 @@ function factory (type, config, load, typed) { var add = load(require('../../function/arithmetic/add')); - function Spa() { + function Spa(length) { if (!(this instanceof Spa)) throw new SyntaxError('Constructor must be called with the new operator'); // allocate vector, TODO use typed arrays - this._values = []; + this._values = new Array(length); + this._index = []; } Spa.prototype.set = function (i, v) { + // check we have a value @ i + if (!this._values[i]) { + // store index + this._index.push(i); + } // set the value @ i this._values[i] = v; }; @@ -22,36 +28,59 @@ function factory (type, config, load, typed) { }; Spa.prototype.accumulate = function (i, v) { + // current value + var value = this._values[i]; + // check we have a value @ i + if (!value) { + // store index + this._index.push(i); + // initialize value + value = 0; + } // accumulate value - this._values[i] = add(this._values[i] || 0, v); + this._values[i] = add(value, v); }; - 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, values = this._values; i <= to; i++) { - // value @ i - var v = values[i]; - // invoke callback if needed - if (v !== undefined) - callback(i, v, this); + Spa.prototype.forEach = function (from, to, callback) { + // loop indexes + for (var i = 0, values = this._values, index = this._index; i < index.length; i++) { + // values index + var k = index[i]; + // check it is in range + if (k >= from && k <=to) { + // value @ k + var v = values[k]; + // invoke callback if needed + if (v !== undefined) + callback(k, v, this); + } } }; Spa.prototype.swap = function (i, j) { - // swap values i <-> j + // values @ i and j var vi = this._values[i]; - this._values[i] = this._values[j]; + var vj = this._values[j]; + // swap values + this._values[i] = vj; this._values[j] = vi; + // chek we need to insert indeces + if (!vi && vj) { + // insert i + this._index.push(i); + // remove j + var k = this._index.indexOf(j); + if (k >=0) + this._index.splice(k, 1); + } + else if (vi && !vj) { + // insert j + this._index.push(j); + // remove i + var k = this._index.indexOf(i); + if (k >=0) + this._index.splice(k, 1); + } }; return Spa; diff --git a/test/type/matrix/CcsMatrix.test.js b/test/type/matrix/CcsMatrix.test.js index 5a449ad5e..0680b71ef 100644 --- a/test/type/matrix/CcsMatrix.test.js +++ b/test/type/matrix/CcsMatrix.test.js @@ -2074,7 +2074,7 @@ describe('CcsMatrix', function() { assert.deepEqual(math.multiply(r.P, m), math.multiply(r.L, r.U)); }); - it.skip('should decompose matrix, n x n', function () { + it('should decompose matrix, n x n', function () { var m = new CcsMatrix( [ [16, -120, 240, -140], @@ -2085,33 +2085,28 @@ describe('CcsMatrix', function() { ); var r = m.lup(); - //assert.deepEqual(math.multiply(r.P, m), math.multiply(r.L, r.U)); // L assert.deepEqual( - r.L, - new CcsMatrix( - [ - [1, 0, 0, 0], - [-0.5, 1, 0, 0], - [-0.5833333, -0.7, 1, 0], - [0.0666667, -0.4, -0.5714286, 1] - ] - )); + r.L.valueOf(), + [ + [1, 0, 0, 0], + [-0.5, 1, 0, 0], + [-0.5833333333333334, -0.7, 1, 0], + [0.06666666666666667, -0.4, -0.5714285714285714, 1] + ]); // U assert.deepEqual( - r.U, - new CcsMatrix( - [ - [240, -2700, 6480, -4200], - [0, -150, 540, -420], - [0, 0, -42, 56], - [0, 0, 0, 4] - ] - )); + r.U.valueOf(), + [ + [240, -2700, 6480, -4200], + [0, -150, 540, -420], + [0, 0, -42, 56], + [0, 0, 0, 4] + ]); // P assert.deepEqual( r.P, - new CcsMatrix( + new DenseMatrix( [ [0, 0, 1, 0], [0, 1, 0, 0], @@ -2119,6 +2114,8 @@ describe('CcsMatrix', function() { [1, 0, 0, 0] ] )); + // verify + assert.deepEqual(math.multiply(r.P, m), math.multiply(r.L, r.U)); }); }); }); \ No newline at end of file