diff --git a/benchmark3.html b/benchmark3.html
new file mode 100644
index 0000000..c0021f6
--- /dev/null
+++ b/benchmark3.html
@@ -0,0 +1,196 @@
+<html>
+<head>
+<link rel="SHORTCUT ICON" href="favicon.ico">
+<link rel="stylesheet" type="text/css" href="resources/style.css">
+<title>Numeric Javascript: Benchmarks</title>
+</head>
+<body>
+<!--#include file="resources/header.html" -->
+
+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.<br><br>
+
+<b>Performance (<a href="http://en.wikipedia.org/wiki/FLOPS">MFLOPS</a>)</b>
+<div style="width:1000px;overflow:hidden;font-size:14px;">
+<div id="placeholder" style="width:700px;height:500px;float:left;"></div>
+<div id="legend" style="width:250px;height:100px;overflow:hidden;"></div>
+</div>
+<div id="meanscore">Geometric mean of scores: </div><br>
+
+<b>Higher is better:</b> 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.<br><br>
+
+<b>What tricks are used to increase performance in Numeric?</b>
+<ul>
+  <li>Replace inner loop <tt>for(j=0;j&lt;n;j++) A[i][j]</tt> by the equivalent <tt>Ai = A[i]; for(j=0;j&lt;n;j++) Ai[j]</tt> ("hoisting" the <tt>[i]</tt> out of the loop</tt>).
+  <li>Preallocate Arrays: <tt>A = new Array(n)</tt> instead of <tt>A = []</tt>.
+  <li>Use <tt>Array</tt> objects directly (abstractions slow you down). Getters and setters are bad.
+  <li>Use <tt>for(j=n-1;j&gt;=0;j--)</tt> if it is faster.
+  <li>Do not put anything in <tt>Array.prototype</tt>. If you modify <tt>Array.prototype</tt>, it slows down everything significantly.
+  <li>Big Matrix*Matrix product: transpose the second matrix and rearrange the loops to exploit locality.
+  <li>Unroll loops.
+  <li>Don't nest functions.
+  <li>Don't call functions, inline them manually. Except...
+  <li>...big functions can confuse the JIT. If a new code path is run in a function, the function can be deoptimized by the JIT.
+  <li>Avoid polymorphism.
+  <li>Localize references. For example: replace <tt>for(i=n-1;i>=0;i--) x[i] = Math.random();</tt> by <tt>rnd = Math.random; for(i=n-1;i>=0;i--) x[i] = rnd();</tt>. (Make sure <tt>rnd</tt> and <tt>x</tt> really are local variables!)
+  <li>Deep lexical scopes are bad. You can create a function without much of a lexical scope by using
+    <tt>new Function('...');</tt>.
+</ul>
+<br>
+
+<b>GC pauses?</b>
+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.
+<br><br><br>
+
+<table id="bench"></table>
+
+<!--[if lte IE 9]><script language="javascript" type="text/javascript" src="tools/excanvas.min.js"></script><![endif]-->
+<!--<script src="lib/numeric.js"></script>
+<script src="tools/sylvester.js"></script>
+<script src="tools/trunk/closure/goog/base.js"></script>
+<script src="tools/jquery-1.7.1.min.js"></script>
+<script src="tools/jquery.flot.min.js"></script>-->
+<script src="./tools/benchlib.js"></script>
+<script src="./src/iterative.js"></script>
+
+<script>
+"use strict";
+
+// Guess which browser needs this.
+if (!('map' in Array.prototype)) {
+    Array.prototype.map= function(mapper, that /*opt*/) {
+        var other= new Array(this.length);
+        for (var i= 0, n= this.length; i<n; i++)
+            if (i in this)
+                other[i]= mapper.call(that, this[i], i, this);
+        return other;
+    };
+}
+
+goog.require('goog.math.Matrix');
+var bench = numeric.bench;
+var geometricmeans = [0,0,0];
+var mkA = function(n) { var A = numeric.ccsScatter(numeric.cdelsq(numeric.cgrid(n))); A.format = 'ccs'; return A;};
+var mkV = function(n) { return numeric.random([(n-2)*(n-2)]); };
+var benchmarks = [
+  [
+  'Solve Ax=b, A is Sparse Laplacian', [4, 6, 9, 14, 21, 32, 48, 72],
+  function(n) { var A = mkA(n); var b = mkV(n); return bench(function() { numeric.ccsLUPSolve(numeric.ccsLUP(A), b); }); },
+  function(n) { var A = mkA(n); var b = mkV(n); return bench(function() { numeric.cg(A, b); }); },
+  function(n) { var A = mkA(n); var b = mkV(n); return bench(function() { numeric.bicgstab(A, b); }); },
+  ],
+  [
+  '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.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); }); },  
+  ],
+];
+
+var pwr = [4, 4];
+
+var k=0,b=0;
+var libs = ['ccsLUPSolve(A,b)','cg(A,b)', 'bicgstab(A,b)'];
+var datasets = [];
+var l;
+var colors = ["#000","#00f","#0f0","#f80","#f0f","#0ff"];
+for(l=0;l<libs.length;l++) { 
+  datasets[l] = { 
+    data: [], 
+    label: libs[l], 
+    color: colors[l], 
+    points: { show: true }, 
+    lines: { show: true } 
+    };
+}
+l=1;
+var l0,xticks = [];
+for(b=0;b<benchmarks.length;b++) {
+  l0 = l;
+  for(k=0;k<benchmarks[b][1].length;k++) {
+    l++;
+  }
+  xticks.push([(l0+l)*0.5,benchmarks[b][0]]);
+}
+k=0;
+b=0;
+var count = 1;
+var MSIE = (/MSIE (\d+\.\d+);/.test(navigator.userAgent));
+numeric.precision = 6;
+var c0 = [[],[],[]];
+var counts = [0,0,0];
+function invbench(b,k,lib,rep) {
+  var ks,sz = benchmarks[b][1][k];
+  var i,j,foo;
+  ks = sz.toString();
+  if(rep>0 && c0[lib][rep-1] < 10) { c0[lib][rep] = c0[lib][rep-1]; }
+  else { foo = benchmarks[b][lib+2]; c0[lib][rep] = foo(sz); }
+  rep++;
+  if(rep === 1) { rep = 0; lib++; }
+  if(lib+2 === benchmarks[b].length) {
+    k++;
+    lib=0;
+    var cps = [];
+    var mi = 1e6/Math.pow(benchmarks[b][1][k-1],pwr[b]);
+    for(i=0;i<c0.length;i++) {
+        cps[i] = 0;
+        for(j=0;j<c0[i].length;j++) cps[i] += c0[i][j];
+      cps[i] /= (c0[i].length*mi);
+    }
+    for(i=0;i<cps.length;i++) {
+      if(MSIE || cps[i]) datasets[i].data.push([count,cps[i]]);
+      if(cps[i]) {
+          counts[i]++;
+          geometricmeans[i] += Math.log(cps[i]);
+      }
+    }
+    var foo = '<td>n='+ks+'</td>';
+    var color = '', uncolor='';
+    for(i=0;i<cps.length;i++) {
+      foo +=  '<td>'+cps[i].toPrecision(8)+'</td>';
+    }
+    if(!MSIE) {
+      var t = document.getElementById('bench');
+      var r = t.insertRow(-1);
+      if(k === 1) { 
+        r.innerHTML = ('<td width=200px><b>'+benchmarks[b][0]+'</b></td>'
+                      +'<td width=125px><b>Numeric</b></td>'
+                        +'<td width=125px><b>Google Closure</b></td>'
+                        +'<td width=125px><b>Sylvester</b></td>');
+        r = t.insertRow(-1); 
+      }
+      r.innerHTML = foo;
+    }
+    $.plot($("#placeholder"), datasets, 
+        { 
+          legend: {container: '#legend'}, 
+          xaxis:  {ticks:xticks, tickLength:0, min:1, max: l-1},
+          yaxis:  {ticks:20}
+        });
+      c0 = [[],[],[]];
+    count++;
+  }
+  if(k === benchmarks[b][1].length) {
+    for(i=0;i<cps.length;i++) datasets[i].data.push(null);
+    k=0; b++;
+  }
+  if(b < benchmarks.length) { 
+    setTimeout('invbench('+b.toString()
+                      +','+k.toString()
+                      +','+lib.toString()
+                      +','+rep.toString()
+                      +')',MSIE?10:0); 
+    } else {
+        geometricmeans = numeric.exp(numeric.div(geometricmeans,counts));
+        $('#meanscore')[0].innerHTML += numeric.prettyPrint(geometricmeans)+'MFLOPS';
+    }
+}
+window.onload = function() { invbench(0,0,0,0); }
+</script>
+
+<br><br><br>
+</body>
diff --git a/src/iterative.js b/src/iterative.js
new file mode 100644
index 0000000..5b5f955
--- /dev/null
+++ b/src/iterative.js
@@ -0,0 +1,380 @@
+// 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<j1;j++) {
+			ri = Aj[j];
+			val = Av[j];
+			ret[ri] += x[k] * val;
+		}
+	}
+	return ret;
+}
+
+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, t, alpha, beta;
+		
+		var bnorm = dot(b, b);
+		
+        while( !converged && iters < maxIters ) {
+            var rDotr = dot(rhat, r);
+			
+            if( Math.abs(rDotr / bnorm) <= residue ) {
+                converged = true;
+                break;
+            }
+			
+            v = mv(A, p);
+            alpha = rDotr / dot(rhat, v);
+			s = axpy(-alpha, v, r);
+			
+			t = mv(A, s);	
+			w = dot(t, s) / dot(t, t);
+
+            x = axpy(w, s, axpy(alpha, p, x));
+            r = axpy(-w, t, s);
+
+            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;
+	
+	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);
+		
+		var bnorm = dot(b, b);
+		var rho = bnorm;
+		var flag = 1;
+		var Ap, pAp, alpha, beta, rho_new;
+        while( !converged && iters < maxIters ) {		
+            if( Math.abs(rho / bnorm) <= residue ) {
+                converged = true;
+                break;
+            }
+			
+            Ap = mv(A, p);
+            pAp = dot(p, Ap);
+            alpha = rho / pAp;
+
+            x = axpy(alpha, p, x);
+            r = axpy(-alpha, Ap, r);
+
+            rho_new = dot(r, r);
+            beta = rho_new / rho;
+            p = axpy(beta, p, r);
+			rho = rho_new;
+            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 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';
+	}
+	}	
+};
+
+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<n;i++) {
+				Ai = A[i];
+				rowsum = 0;
+				for(j=0;j<i;j++) {
+					rowsum += Ai[j] * x[j];
+				}
+				
+				for(j=i+1;j<m;j++) {
+					rowsum += Ai[j] * x[j];					
+				}
+				
+				rowdiff = b[i] - rowsum;
+				res[i] += (rowdiff - x[i] * Ai[i]);
+				x[i] = (1-relax) * x[i] + relax / Ai[i] * rowdiff;
+			}
+			
+			iters++;
+			// check convergence
+			r2 = dot(res, res);
+			converged = ((r2 / bnorm) < residue);
+		}
+		
+		if( converged )
+			console.log('converged in ' + iters + ' iterations');
+		else
+			console.log('not converged in ' + iters + ' iterations');
+			
+		return x;
+	}
+	
+	var sor_ccs = function(A, b) {
+		// initialization
+        var iters = 0;
+        var converged = false;
+        var rows = b.length;
+		
+		var Ai = A[0], Aj = A[1], Av = A[2];
+		var sA = numeric.ccsDim(A);
+		var m = sA[0], n = sA[1];
+		
+		if( m !== rows )
+			throw 'Matrix dimension does not match input vector.';
+		
+		var rep = numeric.rep;
+		var dot = numeric.dot;
+		var rowsum = rep([rows], 0);
+		var r = rep([rows], 0);
+		var i, j, k, j0, j1;
+		var bnorm = dot(b, b);
+
+		// initialize x
+        var x = rep([rows], 0);
+		
+		// find out diagonal
+		var Aii = rep([rows], 0);
+		for(k=0;k!==n;k++) {
+			j0 = Ai[k];
+			j1 = Ai[k+1];
+			for(j=j0;j<j1;j++) {
+				ri = Aj[j];
+				if( ri == k ) {
+					Aii[k] = Av[j];
+					break;
+				}
+			}
+		}
+			
+		while( !converged && iters < maxIters ) {
+			// collect row sum
+			for(k=0;k!==n;k++) {
+				j0 = Ai[k];
+				j1 = Ai[k+1];
+				for(j=j0;j<j1;j++) {
+					ri = Aj[j];
+					val = Av[j];
+					rowsum[ri] += x[k] * val;
+				}
+			}
+			
+			// subtract diagonal elements from row sum and update x
+			for(i=0;i<m;i++) {
+				r[i] = b[i] - rowsum[i];
+				x[i] = (1-relax) * x[i] + relax / Aii[i] * (r[i] + x[i] * Aii[i]);
+				rowsum[i] = 0;
+			}
+		
+			iters++;			
+			r2 = dot(r, r);
+			converged = (r2/bnorm <= residue);
+		}
+		
+		if( converged )
+			console.log('converged in ' + iters + ' iterations');
+		else
+			console.log('not converged in ' + iters + ' iterations');
+			
+		return x;
+	}
+	
+	var sor_crs = function(A, b) {
+		// initialization
+        var iters = 0;
+        var converged = false;
+        var rows = b.length;
+		
+		var Ai = A[0], Aj = A[1], Av = A[2];
+		var nelems = Ai.length;
+		var n = numeric.sup(Aj) + 1;
+		var m = numeric.sup(Ai) + 1;
+		if( m !== rows )
+			throw 'Matrix dimension does not match input vector.';
+		
+		var rep = numeric.rep;
+		var dot = numeric.dot;
+		var rowsum = rep([rows], 0);
+		var r = rep([rows], 0);
+		var i, j, k;
+		var bnorm = dot(b, b);
+
+		// initialize x
+        var x = rep([rows], 0);
+		
+		// find out diagonal
+		var Aii = rep([rows], 0);
+		for(k=0;k!==nelems;k++) {
+			i = Ai[k];
+			j = Aj[k];
+			if( i == j ) {
+				Aii[i] = Av[k];
+			}
+		}
+			
+		while( !converged && iters < maxIters ) {
+			// collect row sum
+			for(k=0;k!==nelems;k++) {
+				i = Ai[k];
+				j = Aj[k];
+				rowsum[i] += Av[k] * x[j];
+			}
+			
+			// subtract diagonal elements from row sum and update x
+			for(i=0;i<m;i++) {
+				r[i] = b[i] - rowsum[i];
+				x[i] = (1-relax) * x[i] + relax / Aii[i] * (r[i] + x[i] * Aii[i]);
+				rowsum[i] = 0;
+			}
+		
+			iters++;			
+			r2 = dot(r, r);
+			converged = (r2/bnorm <= residue);
+		}
+		
+		if( converged )
+			console.log('converged in ' + iters + ' iterations');
+		else
+			console.log('not converged in ' + iters + ' iterations');
+			
+		return x;
+	}
+
+	switch(A.format) {
+	case 'full': {	
+		return sor_full(A, b);
+	}
+	case 'ccs': {
+		return sor_ccs(A, b);
+	}
+	case 'crs': {
+		return sor_crs(A, b);
+	}
+	default: {
+		throw 'Not supported matrix format';
+	}
+	}	
+
+}
\ No newline at end of file
diff --git a/src/numeric.js b/src/numeric.js
index 537b68f..b613221 100644
--- a/src/numeric.js
+++ b/src/numeric.js
@@ -1879,6 +1879,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<j1;j++) {
+			ri = Aj[j];
+			val = Av[j];
+			ret[ri] += x[k] * val;
+		}
+	}
+	return ret;
+}
+
 numeric.ccsDot = function ccsDot(A,B) {
     var Ai = A[0], Aj = A[1], Av = A[2];
     var Bi = B[0], Bj = B[1], Bv = B[2];
