mirror of
https://github.com/josdejong/mathjs.git
synced 2026-01-18 14:59:29 +00:00
CCS lup - partial
This commit is contained in:
rjbaucells
parent
c2de167a17
commit
2cef2ae244
@ -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 {
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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]
|
||||
]
|
||||
));
|
||||
});
|
||||
});
|
||||
});
|
||||
Loading…
x
Reference in New Issue
Block a user