From 1a68545416ba52cdd6b9531a44de333d7c93b944 Mon Sep 17 00:00:00 2001 From: Peihong Guo Date: Thu, 7 Nov 2013 10:15:33 -0600 Subject: [PATCH 1/5] Added conjugate gradient solver. Supports both full matrix and sparse matrix. --- src/iterative.js | 68 ++++++++++++++++++++++++++++++++++++++++++++++++ src/numeric.js | 22 ++++++++++++++++ 2 files changed, 90 insertions(+) create mode 100644 src/iterative.js diff --git a/src/iterative.js b/src/iterative.js new file mode 100644 index 0000000..be7f734 --- /dev/null +++ b/src/iterative.js @@ -0,0 +1,68 @@ +numeric.ConjugateGradient = function ConjugateGradient(A, b, maxIters, residue) { + var maxIters = maxIters || 1024; + var residue = residue || 1e-6; + + var CG = function( A, b, dotVV, dotMV ) { + // initialization + var iters = 0; + var converged = false; + var rows = b.length; + + var dot = dotVV; + var mv = dotMV; + var add = numeric.add; + var mul = numeric.mul; + var axpy = function(alpha, x, y) { + return add(mul(x, alpha), y); + }; + + // initialize x + var x = numeric.rep([rows], 0); + + // r = b - A x + // x is zero initially + var r = mul(b, 1); + var p = mul(b, 1); + + while( !converged && iters < maxIters ) { + var rDotr = dot(r, r); + var Ap = mv(A, p); + var pAp = dot(p, Ap); + var alpha = rDotr / pAp; + + x = axpy(alpha, p, x); + r = axpy(-alpha, Ap, r); + + var rDotr_new = dot(r, r); + + if( rDotr_new < residue ) { + converged = true; + break; + } + + var beta = rDotr_new / rDotr; + p = axpy(beta, p, r); + + iters++; + } + + //console.log('converged in ' + iters + ' iterations'); + + return x; + } + + switch(A.format) { + case 'full': { + return CG(A, b, numeric.dot, numeric.dot); + } + case 'ccs': { + return CG(A, b, numeric.dot, numeric.ccsMV); + } + case 'crs': { + return CG(A, b, numeric.dot, numeric.cdotMV); + } + default: { + throw 'Not supported matrix format'; + } + } +}; \ No newline at end of file diff --git a/src/numeric.js b/src/numeric.js index d401f61..be0e32f 100644 --- a/src/numeric.js +++ b/src/numeric.js @@ -1764,6 +1764,28 @@ numeric.ccsGetBlock = function ccsGetBlock(A,i,j) { return B; } +// CCS matrix mul dense vector +numeric.ccsMV = function ccsMV(A, x) { + var Ai = A[0], Aj = A[1], Av = A[2]; + var sA = numeric.ccsDim(A); + var m = sA[0], n = sA[1]; + var L = x.length; + var ret = numeric.rep([L],0); + if( n !== L ) throw 'Matrix dimension does not match input vector.'; + var i, j, k, j0, j1; + var ri, val; + for(k=0;k!==n;k++) { + j0 = Ai[k]; + j1 = Ai[k+1]; + for(j=j0;j Date: Thu, 7 Nov 2013 11:08:37 -0600 Subject: [PATCH 2/5] added bicgstab algorithm --- src/iterative.js | 101 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 93 insertions(+), 8 deletions(-) diff --git a/src/iterative.js b/src/iterative.js index be7f734..c0ca2b6 100644 --- a/src/iterative.js +++ b/src/iterative.js @@ -1,4 +1,84 @@ -numeric.ConjugateGradient = function ConjugateGradient(A, b, maxIters, residue) { +numeric.bicgstab = function bicgstab(A, b, maxIters, residue) { + var maxIters = maxIters || 1024; + var residue = residue || 1e-6; + + var BiCG = function( A, b, dotVV, dotMV ) { + // initialization + var iters = 0; + var converged = false; + var rows = b.length; + + var dot = dotVV; + var mv = dotMV; + var add = numeric.add; + var mul = numeric.mul; + var axpy = function(alpha, x, y) { + return add(mul(x, alpha), y); + }; + + // initialize x + var x = numeric.rep([rows], 0); + + // r = b - A x + // x is zero initially + var r = mul(b, 1); + var rhat = mul(b, 1); + var rho = 1, alpha = 1, w = 1; + var p = mul(b, 1); + var v = numeric.rep([rows], 0); + var s; + + var bnorm = dot(b, b); + + while( !converged && iters < maxIters ) { + var rDotr = dot(rhat, r); + + if( (rDotr / bnorm) <= residue ) { + converged = true; + break; + } + + var v = mv(A, p); + var alpha = rDotr / dot(rhat, v); + s = axpy(-alpha, v, r); + + var t = mv(A, s); + w = dot(t, s) / dot(t, t); + + x = axpy(w, s, axpy(alpha, p, x)); + r = axpy(-w, t, s); + + var beta = dot(rhat, r) / rDotr * (alpha/w); + p = axpy(beta, axpy(-w, v, p), r); + + iters++; + } + + if( converged ) + console.log('converged in ' + iters + ' iterations'); + else + console.log('not converged in ' + iters + ' iterations'); + + return x; + } + + switch(A.format) { + case 'full': { + return BiCG(A, b, numeric.dot, numeric.dot); + } + case 'ccs': { + return BiCG(A, b, numeric.dot, numeric.ccsMV); + } + case 'crs': { + return BiCG(A, b, numeric.dot, numeric.cdotMV); + } + default: { + throw 'Not supported matrix format'; + } + } +} + +numeric.cg = function cg(A, b, maxIters, residue) { var maxIters = maxIters || 1024; var residue = residue || 1e-6; @@ -24,8 +104,16 @@ numeric.ConjugateGradient = function ConjugateGradient(A, b, maxIters, residue) var r = mul(b, 1); var p = mul(b, 1); + var bnorm = dot(b, b); + while( !converged && iters < maxIters ) { var rDotr = dot(r, r); + + if( (rDotr / bnorm) <= residue ) { + converged = true; + break; + } + var Ap = mv(A, p); var pAp = dot(p, Ap); var alpha = rDotr / pAp; @@ -34,19 +122,16 @@ numeric.ConjugateGradient = function ConjugateGradient(A, b, maxIters, residue) r = axpy(-alpha, Ap, r); var rDotr_new = dot(r, r); - - if( rDotr_new < residue ) { - converged = true; - break; - } - var beta = rDotr_new / rDotr; p = axpy(beta, p, r); iters++; } - //console.log('converged in ' + iters + ' iterations'); + if( converged ) + console.log('converged in ' + iters + ' iterations'); + else + console.log('not converged in ' + iters + ' iterations'); return x; } From a1284cdf98ef18170b81e989639f6779284e3dae Mon Sep 17 00:00:00 2001 From: Peihong Guo Date: Thu, 7 Nov 2013 15:27:50 -0600 Subject: [PATCH 3/5] added test cases for iterative solvers. --- benchmark3.html | 196 +++++++++++++++++++++++++++++++++++++++++++++++ src/iterative.js | 54 +++++++++---- 2 files changed, 234 insertions(+), 16 deletions(-) create mode 100644 benchmark3.html diff --git a/benchmark3.html b/benchmark3.html new file mode 100644 index 0000000..2646f01 --- /dev/null +++ b/benchmark3.html @@ -0,0 +1,196 @@ + + + + +Numeric Javascript: Benchmarks + + + + +We are now running a linear algebra performance benchmark in your browser! Please ensure that your seatbelt +is fastened and your tray table is upright while we invert 100x100 matrices.

+ +Performance (MFLOPS) +
+
+
+
+
Geometric mean of scores:

+ +Higher is better: For each benchmark and library, a function is called repeatedly for a certain amount of time. +The number of function calls per seconds is converted into a FLOP rate. As we move right within each test, the matrix size increases.

+ +What tricks are used to increase performance in Numeric? +
    +
  • Replace inner loop for(j=0;j<n;j++) A[i][j] by the equivalent Ai = A[i]; for(j=0;j<n;j++) Ai[j] ("hoisting" the [i] out of the loop). +
  • Preallocate Arrays: A = new Array(n) instead of A = []. +
  • Use Array objects directly (abstractions slow you down). Getters and setters are bad. +
  • Use for(j=n-1;j>=0;j--) if it is faster. +
  • Do not put anything in Array.prototype. If you modify Array.prototype, it slows down everything significantly. +
  • Big Matrix*Matrix product: transpose the second matrix and rearrange the loops to exploit locality. +
  • Unroll loops. +
  • Don't nest functions. +
  • Don't call functions, inline them manually. Except... +
  • ...big functions can confuse the JIT. If a new code path is run in a function, the function can be deoptimized by the JIT. +
  • Avoid polymorphism. +
  • Localize references. For example: replace for(i=n-1;i>=0;i--) x[i] = Math.random(); by rnd = Math.random; for(i=n-1;i>=0;i--) x[i] = rnd();. (Make sure rnd and x really are local variables!) +
  • Deep lexical scopes are bad. You can create a function without much of a lexical scope by using + new Function('...');. +
+
+ +GC pauses? +If you reload the page, the benchmark will run again and will give slightly different results. +This could be due to GC pauses or other background tasks, DOM updates, etc... +To get an idea of the impact of this, load this page in two or more different tabs (not at the same time, +one after the other) and switch between the tabs and see the differences in the performance chart. +


+ +
+ + + + + + + + +


+ diff --git a/src/iterative.js b/src/iterative.js index c0ca2b6..12f6939 100644 --- a/src/iterative.js +++ b/src/iterative.js @@ -1,3 +1,25 @@ +// CCS matrix mul dense vector +numeric.ccsMV = numeric.ccsMV || function ccsMV(A, x) { + var Ai = A[0], Aj = A[1], Av = A[2]; + var sA = numeric.ccsDim(A); + var m = sA[0], n = sA[1]; + var L = x.length; + var ret = numeric.rep([L],0); + if( n !== L ) throw 'Matrix dimension does not match input vector.'; + var i, j, k, j0, j1; + var ri, val; + for(k=0;k!==n;k++) { + j0 = Ai[k]; + j1 = Ai[k+1]; + for(j=j0;j Date: Thu, 7 Nov 2013 16:38:38 -0600 Subject: [PATCH 4/5] 1. minor change in the solvers 2. bug fixed in the benchmark page --- benchmark3.html | 2 +- src/iterative.js | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmark3.html b/benchmark3.html index 2646f01..c0021f6 100644 --- a/benchmark3.html +++ b/benchmark3.html @@ -85,7 +85,7 @@ ], [ 'Solve Ax=b, A is Banded Laplacian', [4, 6, 9, 14, 21, 32, 48, 72], - function(n) { var A = numeric.cdelsq(numeric.cgrid(n)); A.format = 'crs'; var b = mkV(n); return bench(function() { numeric.ccsLUPSolve(numeric.ccsLUP(A), b); }); }, + function(n) { var A = numeric.cdelsq(numeric.cgrid(n)); A.format = 'crs'; var b = mkV(n); return bench(function() { numeric.cLUsolve(numeric.cLU(A), b); }); }, function(n) { var A = numeric.cdelsq(numeric.cgrid(n)); A.format = 'crs'; var b = mkV(n); return bench(function() { numeric.cg(A, b); }); }, function(n) { var A = numeric.cdelsq(numeric.cgrid(n)); A.format = 'crs'; var b = mkV(n); return bench(function() { numeric.bicgstab(A, b); }); }, ], diff --git a/src/iterative.js b/src/iterative.js index 12f6939..0223b87 100644 --- a/src/iterative.js +++ b/src/iterative.js @@ -61,7 +61,7 @@ numeric.bicgstab = function bicgstab(A, b, maxIters, residue) { } v = mv(A, p); - var alpha = rDotr / dot(rhat, v); + alpha = rDotr / dot(rhat, v); s = axpy(-alpha, v, r); t = mv(A, s); @@ -70,7 +70,7 @@ numeric.bicgstab = function bicgstab(A, b, maxIters, residue) { x = axpy(w, s, axpy(alpha, p, x)); r = axpy(-w, t, s); - var beta = dot(rhat, r) / rDotr * (alpha/w); + beta = dot(rhat, r) / rDotr * (alpha/w); p = axpy(beta, axpy(-w, v, p), r); iters++; From b6789e4ecd788caaa60306a8f667a701b072227b Mon Sep 17 00:00:00 2001 From: Peihong Guo Date: Fri, 8 Nov 2013 23:25:18 -0600 Subject: [PATCH 5/5] added SOR solver. the performance of SOR solver is not as good as CG and BiCGSTAB, but similar to ccsLUPSolve. in practice, CG and BiCGSTAB should always be the best options. --- src/iterative.js | 207 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 206 insertions(+), 1 deletion(-) diff --git a/src/iterative.js b/src/iterative.js index 0223b87..5b5f955 100644 --- a/src/iterative.js +++ b/src/iterative.js @@ -172,4 +172,209 @@ numeric.cg = function cg(A, b, maxIters, residue) { throw 'Not supported matrix format'; } } -}; \ No newline at end of file +}; + +numeric.sor = function sor(A, b, relax, maxIters, residue) { + var maxIters = maxIters || 1024; + var residue = residue || 1e-6; + var relax = relax || 1.0; // no relaxation by default, fall back to Gauss-Seidel + + var sor_full = function(A, b) { + // initialization + var iters = 0; + var converged = false; + var rows = b.length; + + var sA = numeric.dim(A); + var n = sA[0], m = sA[1]; + if( n != rows ) { + throw 'Matrix dimension does not match input vector.'; + } + + // initialize x + var x = numeric.rep([rows], 0); + var dot = numeric.dot; + var bnorm = dot(b, b); + var rowsum = 0; + var i, j; + var Ai; + var rowdiff; + var res = numeric.rep([rows], 0), r2; + + while(!converged && iters < maxIters) { + for(i=0;i