Backup_Repos/jsketcher
1
0
Fork 0
mirror of https://github.com/xibyte/jsketcher synced 2025-12-10 18:36:30 +01:00
Code Issues Projects Releases Wiki Activity

conjugate gradient for linear systems

Browse source
This commit is contained in:
Val Erastov 2014-11-01 00:09:25 -07:00
parent acacd8626e
commit 178ea85934
3 changed files with 44 additions and 4 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

39
web/app/math/noptim.js
View file

@ -205,6 +205,7 @@ optim.inv = function inv(x) {
Ij = I[i0]; I[i0] = I[j]; I[j] = Ij;
x = Aj[j];
if (x === 0) {
console.log("CAN' INVERSE MATRIX");
x = 1e-32
}
for(k=j;k!==n;++k) Aj[k] /= x;
@ -338,7 +339,9 @@ optim.dog_leg = function(subsys) {
// get the gauss-newton step
// h_gn = n.solve(Jx, n.mul(fx, -1));
h_gn = lsolve(Jx, n.mul(fx, -1));
//h_gn = lsolve(Jx, n.mul(fx, -1));
h_gn = optim.cg(Jx, h_gn, n.mul(fx, -1), 1e-8, maxIterNumber);
// solve linear problem using svd formula to get the gauss-newton step
// h_gn = lls(Jx, n.mul(fx, -1));
@ -436,3 +439,37 @@ optim.dog_leg = function(subsys) {
};
optim.cg = function(A, x, b, tol, maxIt) {
var _ = numeric;
var tr = _.transpose;
var At = tr(A);
if (A.length != A[0].length) {
var A = _.dot(At, A);
var b = _.dot(At, b);
At = tr(A);
}
var r = _.sub(_.dot(A, x), b);
var p = _.mul(r, -1);
var rr = _.dotVV(r, r);
var a;
var _rr;
var beta;
for (var i = 0; i < maxIt; ++i) {
if (_.norm2(r) <= tol) break;
var Axp =_.dot(A, p);
a = rr / _.dotVV(Axp, p);
x = _.add(x, _.mul(p, a));
r = _.add(r, _.mul(Axp, a));
_rr = rr;
rr = _.dotVV(r, r);
beta = rr / _rr;
p = _.add(_.mul(r, -1), _.mul(p, beta));
}
// console.log("liner problem solved in " + i);
return x;
};

7
web/app/math/test.js
View file

@ -57,8 +57,11 @@ function lsolve() {
var n = numeric;
A = [[1, 0, 0], [2,3,0]]
b = [10,25]
var At = n.transpose(A);
var res = n.dot(n.dot(At, n.inv(n.dot(A, At)) ), b);
// var At = n.transpose(A);
// var A = n.dot(At, A);
// var b = n.dot(At, b);
// var res = n.dot(n.dot(At, n.inv(n.dot(A, At)) ), b);
var res = optim.cg(A, [100, 100, 0], b, 1e-8, 800);
console.log(res);
}

2
web/test.html
View file

@ -31,7 +31,7 @@
new TCAD.App2D();
}
window.onload = function() {
testCompare();
// testCompare();
lsolve();
}
</script>