diff --git a/lib/type/matrix/CcsMatrix.js b/lib/type/matrix/CcsMatrix.js
index a0436bd8d..e95482a8e 100644
--- a/lib/type/matrix/CcsMatrix.js
+++ b/lib/type/matrix/CcsMatrix.js
@@ -1291,46 +1291,16 @@ function factory (type, config, load, typed) {
       // skip first column in upper triangular matrix
       if (j > 0) {
         // 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) {
+        spa.forEach(0, j, function (k, vkj) {
+          // loop rows in column k (L)
+          _forEachRow(k, lvalues, lindex, lptr, function (i, vik) {
+            // check row is below k
+            if (i > k) {
               // update spa value
-              spa.accumulate(i, unaryMinus(multiply(vix, vxj)));
+              spa.accumulate(i, unaryMinus(multiply(vik, vkj)));
             }
           });
-        });*/
-        
-        
-        // 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;
diff --git a/lib/type/matrix/Spa.js b/lib/type/matrix/Spa.js
index 9c418edc8..6f2eb4368 100644
--- a/lib/type/matrix/Spa.js
+++ b/lib/type/matrix/Spa.js
@@ -3,7 +3,14 @@
 function factory (type, config, load, typed) {
   
   var add = load(require('../../function/arithmetic/add'));
+  var equal = load(require('../../function/relational/equal'));
   
+  /**
+   * An ordered Sparse Accumulator is a representation for a sparse vector that includes a dense array 
+   * of the vector elements and an ordered list of non-zero elements.
+   *
+   * @param {Integer} length          The initial dense vector length.
+   */
   function Spa(length) {
     if (!(this instanceof Spa))
       throw new SyntaxError('Constructor must be called with the new operator');
@@ -17,65 +24,98 @@ function factory (type, config, load, typed) {
     // check we have a value @ i
     if (!this._values[i]) {
       // insert in heap
-      this._heap.insert(i);
+      var node = this._heap.insert(i, v);
+      // set the value @ i
+      this._values[i] = node;
+    }
+    else {
+      // update the value @ i
+      this._values[i].value = v;
     }
-    // set the value @ i
-    this._values[i] = v;
   };
   
   Spa.prototype.get = function (i) {
-    return this._values[i] || 0;
+    var node = this._values[i];
+    if (node)
+      return node.value;
+    return 0;
   };
   
   Spa.prototype.accumulate = function (i, v) {
-    // current value
-    var value = this._values[i];
-    // check we have a value @ i
-    if (!value) {
+    // node @ i
+    var node = this._values[i];
+    if (!node) {
       // insert in heap
-      this._heap.insert(i);
+      node = this._heap.insert(i);
       // initialize value
-      value = 0;
+      this._values[i] = node;
+    }
+    else {
+      // accumulate value
+      node.value = add(node.value, v);
     }
-    // accumulate value
-    this._values[i] = add(value, v);
   };
   
   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];
+    // references
+    var heap = this._heap;
+    var values = this._values;
+    // nodes
+    var nodes = [];
+    // current node
+    var node = heap.extractMinimum();
+    // extract all nodes from heap (ordered)
+    while (node) {
+      // save node
+      nodes.push(node);
       // 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);
+      if (node.key >= from && node.key <= to) {
+        // check value is not zero
+        if (!equal(node.value, 0)) {
+          // invoke callback
+          callback(node.key, node.value, this);
+        }
       }
+      // extract next node
+      node = heap.extractMinimum();
+    }
+    // reinsert all nodes in heap
+    for (var i = 0; i < nodes.length; i++) {
+      // current node
+      var n = nodes[i];
+      // insert node in heap
+      node = heap.insert(n.key, n.value);
+      // update values
+      values[n.key] = node;
     }
   };
   
   Spa.prototype.swap = function (i, j) {
-    // values @ i and j
-    var vi = this._values[i];
-    var vj = this._values[j];
-    // swap values
-    this._values[i] = vj;
-    this._values[j] = vi;
-    // chek we need to insert indeces
-    if (!vi && vj) {
+    // node @ i and j
+    var nodei = this._values[i];
+    var nodej = this._values[j];
+    // check we need to insert indeces
+    if (!nodei && nodej) {
       // insert in heap
-      this._heap.insert(i);
+      nodei = this._heap.insert(i, nodej.value);
       // remove from heap
-      this._heap.remove(j);
+      this._heap.remove(nodej);
+      // set values
+      this._values[i] = nodei;
     }
-    else if (vi && !vj) {
+    else if (nodei && !nodej) {
       // insert in heap
-      this._heap.insert(j);
+      nodej = this._heap.insert(j, nodei.value);
       // remove from heap
-      this._heap.remove(i);
+      this._heap.remove(nodei);
+      // set values
+      this._values[j] = nodej;
+    }
+    else if (nodei && nodej) {
+      // swap values
+      var v = nodei.value;
+      nodei.value = nodej.value;
+      nodej.value = v;
     }
   };