mirror of
https://github.com/josdejong/mathjs.git
synced 2026-02-01 16:07:46 +00:00
lup() - partial
This commit is contained in:
rjbaucells
parent
314b705ef6
commit
a478198e12
@ -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<Matrix>} 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++) {
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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));
|
||||
});
|
||||
});
|
||||
});
|
||||
Loading…
x
Reference in New Issue
Block a user