Repository: josdejong/mathjs
Pull request #321 (2015-04-09 21:43:21)

Matrix LU decomposition with partial pivoting and LU based linear equations solver
This pull request adds Matrix LU decomposition with partial pivoting and a LU based linear equations solver.

Syntax (CLI): 

`m = [1, 2, -1; 2, 1, 1; 1, 2, 1];`
`x = lusolve(m, [4; -2; 2])`

Using explicit LUP decomposition:

`m = [1, 2, -1; 2, 1, 1; 1, 2, 1];`
`d = lup(m);`
`x = lusolve(d, [4; -2; 2])` --> Solving linear equations system 1
`y = lusolve(d, [1; 1; 8])`  --> Solving linear equations system 2 (same LUP decomposition)

Algorithms have been optimized for Dense and Sparse Matrices:

`m = matrix([1, 2, -1; 2, 1, 1; 1, 2, 1], "ccs");`
`x = lusolve(m, [4; -2; 2])`





Patch:
diff --git a/index.js b/index.js
index e2d4b7a18..787701a51 100644
--- a/index.js
+++ b/index.js
@@ -38,23 +38,32 @@ function create (config) {
   math.type.Help = require('./lib/type/Help');
   math.type.ResultSet = require('./lib/type/ResultSet');
   math.import(require('./lib/type/BigNumber'));
-
+  math.import(require('./lib/type/FibonacciHeap'));
+  
+  // sparse accumulator
+  math.import(require('./lib/type/matrix/Spa'));
+  
   // matrix storage formats
   math.import(require('./lib/type/matrix/CcsMatrix'));
   math.import(require('./lib/type/matrix/CrsMatrix'));
   math.import(require('./lib/type/matrix/DenseMatrix'));
-
+  
   // matrix storage format registry
   math.type.Matrix._storage.ccs = math.type.CcsMatrix;
   math.type.Matrix._storage.crs = math.type.CrsMatrix;
   math.type.Matrix._storage.dense = math.type.DenseMatrix;
   math.type.Matrix._storage['default'] = math.type.DenseMatrix;
 
-
   // serialization utilities
   // math.json.*
   math.import(require('./lib/json/reviver'));
 
+  // functions - algebra/decomposition
+  math.import(require('./lib/function/algebra/decomposition/lup'));
+  
+  // functions - algebra/solver
+  math.import(require('./lib/function/algebra/solver/lusolve'));
+  
   // functions - arithmetic
   math.import(require('./lib/function/arithmetic/abs'));
   math.import(require('./lib/function/arithmetic/add'));
diff --git a/lib/function/algebra/decomposition/lup.js b/lib/function/algebra/decomposition/lup.js
new file mode 100644
index 000000000..f12168b75
--- /dev/null
+++ b/lib/function/algebra/decomposition/lup.js
@@ -0,0 +1,402 @@
+'use strict';
+
+var util = require('../../../util/index');
+
+var object = util.object;
+
+function factory (type, config, load, typed) {
+
+  var matrix = load(require('../../construction/matrix'));
+  var abs = load(require('../../arithmetic/abs'));
+  var add = load(require('../../arithmetic/add'));
+  var divideScalar = load(require('../../arithmetic/divideScalar'));
+  var multiply = load(require('../../arithmetic/multiply'));
+  var subtract = load(require('../../arithmetic/subtract'));
+  var larger = load(require('../../relational/larger'));
+  var equal = load(require('../../relational/equal'));
+  var unaryMinus = load(require('../../arithmetic/unaryMinus'));
+  
+  var CcsMatrix = type.CcsMatrix;
+  var CrsMatrix = type.CrsMatrix;
+  var DenseMatrix = type.DenseMatrix;
+  var Spa = type.Spa;
+  
+  /**
+   * Calculate the Matrix LU decomposition with pivoting. Matrix A is decomposed in three matrices (L, U, P) where
+   * A * P = L * U
+   *
+   * @param {Matrix | Array} A two dimensional matrix or array for which to get the LUP decomposition.
+   *
+   * @return {Array<Matrix>} The lower triangular matrix, the upper triangular matrix and the permutation matrix.
+   */
+  var lup = typed('lup', {
+    'Matrix': function (m) {
+      // process matrix storage format
+      switch (m.storage()) {
+        case 'dense':
+          return _denseLUP(m);
+        case 'ccs':
+          return _ccsLUP(m);
+        case 'crs':
+          return _crsLUP(m);
+      }
+    },
+    'Array': function (a) {
+      // create dense matrix from array
+      var m = matrix(a);
+      // lup, use matrix implementation
+      var r = lup(m);
+      // result
+      return {
+        L: r.L.valueOf(),
+        U: r.U.valueOf(),
+        P: r.P.valueOf()
+      };
+    }
+  });
+
+  var _denseLUP = function (m) {
+    // rows & columns
+    var rows = m._size[0];
+    var columns = m._size[1];
+    // minimum rows and columns
+    var n = Math.min(rows, columns);
+    // matrix array, clone original data
+    var data = object.clone(m._data);
+    // l matrix arrays
+    var ldata = [];
+    var lsize = [rows, n];
+    // u matrix arrays
+    var udata = [];
+    var usize = [n, columns];
+    // p (eye [n, n])
+    var p = DenseMatrix.diagonal([rows, rows], 1, 0);
+    // vars
+    var i, j, k;
+    // loop columns
+    for (j = 0; j < columns; j++) {
+      // skip first column in upper triangular matrix
+      if (j > 0) {
+        // loop rows
+        for (i = 0; i < rows; i++) {
+          // 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++) {
+            // s = l[i, k] - data[k, j]
+            s = add(s, multiply(data[i][k], data[k][j]));
+          }
+          data[i][j] = subtract(data[i][j], s);
+        }
+      }      
+      // row with larger value in cvector, row >= j
+      var pi = j;
+      var pabsv = 0;
+      var vjj = 0;
+      // loop rows
+      for (i = j; i < rows; i++) {
+        // data @ i, j
+        var v = data[i][j];
+        // absolute value
+        var absv = abs(v);
+        // value is greater than pivote value
+        if (larger(absv, pabsv)) {
+          // store row
+          pi = i;
+          // update max value
+          pabsv = absv;
+          // value @ [j, j]
+          vjj = v;
+        }
+      }
+      // swap rows (j <-> pi)
+      if (j !== pi) {
+        // swap values j <-> pi in P
+        DenseMatrix._swapRows(j, pi, p._data);
+        // swap j <-> pi in data
+        DenseMatrix._swapRows(j, pi, data);
+      }
+      // check column is in lower triangular matrix
+      if (j < rows) {
+        // loop rows (lower triangular matrix)
+        for (i = j + 1; i < rows; i++) {
+          // value @ i, j
+          var vij = data[i][j];
+          if (!equal(vij, 0)) {
+            // update data
+            data[i][j] = divideScalar(data[i][j], vjj);
+          }
+        }
+      }
+    }
+    // loop columns
+    for (j = 0; j < columns; j++) {
+      // loop rows
+      for (i = 0; i < rows; i++) {
+        // initialize row in arrays
+        if (j === 0) {
+          // check row exists in upper triangular matrix
+          if (i < columns) {
+            // U
+            udata[i] = [];
+          }
+          // L
+          ldata[i] = [];
+        }
+        // check we are in the upper triangular matrix
+        if (i < j) {
+          // check row exists in upper triangular matrix
+          if (i < columns) {
+            // U
+            udata[i][j] = data[i][j];
+          }
+          // check column exists in lower triangular matrix
+          if (j < rows) {
+            // L
+            ldata[i][j] = 0;
+          }
+          continue;
+        }
+        // diagonal value
+        if (i === j) {
+          // check row exists in upper triangular matrix
+          if (i < columns) {
+            // U
+            udata[i][j] = data[i][j];
+          }
+          // check column exists in lower triangular matrix
+          if (j < rows) {
+            // L
+            ldata[i][j] = 1;
+          }
+          continue;
+        }
+        // check row exists in upper triangular matrix
+        if (i < columns) {
+          // U
+          udata[i][j] = 0;
+        }
+        // check column exists in lower triangular matrix
+        if (j < rows) {
+          // L
+          ldata[i][j] = data[i][j];
+        }
+      }
+    }
+    // l matrix
+    var l =  new DenseMatrix({
+      data: ldata,
+      size: lsize
+    });
+    // u matrix
+    var u =  new DenseMatrix({
+      data: udata,
+      size: usize
+    });    
+    // return matrices
+    return { 
+      L: l, 
+      U: u, 
+      P: p, 
+      toString: function () {
+        return 'L: ' + this.L.toString() + '\nU: ' + this.U.toString() + '\nP: ' + this.P.toString();
+      }
+    };
+  };
+  
+  var _ccsLUP = function (m) {
+    // rows & columns
+    var rows = m._size[0];
+    var columns = m._size[1];
+    // minimum rows and columns
+    var n = Math.min(rows, columns);
+    // matrix arrays (will not be modified, thanks to permutation vector)
+    var values = m._values;
+    var index = m._index;
+    var ptr = m._ptr;
+    // l matrix arrays
+    var lvalues = [];
+    var lindex = [];
+    var lptr = [];
+    var lsize = [rows, n];
+    // u matrix arrays
+    var uvalues = [];
+    var uindex = [];
+    var uptr = [];
+    var usize = [n, columns];
+    // p matrix arrays
+    var pvalues = [];
+    var pindex = [];
+    var pptr = [];
+    var psize = [rows, rows];
+    // vars
+    var i, j, k;
+    // permutation vectors, (current index -> original index) and (original index -> current index)
+    var pv_co = new Array(rows);
+    var pv_oc = new Array(rows);
+    for (i = 0; i < rows; i++) {
+      pv_co[i] = i;
+      pv_oc[i] = i;
+    }
+    // swap indices in permutation vectors (condition x < y)!
+    var swapIndeces = function (x, y) {      
+      // find pv indeces getting data from x and y
+      var kx = pv_oc[x];
+      var ky = pv_oc[y];
+      // update permutation vector current -> original
+      pv_co[kx] = y;
+      pv_co[ky] = x;
+      // update permutation vector original -> current
+      pv_oc[x] = ky;
+      pv_oc[y] = kx;
+    };
+    // loop columns
+    for (j = 0; j < columns; j++) {
+      // sparse accumulator
+      var spa = new Spa(rows);            
+      // check lower triangular matrix has a value @ column j
+      if (j < rows) {
+        // update ptr
+        lptr.push(lvalues.length);
+        // first value in j column for lower triangular matrix
+        lvalues.push(1);
+        lindex.push(j);
+      }
+      // update ptr
+      uptr.push(uvalues.length);
+      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+      var k0 = ptr[j];
+      var k1 = ptr[j + 1];
+      // copy column j into sparse accumulator
+      for (k = k0; k < k1; k++) {
+        // row
+        i = index[k];
+        // copy column values into sparse accumulator (use permutation vector)
+        spa.set(pv_co[i], values[k]);
+      }
+      // skip first column in upper triangular matrix
+      if (j > 0) {
+        // loop rows in column j (above diagonal)
+        spa.forEach(0, j - 1, function (k, vkj) {
+          // loop rows in column k (L)
+          CcsMatrix._forEachRow(k, lvalues, lindex, lptr, function (i, vik) {
+            // check row is below k
+            if (i > k) {
+              // update spa value
+              spa.accumulate(i, unaryMinus(multiply(vik, vkj)));
+            }
+          });
+        });        
+      }
+      // row with larger value in spa, row >= j
+      var pi = j;
+      var vjj = spa.get(j);
+      var pabsv = abs(vjj);      
+      // loop values in spa (order by row, below diagonal)
+      spa.forEach(j + 1, rows - 1, function (x, v) {
+        // absolute value
+        var absv = abs(v);
+        // value is greater than pivote value
+        if (larger(absv, pabsv)) {
+          // store row
+          pi = x;
+          // update max value
+          pabsv = absv;
+          // value @ [j, j]
+          vjj = v;
+        }
+      });
+      // swap rows (j <-> pi)
+      if (j !== pi) {
+        // swap values j <-> pi in L
+        CcsMatrix._swapRows(j, pi, lsize[1], lvalues, lindex, lptr);
+        // swap values j <-> pi in U
+        CcsMatrix._swapRows(j, pi, usize[1], uvalues, uindex, uptr);
+        // swap values in spa
+        spa.swap(j, pi);
+        // update permutation vector (swap values @ j, pi)
+        swapIndeces(j, pi);
+      }
+      // loop values in spa (order by row)
+      spa.forEach(0, rows - 1, function (x, v) {
+        // check we are above diagonal
+        if (x <= j) {
+          // update upper triangular matrix
+          uvalues.push(v);
+          uindex.push(x);
+        }
+        else {
+          // update value
+          v = divideScalar(v, vjj);
+          // check value is non zero
+          if (!equal(v, 0)) {
+            // update lower triangular matrix
+            lvalues.push(v);
+            lindex.push(x);
+          }
+        }
+      });
+    }
+    // update ptrs
+    lptr.push(lvalues.length);
+    uptr.push(uvalues.length);
+    // generate P matrix
+    for (i = 0; i < rows; i++) {
+      // ptr
+      pptr.push(pvalues.length);
+      // value
+      pvalues.push(1);
+      // row
+      pindex.push(pv_co[i]);
+    }
+    // close ptr
+    pptr.push(pvalues.length);
+    // return matrices
+    return {
+      L: new CcsMatrix({
+        values: lvalues,
+        index: lindex,
+        ptr: lptr,
+        size: lsize
+      }), 
+      U: new CcsMatrix({
+        values: uvalues,
+        index: uindex,
+        ptr: uptr,
+        size: usize
+      }), 
+      P: new CcsMatrix({
+        values: pvalues,
+        index: pindex,
+        ptr: pptr,
+        size: psize
+      }), 
+      toString: function () {
+        return 'L: ' + this.L.toString() + '\nU: ' + this.U.toString() + '\nP: ' + this.P.toString();
+      }
+    };
+  };
+  
+  var _crsLUP = function (m) {
+    // none of the LUP decomposition algorithms are optimized for CRS. Converting from CRS to CCS performs very whell...
+    var ccs = new type.CcsMatrix(m);
+    // use lup from ccs
+    var r = lup(ccs);
+    // return LUP
+    return {
+      L: new CrsMatrix(r.L), 
+      U: new CrsMatrix(r.U), 
+      P: new CrsMatrix(r.P), 
+      toString: function () {
+        return 'L: ' + this.L.toString() + '\nU: ' + this.U.toString() + '\nP: ' + this.P.toString();
+      }
+    };
+  };
+  
+  return lup;
+}
+
+exports.name = 'lup';
+exports.factory = factory;
diff --git a/lib/function/algebra/solver/backwardSubstitution.js b/lib/function/algebra/solver/backwardSubstitution.js
new file mode 100644
index 000000000..476721428
--- /dev/null
+++ b/lib/function/algebra/solver/backwardSubstitution.js
@@ -0,0 +1,246 @@
+'use strict';
+
+function factory (type, config, load, typed) {
+
+  var matrix = load(require('../../construction/matrix'));
+  var divideScalar = load(require('../../arithmetic/divideScalar'));
+  var multiply = load(require('../../arithmetic/multiply'));
+  var subtract = load(require('../../arithmetic/subtract'));
+  var equal = load(require('../../relational/equal'));
+
+  var substitutionValidation = load(require('./substitutionValidation'));
+  
+  var CcsMatrix = type.CcsMatrix;
+  var CrsMatrix = type.CrsMatrix;
+  var DenseMatrix = type.DenseMatrix;
+
+  /**
+   * Solves the linear equation system by backward substitution. Matrix must be an upper triangular matrix.
+   *
+   * M * x = b
+   *
+   * @param {Matrix, Array}         A N x N matrix or array
+   * @param {Matrix, Array}         A column vector with the b values
+   *
+   * @return {Matrix}               A column vector with the linear system solution
+   */
+  var backwardSubstitution = typed('backwardSubstitution', {
+    'Matrix, Array | Matrix': function (m, b) {
+      // process matrix storage format
+      switch (m.storage()) {
+        case 'dense':
+          return _denseBackwardSubstitution(m, b);
+        case 'ccs':
+          return _ccsBackwardSubstitution(m, b);
+        case 'crs':
+          return _crsBackwardSubstitution(m, b);
+      }
+    },
+    'Array, Array | Matrix': function (a, b) {
+      // create dense matrix from array
+      var m = matrix(a);
+      // use matrix implementation
+      var r = backwardSubstitution(m, b);
+      // result
+      return r.valueOf();
+    }
+  });
+
+  var _denseBackwardSubstitution = function (m, b) {
+    // validate matrix and vector, return copy of column vector b
+    b = substitutionValidation(m, b);
+    // rows & columns
+    var rows = m._size[0];
+    var columns = m._size[1];
+    // result
+    var x = new Array(rows);
+    // arrays
+    var data = m._data;
+    // backward solve m * x = b, loop columns (backwards)
+    for (var j = columns - 1; j >= 0 ; j--) {
+      // b[j]
+      var bj = b[j] || 0;
+      // x[j]
+      var xj;
+      // backward substitution (outer product) avoids inner looping when bj == 0
+      if (!equal(bj, 0)) {
+        // value @ [j, j]
+        var vjj = data[j][j];
+        // check vjj
+        if (equal(vjj, 0)) {
+          // system cannot be solved
+          throw new Error('Linear system cannot be solved since matrix is singular');
+        }
+        // calculate xj
+        xj = divideScalar(bj, vjj);        
+        // loop rows
+        for (var i = j - 1; i >= 0; i--) {
+          // update copy of b
+          b[i] = subtract(b[i] || 0, multiply(xj, data[i][j]));
+        }
+      }
+      else {
+        // zero value @ j
+        xj = 0;
+      }
+      // update x
+      x[j] = [xj];
+    }
+    // return column vector
+    return new DenseMatrix({
+      data: x,
+      size: [rows, 1]
+    });
+  };
+  
+  var _ccsBackwardSubstitution = function (m, b) {
+    // validate matrix and vector, return copy of column vector b
+    b = substitutionValidation(m, b);
+    // rows & columns
+    var rows = m._size[0];
+    var columns = m._size[1];
+    // matrix arrays
+    var values = m._values;
+    var index = m._index;
+    var ptr = m._ptr;
+    // result arrays
+    var xvalues = [];
+    var xindex = [];
+    var xptr = [];
+    // vars
+    var i, k;
+    // init ptr
+    xptr.push(0);
+    // backward solve m * x = b, loop columns (backwards)
+    for (var j = columns - 1; j >= 0 ; j--) {
+      // b[j]
+      var bj = b[j] || 0;
+      // backward substitution (outer product) avoids inner looping when bj == 0
+      if (!equal(bj, 0)) {
+        // value @ [j, j]
+        var vjj = 0;
+        // first & last indeces in column
+        var f = ptr[j];
+        var l = ptr[j + 1];
+        // values in column, find value @ [j, j], loop backwards
+        for (k = l - 1; k >= f; k--) {
+          // row
+          i = index[k];
+          // check row
+          if (i === j) {
+            // update vjj
+            vjj = values[k];
+          }
+          else if (i < j) {
+            // exit loop
+            break;
+          }
+        }
+        // at this point we must have a value @ [j, j]
+        if (equal(vjj, 0)) {
+          // system cannot be solved, there is no value @ [j, j]
+          throw new Error('Linear system cannot be solved since matrix is singular');
+        }
+        // calculate xj
+        var xj = divideScalar(bj, vjj);
+        // values in column, continue from last loop
+        for (; k >= f; k--) {
+          // row
+          i = index[k];
+          // update copy of b
+          b[i] = subtract(b[i] || 0, multiply(xj, values[k]));
+        }
+        // check for non zero
+        if (!equal(xj, 0)) {
+          // push to values (at the beginning since we are looping backwards)
+          xvalues.unshift(xj);
+          // row
+          xindex.unshift(j);
+        }
+      }
+    }
+    // update ptr
+    xptr.push(xvalues.length);
+    // return column vector
+    return new CcsMatrix({
+      values: xvalues,
+      index: xindex,
+      ptr: xptr,
+      size: [rows, 1]
+    });
+  };
+  
+  var _crsBackwardSubstitution = function (m, b) {
+    // validate matrix and vector, initialize x as copy of column vector b
+    var x = substitutionValidation(m, b);
+    // rows & columns
+    var rows = m._size[0];
+    // matrix arrays
+    var values = m._values;
+    var index = m._index;
+    var ptr = m._ptr;
+    // result arrays
+    var xvalues = [];
+    var xindex = [];
+    var xptr = [];
+    // backward solve m * x = b
+    for (var i = rows - 1; i >= 0; i--) {
+      // update ptr
+      xptr.push(xvalues.length);
+      // initialize x
+      var xi = x[i] || 0;
+      // value @ [i, i]
+      var vii = 0;
+      // values in row (backwards, upper triangular)
+      for (var k = ptr[i + 1] - 1, f = ptr[i]; k >= f; k--) {
+        // column
+        var j = index[k];
+        // check column
+        if (j > i) {
+          // update x
+          xi = subtract(xi, multiply(values[k], x[j]));
+        }
+        else if (j === i) {
+          // set vii value
+          vii = values[k];
+          // exit loop
+          break; 
+        }
+        else {
+          // exit loop, singular matrix
+          break;
+        }
+      }
+      // check we have a value @ [i, i]
+      if (equal(vii, 0)) {
+        // system cannot be solved, there is no value @ [i, i]
+        throw new Error('Linear system cannot be solved since matrix is singular');
+      }
+      // update xi
+      xi = divideScalar(xi, vii);
+      // check for zero
+      if (!equal(xi, 0)) {
+        // push to values (at the beginning since we are looping backwards)
+        xvalues.unshift(xi);
+        // column vector
+        xindex.push(0);
+      }
+      // update dense vector
+      x[i] = xi;
+    }
+    // update ptr
+    xptr.push(xvalues.length);
+    // return column vector
+    return new CrsMatrix({
+      values: xvalues,
+      index: xindex,
+      ptr: xptr,
+      size: [rows, 1]
+    });
+  };
+  
+  return backwardSubstitution;
+}
+
+exports.name = 'backwardSubstitution';
+exports.factory = factory;
diff --git a/lib/function/algebra/solver/forwardSubstitution.js b/lib/function/algebra/solver/forwardSubstitution.js
new file mode 100644
index 000000000..ade98f3af
--- /dev/null
+++ b/lib/function/algebra/solver/forwardSubstitution.js
@@ -0,0 +1,245 @@
+'use strict';
+
+function factory (type, config, load, typed) {
+
+  var matrix = load(require('../../construction/matrix'));
+  var divideScalar = load(require('../../arithmetic/divideScalar'));
+  var multiply = load(require('../../arithmetic/multiply'));
+  var subtract = load(require('../../arithmetic/subtract'));
+  var equal = load(require('../../relational/equal'));
+
+  var substitutionValidation = load(require('./substitutionValidation'));
+
+  var CcsMatrix = type.CcsMatrix;
+  var CrsMatrix = type.CrsMatrix;
+  var DenseMatrix = type.DenseMatrix;
+
+  /** 
+   * Solves the linear equation system by forwards substitution. Matrix must be a lower triangular matrix.
+   *
+   * M * x = b
+   *
+   * @param {Matrix, Array}         A N x N matrix or array
+   * @param {Matrix, Array}         A column vector with the b values
+   *
+   * @return {Matrix}               A column vector with the linear system solution
+   */
+  var forwardSubstitution = typed('forwardSubstitution', {
+    'Matrix, Array | Matrix': function (m, b) {
+      // process matrix storage format
+      switch (m.storage()) {
+        case 'dense':
+          return _denseForwardSubstitution(m, b);
+        case 'ccs':
+          return _ccsForwardSubstitution(m, b);
+        case 'crs':
+          return _crsForwardSubstitution(m, b);
+      }
+    },
+    'Array, Array | Matrix': function (a, b) {
+      // create dense matrix from array
+      var m = matrix(a);
+      // use matrix implementation
+      var r = forwardSubstitution(m, b);
+      // result
+      return r.valueOf();
+    }
+  });
+
+  var _denseForwardSubstitution = function (m, b) {
+    // validate matrix and vector, return copy of column vector b
+    b = substitutionValidation(m, b);
+    // rows & columns
+    var rows = m._size[0];
+    var columns = m._size[1];
+    // result
+    var x = new Array(rows);
+    // data
+    var data = m._data;
+    // forward solve m * x = b, loop columns
+    for (var j = 0; j < columns; j++) {
+      // b[j]
+      var bj = b[j] || 0;
+      // x[j]
+      var xj;
+      // forward substitution (outer product) avoids inner looping when bj == 0
+      if (!equal(bj, 0)) {
+        // value @ [j, j]
+        var vjj = data[j][j];
+        // check vjj
+        if (equal(vjj, 0)) {
+          // system cannot be solved
+          throw new Error('Linear system cannot be solved since matrix is singular');
+        }
+        // calculate xj
+        xj = divideScalar(bj, vjj);
+        // loop rows
+        for (var i = j + 1; i < rows; i++) {
+          // update copy of b
+          b[i] = subtract(b[i] || 0, multiply(xj, data[i][j]));
+        }
+      }
+      else {
+        // zero @ j
+        xj = 0;
+      }
+      // update x
+      x[j] = [xj];
+    }
+    // return vector
+    return new DenseMatrix({
+      data: x,
+      size: [rows, 1]
+    });
+  };
+
+  var _ccsForwardSubstitution = function (m, b) {
+    // validate matrix and vector, return copy of column vector b
+    b = substitutionValidation(m, b);
+    // rows & columns
+    var rows = m._size[0];
+    var columns = m._size[1];
+    // matrix arrays
+    var values = m._values;
+    var index = m._index;
+    var ptr = m._ptr;
+    // result arrays
+    var xvalues = [];
+    var xindex = [];
+    var xptr = [];
+    // vars
+    var i, k;
+    // init ptr
+    xptr.push(0);
+    // forward solve m * x = b, loop columns
+    for (var j = 0; j < columns; j++) {
+      // b[j]
+      var bj = b[j] || 0;
+      // forward substitution (outer product) avoids inner looping when bj == 0
+      if (!equal(bj, 0)) {
+        // value @ [j, j]
+        var vjj = 0;
+        // last index in column
+        var l = ptr[j + 1];
+        // values in column, find value @ [j, j]
+        for (k = ptr[j]; k < l; k++) {
+          // row
+          i = index[k];
+          // check row
+          if (i === j) {
+            // update vjj
+            vjj = values[k];
+          }
+          else if (i > j) {
+            // exit loop
+            break;
+          }
+        }
+        // at this point we must have a value @ [j, j]
+        if (equal(vjj, 0)) {
+          // system cannot be solved, there is no value @ [j, j]
+          throw new Error('Linear system cannot be solved since matrix is singular');
+        }
+        // calculate xj
+        var xj = divideScalar(bj, vjj);
+        // values in column, continue from last loop
+        for (; k < l; k++) {
+          // row
+          i = index[k];
+          // update copy of b
+          b[i] = subtract(b[i] || 0, multiply(xj, values[k]));
+        }
+        // check for non zero
+        if (!equal(xj, 0)) {
+          // value @ row i
+          xvalues.push(xj);
+          // row
+          xindex.push(j);
+        }
+      }
+    }
+    // update ptr
+    xptr.push(xvalues.length);
+    // return column vector
+    return new CcsMatrix({
+      values: xvalues,
+      index: xindex,
+      ptr: xptr,
+      size: [rows, 1]
+    });
+  };
+
+  var _crsForwardSubstitution = function (m, b) {
+    // validate matrix and vector, initialize x as copy of column vector b
+    var x = substitutionValidation(m, b);
+    // rows & columns
+    var rows = m._size[0];
+    // matrix arrays
+    var values = m._values;
+    var index = m._index;
+    var ptr = m._ptr;
+    // result arrays
+    var xvalues = [];
+    var xindex = [];
+    var xptr = [];
+    // forward solve m * x = b
+    for (var i = 0; i < rows; i++) {
+      // update ptr
+      xptr.push(xvalues.length);
+      // initialize x
+      var xi = x[i] || 0;
+      // value @ [i, i]
+      var vii = 0;
+      // values in row
+      for (var k = ptr[i], l = ptr[i + 1]; k < l; k++) {
+        // column
+        var j = index[k];
+        // check column
+        if (j < i) {
+          // update x
+          xi = subtract(xi, multiply(values[k], x[j]));
+        }
+        else if (j === i) {
+          // set vii value
+          vii = values[k];
+          // exit loop
+          break;          
+        }
+        else {
+          // exit loop, singular matrix
+          break;
+        }
+      }
+      // check we have a value @ [i, i]
+      if (equal(vii, 0)) {
+        // system cannot be solved, there is no value @ [i, i]
+        throw new Error('Linear system cannot be solved since matrix is singular');
+      }
+      // update xi
+      xi = divideScalar(xi, vii);
+      // check for zero
+      if (!equal(xi, 0)) {
+        // push to values
+        xvalues.push(xi);
+        // column vector
+        xindex.push(0);
+      }
+      // update dense vector
+      x[i] = xi;
+    }
+    // update ptr
+    xptr.push(xvalues.length);
+    // return column vector
+    return new CrsMatrix({
+      values: xvalues,
+      index: xindex,
+      ptr: xptr,
+      size: [rows, 1]
+    });
+  };
+
+  return forwardSubstitution;
+}
+
+exports.name = 'forwardSubstitution';
+exports.factory = factory;
diff --git a/lib/function/algebra/solver/lusolve.js b/lib/function/algebra/solver/lusolve.js
new file mode 100644
index 000000000..f68200012
--- /dev/null
+++ b/lib/function/algebra/solver/lusolve.js
@@ -0,0 +1,92 @@
+'use strict';
+
+var isArray = Array.isArray;
+
+function factory (type, config, load, typed) {
+  
+  var matrix = load(require('../../construction/matrix'));
+  var lup = load(require('../decomposition/lup'));
+
+  var backwardSubstitution = load(require('./backwardSubstitution'));
+  var forwardSubstitution = load(require('./forwardSubstitution'));
+
+  /**
+   * Solves the linear system A * x = b where A is an [n x n] matrix and b is a [n] column vector.
+   *
+   * Syntax:
+   *
+   *    var x = math.lusolve(matrix, b)  // returns column vector with the solution to the linear system A * x = b
+   *
+   *    var x = math.lusolve(lup, b)     // returns column vector with the solution to the linear system A * x = b, lup = math.lup(A)
+   *
+   * Examples:
+   *
+   *    var m = [[1, 0, 0, 0], [0, 2, 0, 0], [0, 0, 3, 0], [0, 0, 0, 4]];
+   *
+   *    var x = math.lusolve(m, [-1, -1, -1, -1]);        // x = [-1, -0.5, -1/3, -0.25]
+   *
+   *    var f = math.lup(m);
+   *    var x1 = math.lusolve(f, [-1, -1, -1, -1]);       // x1 = [-1, -0.5, -1/3, -0.25]
+   *    var x2 = math.lusolve(f, [1, 2, 1, -1]);          // x2 = [1, 1, 1/3, -0.25]
+   *
+   * See also:
+   *
+   *    lup
+   *
+   * @param {Matrix | Array | Object} a      Invertible Matrix or the Matrix LUP decomposition
+   * @param {Matrix | Array} b               Column Vector
+   *
+   * @return {Matrix | Array}       Column vector with the solution to the linear system A * x = b
+   */
+  var lusolve = typed('lusolve', {
+    'Array, Array | Matrix': function (a, b) {
+      // convert a to matrix
+      a = matrix(a);
+      // matrix lup decomposition
+      var d = lup(a);
+      // solve
+      var x = _lusolve(d.L, d.U, d.P, b);
+      // convert result to array
+      return x.valueOf();
+    },
+    'Matrix, Array | Matrix': function (a, b) {
+      // matrix lup decomposition
+      var d = lup(a);
+      // solve
+      return _lusolve(d.L, d.U, d.P, b);
+    },
+    'Object, Array | Matrix': function (d, b) {
+      // solve
+      return _lusolve(d.L, d.U, d.P, b);
+    }
+  });
+  
+  var _toMatrix = function (a) {
+    // check it is a matrix
+    if (a instanceof type.Matrix)
+      return a;
+    // check array
+    if (isArray(a))
+      return matrix(a);
+    // throw
+    throw new TypeError('Invalid Matrix LUP decomposition');
+  };
+  
+  var _lusolve = function (l, u, p, b) {
+    // verify L, U, P
+    l = _toMatrix(l);
+    u = _toMatrix(u);
+    p = _toMatrix(p);
+    // modify column vector applying permutations
+    b = p.multiply(b);
+    // use forward substitution to resolve L * y = b
+    var y = forwardSubstitution(l, b);
+    // use backward substitution to resolve U * x = y
+    return backwardSubstitution(u, y);
+  };
+
+  return lusolve;
+}
+
+exports.name = 'lusolve';
+exports.factory = factory;
diff --git a/lib/function/algebra/solver/substitutionValidation.js b/lib/function/algebra/solver/substitutionValidation.js
new file mode 100644
index 000000000..95014ff5a
--- /dev/null
+++ b/lib/function/algebra/solver/substitutionValidation.js
@@ -0,0 +1,93 @@
+'use strict';
+
+var util = require('../../../util/index');
+
+var string = util.string;
+var array = util.array;
+var object = util.object;
+
+var isArray = Array.isArray;
+
+function factory (type, config, load, typed) {
+  
+  var Matrix = type.Matrix;
+  
+  /**
+   * Validates matrix and column vector b for backward/forward substitution algorithms.
+   *
+   * @param {Matrix} m            An N x N matrix
+   * @param {Array | Matrix} b    A column vector
+   *
+   * @return {Array}    Dense copy of column vector b
+   */
+  var substitutionValidation = function (m, b) {
+    // matrix size
+    var size = m.size();
+    // validate matrix dimensions
+    if (size.length !== 2)
+      throw new RangeError('Matrix must be two dimensional (size: ' + string.format(size) + ')');
+    // rows & columns
+    var rows = size[0];
+    var columns = size[1];
+    // validate rows & columns
+    if (rows !== columns) 
+      throw new RangeError('Matrix must be square (size: ' + string.format(size) + ')');
+    // check b is matrix
+    if (b instanceof Matrix) {
+      // matrix size
+      var msize = b.size();
+      // check matrix dimensions, vector
+      if (msize.length === 1) {
+        // check vector length
+        if (msize[0] !== rows)
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.');
+      }
+      else if (msize.length === 2) {
+        // array must be a column vector
+        if (msize[0] !== rows || msize[1] !== 1)
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.');
+      }
+      else {
+        // throw error
+        throw new RangeError('Dimension mismatch. Matrix columns must match vector length.');
+      }
+      // dense copy of column vector
+      var x = new Array(rows);
+      // copy values (skip zeros)
+      b.forEach(function (v, index) {
+        x[index[0]] = object.clone(v);
+      }, true);
+      // return vector
+      return x;
+    }
+    // check b is array
+    if (isArray(b)) {
+      // size
+      var asize = array.size(b);
+      // check matrix dimensions, vector
+      if (asize.length === 1) {
+        // check vector length
+        if (asize[0] !== rows)
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.');
+        // return copy of b
+        return object.clone(b);
+      }
+      if (asize.length === 2) {
+        // array must be a column vector
+        if (asize[0] !== rows || asize[1] !== 1)
+          throw new RangeError('Dimension mismatch. Matrix columns must match vector length.');
+        // loop rows in array, return copy of column vector
+        return b.map(function (v) {
+          return object.clone(v[0]);
+        });
+      }
+      // throw error
+      throw new RangeError('Dimension mismatch. Matrix columns must match vector length.');      
+    }
+  };
+  
+  return substitutionValidation;
+}
+
+exports.name = 'substitutionValidation';
+exports.factory = factory;
\ No newline at end of file
diff --git a/lib/function/construction/matrix.js b/lib/function/construction/matrix.js
index cb3e8bc3d..c654cf29d 100644
--- a/lib/function/construction/matrix.js
+++ b/lib/function/construction/matrix.js
@@ -59,10 +59,10 @@ function factory (type, config, load, typed) {
    */
   function _create(data, format) {
     // get storage format constructor
-    var constructor = type.Matrix.storage(format || 'default');
+    var M = type.Matrix.storage(format || 'default');
 
     // create instance
-    return new constructor(data);
+    return new M(data);
   }
 }
 
diff --git a/lib/type/FibonacciHeap.js b/lib/type/FibonacciHeap.js
new file mode 100644
index 000000000..20d6f3064
--- /dev/null
+++ b/lib/type/FibonacciHeap.js
@@ -0,0 +1,328 @@
+'use strict';
+
+function factory (type, config, load, typed) {
+  
+  var smaller = load(require('../function/relational/smaller'));
+  var larger = load(require('../function/relational/larger'));
+  
+  var oneOverLogPhi = 1.0 / Math.log((1.0 + Math.sqrt(5.0)) / 2.0);
+  
+  function FibonacciHeap() {
+    if (!(this instanceof FibonacciHeap))
+      throw new SyntaxError('Constructor must be called with the new operator');
+
+    // initialize fields
+    this._minimum = null;
+    this._size = 0;
+  }
+  
+  /**
+   * Inserts a new data element into the heap. No heap consolidation is
+   * performed at this time, the new node is simply inserted into the root
+   * list of this heap. Running time: O(1) actual.
+   */
+  FibonacciHeap.prototype.insert = function (key, value) {
+    // create node
+    var node = {
+      key: key,
+      value: value,
+      degree: 0
+    };
+    // check we have a node in the minimum
+    if (this._minimum) {
+      // minimum node
+      var minimum = this._minimum;
+      // update left & right of node
+      node.left = minimum;
+      node.right = minimum.right;
+      minimum.right = node;
+      node.right.left = node;
+      // update minimum node in heap if needed
+      if (smaller(key, minimum.key)) {
+        // node has a smaller key, use it as minimum
+        this._minimum = node;
+      }
+    }
+    else {
+      // set left & right
+      node.left = node;
+      node.right = node;
+      // this is the first node
+      this._minimum = node;
+    }
+    // increment number of nodes in heap
+    this._size++;
+    // return node
+    return node;
+  };
+
+  /**
+   * Returns the number of nodes in heap. Running time: O(1) actual.
+   */
+  FibonacciHeap.prototype.size = function () {
+    return this._size;
+  };
+
+  /**
+   * Removes all elements from this heap.
+   */
+  FibonacciHeap.prototype.clear = function () {
+    this._minimum = null;
+    this._size = 0;
+  };
+
+  /**
+   * Returns true if the heap is empty, otherwise false.
+   */
+  FibonacciHeap.prototype.isEmpty = function () {
+    return !!this._minimum;
+  };
+  
+  /**
+   * Extracts the node with minimum key from heap. Amortized running 
+   * time: O(log n).
+   */
+  FibonacciHeap.prototype.extractMinimum = function () {
+    // node to remove
+    var node = this._minimum;
+    // check we have a minimum
+    if (node === null)
+      return node;
+    // current minimum
+    var minimum = this._minimum;
+    // get number of children
+    var numberOfChildren = node.degree;
+    // pointer to the first child
+    var x = node.child;
+    // for each child of node do...
+    while (numberOfChildren > 0) {
+      // store node in right side
+      var tempRight = x.right;
+      // remove x from child list
+      x.left.right = x.right;
+      x.right.left = x.left;
+      // add x to root list of heap
+      x.left = minimum;
+      x.right = minimum.right;
+      minimum.right = x;
+      x.right.left = x;
+      // set Parent[x] to null
+      x.parent = null;
+      x = tempRight;
+      numberOfChildren--;
+    }
+    // remove node from root list of heap
+    node.left.right = node.right;
+    node.right.left = node.left;
+    // update minimum
+    if (node == node.right) {
+      // empty
+      minimum = null;
+    }
+    else {
+      // update minimum
+      minimum = node.right;
+      // we need to update the pointer to the root with minimum key
+      minimum = _findMinimumNode(minimum, this._size);
+    }
+    // decrement size of heap
+    this._size--;
+    // update minimum
+    this._minimum = minimum;
+    // return node
+    return node;
+  };
+  
+  /**
+   * Removes a node from the heap given the reference to the node. The trees
+   * in the heap will be consolidated, if necessary. This operation may fail
+   * to remove the correct element if there are nodes with key value -Infinity.
+   * Running time: O(log n) amortized.
+   */
+  FibonacciHeap.prototype.remove = function (node) {
+    // decrease key value
+    this._minimum = _decreaseKey(this._minimum, node, null);
+    // remove the smallest
+    this.extractMinimum();
+  };
+  
+  /**
+   * Decreases the key value for a heap node, given the new value to take on.
+   * The structure of the heap may be changed and will not be consolidated. 
+   * Running time: O(1) amortized.
+   */
+  var _decreaseKey = function (minimum, node, key) {
+    // set node key
+    node.key = key;
+    // get parent node
+    var parent = node.parent;
+    if (parent && smaller(node.key, parent.key)) {
+      // remove node from parent
+      _cut(minimum, node, parent);
+      // remove all nodes from parent to the root parent
+      _cascadingCut(minimum, parent);
+    }
+    // update minimum node if needed
+    if (smaller(node.key, minimum.key))
+      minimum = node;
+    // return minimum
+    return minimum;
+  };
+  
+  /**
+   * The reverse of the link operation: removes node from the child list of parent.
+   * This method assumes that min is non-null. Running time: O(1).
+   */
+  var _cut = function (minimum, node, parent) {
+    // remove node from parent children and decrement Degree[parent]
+    node.left.right = node.right;
+    node.right.left = node.left;
+    parent.degree--;
+    // reset y.child if necessary
+    if (parent.child == node)
+      parent.child = node.right;
+    // remove child if degree is 0
+    if (parent.degree === 0)
+      parent.child = null;
+    // add node to root list of heap
+    node.left = minimum;
+    node.right = minimum.right;
+    minimum.right = node;
+    node.right.left = node;
+    // set parent[node] to null
+    node.parent = null;
+    // set mark[node] to false
+    node.mark = false;
+  };
+  
+  /**
+   * Performs a cascading cut operation. This cuts node from its parent and then
+   * does the same for its parent, and so on up the tree.
+   * Running time: O(log n); O(1) excluding the recursion.
+   */
+  var _cascadingCut= function (minimum, node) {
+    // store parent node
+    var parent = node.parent;
+    // if there's a parent...
+    if (!parent)
+      return;
+    // if node is unmarked, set it marked
+    if (!node.mark) {
+      node.mark = true;
+    }
+    else {
+      // it's marked, cut it from parent
+      _cut(minimum, node, parent);
+      // cut its parent as well
+      _cascadingCut(parent);
+    }
+  };
+  
+  /**
+   * Make the first node a child of the second one. Running time: O(1) actual.
+   */
+  var _linkNodes = function (node, parent) {
+    // remove node from root list of heap
+    node.left.right = node.right;
+    node.right.left = node.left;
+    // make node a Child of parent
+    node.parent = parent;
+    if (!parent.child) {
+      parent.child = node;
+      node.right = node;
+      node.left = node;
+    }
+    else {
+      node.left = parent.child;
+      node.right = parent.child.right;
+      parent.child.right = node;
+      node.right.left = node;
+    }
+    // increase degree[parent]
+    parent.degree++;
+    // set mark[node] false
+    node.mark = false;
+  };
+  
+  var _findMinimumNode = function (minimum, size) {
+    // to find trees of the same degree efficiently we use an array of length O(log n) in which we keep a pointer to one root of each degree
+    var arraySize = Math.floor(Math.log(size) * oneOverLogPhi) + 1;
+    // create list with initial capacity
+    var array = new Array(arraySize);
+    // find the number of root nodes.
+    var numRoots = 0;
+    var x = minimum;
+    if (x) {
+      numRoots++;
+      x = x.right;
+      while (x !== minimum) {
+        numRoots++;
+        x = x.right;
+      }
+    }
+    // vars
+    var y;
+    // For each node in root list do...
+    while (numRoots > 0) {
+      // access this node's degree..
+      var d = x.degree;
+      // get next node
+      var next = x.right;
+      // check if there is a node already in array with the same degree
+      while (true) {
+        // get node with the same degree is any
+        y = array[d];
+        if (!y)
+          break;
+        // make one node with the same degree a child of the other, do this based on the key value.
+        if (larger(x.key, y.key)) {
+          var temp = y;
+          y = x;
+          x = temp;
+        }
+        // make y a child of x
+        _linkNodes(y, x);
+        // we have handled this degree, go to next one.
+        array[d] = null;
+        d++;
+      }
+      // save this node for later when we might encounter another of the same degree.
+      array[d] = x;
+      // move forward through list.
+      x = next;
+      numRoots--;
+    }
+    // Set min to null (effectively losing the root list) and reconstruct the root list from the array entries in array[].
+    minimum = null;
+    // loop nodes in array
+    for (var i = 0; i < arraySize; i++) {
+      // get current node
+      y = array[i];
+      if (!y)
+        continue;
+      // check if we have a linked list
+      if (minimum) {
+        // First remove node from root list.
+        y.left.right = y.right;
+        y.right.left = y.left;
+        // now add to root list, again.
+        y.left = minimum;
+        y.right = minimum.right;
+        minimum.right = y;
+        y.right.left = y;
+        // check if this is a new min.
+        if (smaller(y.key, minimum.key))
+          minimum = y;
+      }
+      else
+        minimum = y;
+    }
+    return minimum;
+  };
+  
+  return FibonacciHeap;
+}
+
+exports.name = 'FibonacciHeap';
+exports.path = 'type';
+exports.factory = factory;
diff --git a/lib/type/Matrix.js b/lib/type/Matrix.js
index a9cc8d801..436eef892 100644
--- a/lib/type/Matrix.js
+++ b/lib/type/Matrix.js
@@ -1,7 +1,13 @@
 'use strict';
 
-var string = require('../util/string');
+var util = require('../util/index');
+
+var string = util.string;
+var array = util.array;
+var object = util.object;
+
 var isString = string.isString;
+var isArray = Array.isArray;
 
 function factory (type, config, load, typed) {
   /**
@@ -256,7 +262,7 @@ function factory (type, config, load, typed) {
     // must be implemented by each of the Matrix implementations
     throw new Error('Cannot invoke multiply on a Matrix interface');
   };
-  
+   
   // exports
   return Matrix;
 }
diff --git a/lib/type/matrix/CcsMatrix.js b/lib/type/matrix/CcsMatrix.js
index 543fde767..fefb2baed 100644
--- a/lib/type/matrix/CcsMatrix.js
+++ b/lib/type/matrix/CcsMatrix.js
@@ -15,10 +15,11 @@ var isInteger = util.number.isInteger;
 var validateIndex = array.validateIndex;
 
 function factory (type, config, load, typed) {
+  
   var add = load(require('../../function/arithmetic/add'));
   var multiply = load(require('../../function/arithmetic/multiply'));
   var equal = load(require('../../function/relational/equal'));
-
+  
   var Index = type.Index;
   var BigNumber = type.BigNumber;
   var Matrix = type.Matrix;
@@ -28,18 +29,8 @@ function factory (type, config, load, typed) {
       throw new SyntaxError('Constructor must be called with the new operator');
 
     if (data instanceof Matrix) {
-      // check data is a CcsMatrix
-      if (data.type === 'CcsMatrix') {
-        // clone arrays
-        this._values = object.clone(data._values);
-        this._index = object.clone(data._index);
-        this._ptr = object.clone(data._ptr);
-        this._size = object.clone(data._size);
-      }
-      else {
-        // build from matrix data
-        _createFromArray(this, data.valueOf());
-      }
+      // create from matrix
+      _createFromMatrix(this, data);
     }
     else if (data && isArray(data.values) && isArray(data.index) && isArray(data.ptr) && isArray(data.size)) {
       // initialize fields
@@ -65,6 +56,71 @@ function factory (type, config, load, typed) {
     }
   }
   
+  var _createFromMatrix = function (matrix, source) {
+    // check matrix type
+    if (source.type === 'CcsMatrix') {
+      // clone arrays
+      matrix._values = object.clone(source._values);
+      matrix._index = object.clone(source._index);
+      matrix._ptr = object.clone(source._ptr);
+      matrix._size = object.clone(source._size);
+    }
+    else if (source.type === 'CrsMatrix') {
+      // crs arrays
+      var svalues = source._values;
+      var sindex = source._index;
+      var sptr = source._ptr;
+      var ssize = source._size;
+      // rows and columns
+      var rows = ssize[0];
+      var columns = ssize[1];
+      // non-zero elements
+      var nzl = svalues.length;
+      // ccs arrays
+      var values = new Array(nzl);
+      var index = new Array(nzl);
+      var ptr = new Array(columns + 1);
+      // vars
+      var i, j, k, n;
+      // initialize ptrs
+      for (j = 0; j <= columns; j++) 
+        ptr[j] = 0;
+      // column lengths
+      for (k = 0; k < nzl; k++)
+        ptr[sindex[k] + 1]++;
+      for (j = 0; j < columns; j++)
+        ptr[j + 1] += ptr[j];
+      // loop rows
+      for (i = 0; i < rows; i++) {
+        // loop values for row
+        for (k = sptr[i], n = sptr[i + 1]; k < n; k++) {
+          // column
+          j = sindex[k];
+          // index in ccs arrays
+          var ccsk = ptr[j]++;
+          // row index
+          index[ccsk] = i;
+          // value
+          values[ccsk] = object.clone(svalues[k]);
+        }
+      }
+      // shift back ptr
+      for (j = columns; j > 0; j--)
+        ptr[j] = ptr[j - 1];
+      // first value (must be zero)
+      ptr[0] = 0;
+      // matrix arrays
+      matrix._values = values;
+      matrix._index = index;
+      matrix._ptr = ptr;
+      matrix._size = [rows, columns];
+    }
+    else {
+      // build from matrix data
+      _createFromArray(matrix, source.valueOf());
+    }
+  };
+  
   var _createFromArray = function (matrix, data) {
     // initialize fields
     matrix._values = [];
@@ -1175,12 +1231,10 @@ function factory (type, config, load, typed) {
     var index = [];
     var ptr = [];
 
-    // create array with rows entries
-    var data = [];
-    for (var x = 0; x < rows; x++)
-      data[x] = 0;
     // loop value columns
     for (var z = 0; z < c; z++) {
+      // create sparse accumulator
+      var spa = new type.Spa(rows);
       // update ptr
       ptr.push(values.length);
       // do not traverse rows in matrix, it is not efficient in CCS
@@ -1193,20 +1247,15 @@ function factory (type, config, load, typed) {
           // row
           var i = matrix._index[k];
           // multiply & aggregate
-          data[i] = add(data[i], multiply(matrix._values[k], get(j, z)));
-        }          
-      }
-      // finished processing column z, compress results
-      for (var y = 0; y < rows; y++) {
-        // check value is different than zero
-        if (!equal(data[y], 0)) {
-          // push value
-          values.push(data[y]);
-          index.push(y);          
+          spa.accumulate(i, multiply(matrix._values[k], get(j, z)));
         }
-        // reset value
-        data[y] = 0;
       }
+      // finished processing column z
+      spa.forEach(0, rows -1, function (x, v) {
+        // push value
+        values.push(v);
+        index.push(x);
+      });     
     }
     // update ptr
     ptr.push(values.length);
@@ -1224,6 +1273,109 @@ function factory (type, config, load, typed) {
     });
   };
   
+  /**
+   * Swap rows i and j in Matrix.
+   *
+   * @param {Number} i       Matrix row index 1
+   * @param {Number} j       Matrix row index 2
+   *
+   * @return {Matrix}        The matrix reference
+   */
+  CcsMatrix.prototype.swapRows = function (i, j) {
+    // check index
+    if (!isNumber(i) || !isInteger(i) || !isNumber(j) || !isInteger(j)) {
+      throw new Error('Row index must be positive integers');
+    }
+    // check dimensions
+    if (this._size.length !== 2) {
+      throw new Error('Only two dimensional matrix is supported');
+    }
+    // validate index
+    validateIndex(i, this._size[0]);
+    validateIndex(j, this._size[0]);
+    
+    // swap rows
+    CcsMatrix._swapRows(i, j, this._size[1], this._values, this._index, this._ptr);
+    // return current instance
+    return this;
+  };
+  
+  /**
+   * Loop rows with data in column j.
+   *
+   * @param (Number} j            Column
+   * @param {Array} values        Matrix values
+   * @param {Array} index         Matrix row indeces
+   * @param {Array} ptr           Matrix column pointers
+   * @param {function} callback   Callback function invoked for every row in column j
+   */
+  CcsMatrix._forEachRow = function (j, values, index, ptr, callback) {
+    // indeces for column j
+    var k0 = ptr[j];
+    var k1 = ptr[j + 1];
+    // loop
+    for (var k = k0; k < k1; k++) {
+      // invoke callback
+      callback(index[k], values[k]);
+    }
+  };
+  
+  /**
+   * Swap rows x and y in CCS Matrix data structures.
+   *
+   * @param {Number} x         Matrix row index 1
+   * @param {Number} y         Matrix row index 2
+   * @param {Number} columns   Number of columns in matrix
+   * @param {Array} values     Matrix values
+   * @param {Array} index      Matrix row indeces
+   * @param {Array} ptr        Matrix column pointers
+   */
+  CcsMatrix._swapRows = function (x, y, columns, values, index, ptr) {
+    // loop columns
+    for (var j = 0; j < columns; j++) {
+      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+      var k0 = ptr[j];
+      var k1 = ptr[j + 1];
+      // find value index @ x
+      var kx = _getValueIndex(x, k0, k1, index);
+      // find value index @ x
+      var ky = _getValueIndex(y, k0, k1, index);
+      // check both rows exist in matrix
+      if (kx < k1 && ky < k1 && index[kx] === x && index[ky] === y) {
+        // swap values
+        var v = values[kx];
+        values[kx] = values[ky];
+        values[ky] = v;
+        // next column
+        continue;
+      }
+      // check x row exist & no y row
+      if (kx < k1 && index[kx] === x && (ky >= k1 || index[ky] !== y)) {
+        // value @ x
+        var vx = values[kx];
+        // insert value @ y
+        values.splice(ky, 0, vx);
+        index.splice(ky, 0, y);
+        // remove value @ x (adjust array index if needed)
+        values.splice(ky <= kx ? kx + 1 : kx, 1);
+        index.splice(ky <= kx ? kx + 1 : kx, 1);
+        // next column
+        continue;
+      }
+      // check y row exist & no x row
+      if (ky < k1 && index[ky] === y && (kx >= k1 || index[kx] !== x)) {
+        // value @ y
+        var vy = values[ky];
+        // insert value @ x
+        values.splice(kx, 0, vy);
+        index.splice(kx, 0, x);
+        // remove value @ y (adjust array index if needed)
+        values.splice(kx <= ky ? ky + 1 : ky, 1);
+        index.splice(kx <= ky ? ky + 1 : ky, 1);
+      }
+    }
+  };
+
   return CcsMatrix;
 }
 
diff --git a/lib/type/matrix/CrsMatrix.js b/lib/type/matrix/CrsMatrix.js
index 471959800..ac6e061a3 100644
--- a/lib/type/matrix/CrsMatrix.js
+++ b/lib/type/matrix/CrsMatrix.js
@@ -15,6 +15,7 @@ var isInteger = util.number.isInteger;
 var validateIndex = array.validateIndex;
 
 function factory (type, config, load, typed) {
+
   var add = load(require('../../function/arithmetic/add'));
   var multiply = load(require('../../function/arithmetic/multiply'));
   var equal = load(require('../../function/relational/equal'));
@@ -28,18 +29,8 @@ function factory (type, config, load, typed) {
       throw new SyntaxError('Constructor must be called with the new operator');
 
     if (data instanceof Matrix) {
-      // check data is a CrsMatrix
-      if (data.type === 'CrsMatrix') {
-        // clone arrays
-        this._values = object.clone(data._values);
-        this._index = object.clone(data._index);
-        this._ptr = object.clone(data._ptr);
-        this._size = object.clone(data._size);
-      }
-      else {
-        // build from matrix data
-        _createFromArray(this, data.valueOf());
-      }
+      // create from matrix
+      _createFromMatrix(this, data);
     }
     else if (data && isArray(data.values) && isArray(data.index) && isArray(data.ptr) && isArray(data.size)) {
       // initialize fields
@@ -65,6 +56,71 @@ function factory (type, config, load, typed) {
     }
   }
 
+  var _createFromMatrix = function (matrix, source) {
+    // check matrix type
+    if (source.type === 'CrsMatrix') {
+      // clone arrays
+      matrix._values = object.clone(source._values);
+      matrix._index = object.clone(source._index);
+      matrix._ptr = object.clone(source._ptr);
+      matrix._size = object.clone(source._size);
+    }
+    else if (source.type === 'CcsMatrix') {
+      // ccs arrays
+      var svalues = source._values;
+      var sindex = source._index;
+      var sptr = source._ptr;
+      var ssize = source._size;
+      // rows and columns
+      var rows = ssize[0];
+      var columns = ssize[1];
+      // non-zero elements
+      var nzl = svalues.length;
+      // ccs arrays
+      var values = new Array(nzl);
+      var index = new Array(nzl);
+      var ptr = new Array(rows + 1);
+      // vars
+      var i, j, k, n;
+      // initialize ptrs
+      for (i = 0; i <= rows; i++) 
+        ptr[i] = 0;
+      // row lengths
+      for (k = 0; k < nzl; k++)
+        ptr[sindex[k] + 1]++;
+      for (i = 0; i < rows; i++)
+        ptr[i + 1] += ptr[i];
+      // loop columns in ccs
+      for (j = 0; j < columns; j++) {
+        // loop values for column
+        for (k = sptr[j], n = sptr[j + 1]; k < n; k++) {
+          // row
+          i = sindex[k];
+          // index in crs arrays
+          var crsk = ptr[i]++;
+          // column index
+          index[crsk] = j;
+          // value
+          values[crsk] = object.clone(svalues[k]);
+        }
+      }
+      // shift back ptr
+      for (i = rows; i > 0; i--)
+        ptr[i] = ptr[i - 1];
+      // first value (must be zero)
+      ptr[0] = 0;
+      // matrix arrays
+      matrix._values = values;
+      matrix._index = index;
+      matrix._ptr = ptr;
+      matrix._size = [rows, columns];
+    }
+    else {
+      // build from matrix data
+      _createFromArray(matrix, source.valueOf());
+    }
+  };
+  
   var _createFromArray = function (matrix, data) {
     // initialize fields
     matrix._values = [];
@@ -1211,6 +1267,76 @@ function factory (type, config, load, typed) {
     });
   };
 
+  /**
+   * Swap rows i and j in Matrix.
+   *
+   * @param {Number} i       Matrix row index 1
+   * @param {Number} j       Matrix row index 2
+   *
+   * @return {Matrix}        The matrix reference
+   */
+  CrsMatrix.prototype.swapRows = function (i, j) {
+    // check index
+    if (!isNumber(i) || !isInteger(i) || !isNumber(j) || !isInteger(j)) {
+      throw new Error('Row index must be positive integers');
+    }
+    // check dimensions
+    if (this._size.length !== 2) {
+      throw new Error('Only two dimensional matrix is supported');
+    }
+    // validate index
+    validateIndex(i, this._size[0]);
+    validateIndex(j, this._size[0]);
+    
+    // swap rows
+    _swapRows(i, j, this._values, this._index, this._ptr);
+    return this;
+  };
+  
+  var _swapRows = function (x, y, values, index, ptr) {
+    // make sure x is before y
+    if (x > y) {
+      var t = x;
+      x = y;
+      y = t;
+    }
+    // find indeces for row x
+    var kx0 = ptr[x];
+    var kx1 = ptr[x + 1];
+    var lx = kx1 - kx0;
+    // find indeces for row y
+    var ky0 = ptr[y];
+    var ky1 = ptr[y + 1];
+    var ly = ky1 - ky0;
+    // vars
+    var i, yvalues, yindex, xvalues, xindex;
+    // splice row y
+    yvalues = values.splice(ky0, ly);
+    yindex = index.splice(ky0, ly);
+    // splice row x
+    xvalues = values.splice(kx0, lx);
+    xindex = index.splice(kx0, lx);
+    // insert row y
+    for (i = 0; i < ly; i++) {
+      values.splice(kx0 + i, 0, yvalues[i]);
+      index.splice(kx0 + i, 0, yindex[i]);
+    }
+    // recalculate x index
+    var k = ky0 - lx + ly;
+    // insert row x
+    for (i = 0; i < lx; i++) {
+      values.splice(k + i, 0, xvalues[i]);
+      index.splice(k + i, 0, xindex[i]);
+    }
+    // update ptr if needed
+    if (lx !== ly) {
+      // loop from x to y rows
+      for (i = x + 1; i <= y; i++) {
+        ptr[i] += -lx + ly;
+      }
+    }
+  };
+  
   return CrsMatrix;
 }
 
diff --git a/lib/type/matrix/DenseMatrix.js b/lib/type/matrix/DenseMatrix.js
index 9b2296ee2..79778ffc9 100644
--- a/lib/type/matrix/DenseMatrix.js
+++ b/lib/type/matrix/DenseMatrix.js
@@ -14,9 +14,13 @@ var isInteger = util.number.isInteger;
 var validateIndex = array.validateIndex;
 
 function factory (type, config, load, typed) {
+  
   var add = load(require('../../function/arithmetic/add'));
+  var divideScalar = load(require('../../function/arithmetic/divideScalar'));
   var multiply = load(require('../../function/arithmetic/multiply'));
-
+  var subtract = load(require('../../function/arithmetic/subtract'));
+  var equal = load(require('../../function/relational/equal'));
+  
   var Index = type.Index;
   var BigNumber = type.BigNumber;
   var Matrix = type.Matrix;
@@ -1058,6 +1062,46 @@ function factory (type, config, load, typed) {
     });
   };
   
+  /**
+   * Swap rows i and j in Matrix.
+   *
+   * @param {Number} i       Matrix row index 1
+   * @param {Number} j       Matrix row index 2
+   *
+   * @return {Matrix}        The matrix reference
+   */
+  DenseMatrix.prototype.swapRows = function (i, j) {
+    // check index
+    if (!isNumber(i) || !isInteger(i) || !isNumber(j) || !isInteger(j)) {
+      throw new Error('Row index must be positive integers');
+    }
+    // check dimensions
+    if (this._size.length !== 2) {
+      throw new Error('Only two dimensional matrix is supported');
+    }
+    // validate index
+    validateIndex(i, this._size[0]);
+    validateIndex(j, this._size[0]);
+
+    // swap rows
+    DenseMatrix._swapRows(i, j, this._data);
+    // return current instance
+    return this;
+  };
+
+  /**
+   * Swap rows i and j in Dense Matrix data structure.
+   *
+   * @param {Number} i       Matrix row index 1
+   * @param {Number} j       Matrix row index 2
+   */
+  DenseMatrix._swapRows = function (i, j, data) {
+    // swap values i <-> j
+    var vi = data[i];
+    data[i] = data[j];
+    data[j] = vi;
+  };
+   
   /**
    * Preprocess data, which can be an Array or DenseMatrix with nested Arrays and
    * Matrices. Replaces all nested Matrices with Arrays
diff --git a/lib/type/matrix/Spa.js b/lib/type/matrix/Spa.js
new file mode 100644
index 000000000..72ff0bf92
--- /dev/null
+++ b/lib/type/matrix/Spa.js
@@ -0,0 +1,137 @@
+'use strict';
+
+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');
+    
+    // allocate vector, TODO use typed arrays
+    this._values = new Array(length);
+    this._heap = new type.FibonacciHeap();
+  }
+  
+  /**
+   * Set the value for index i.
+   *
+   * @param {Number} i                       The index
+   * @param {Number | BigNumber | Complex}   The value at index i
+   */
+  Spa.prototype.set = function (i, v) {
+    // check we have a value @ i
+    if (!this._values[i]) {
+      // insert in heap
+      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;
+    }
+  };
+  
+  Spa.prototype.get = function (i) {
+    var node = this._values[i];
+    if (node)
+      return node.value;
+    return 0;
+  };
+  
+  Spa.prototype.accumulate = function (i, v) {
+    // node @ i
+    var node = this._values[i];
+    if (!node) {
+      // insert in heap
+      node = this._heap.insert(i, v);
+      // initialize value
+      this._values[i] = node;
+    }
+    else {
+      // accumulate value
+      node.value = add(node.value, v);
+    }
+  };
+  
+  Spa.prototype.forEach = function (from, to, callback) {
+    // references
+    var heap = this._heap;
+    var values = this._values;
+    // nodes
+    var nodes = [];
+    // node with minimum key, save it
+    var node = heap.extractMinimum();
+    if (node)
+      nodes.push(node);
+    // extract nodes from heap (ordered)
+    while (node && node.key <= to) {
+      // check it is in range
+      if (node.key >= from) {
+        // check value is not zero
+        if (!equal(node.value, 0)) {
+          // invoke callback
+          callback(node.key, node.value, this);
+        }
+      }
+      // extract next node, save it
+      node = heap.extractMinimum();
+      if (node)
+        nodes.push(node);
+    }
+    // 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[node.key] = node;
+    }
+  };
+  
+  Spa.prototype.swap = function (i, j) {
+    // 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
+      nodei = this._heap.insert(i, nodej.value);
+      // remove from heap
+      this._heap.remove(nodej);
+      // set values
+      this._values[i] = nodei;
+      this._values[j] = undefined;
+    }
+    else if (nodei && !nodej) {
+      // insert in heap
+      nodej = this._heap.insert(j, nodei.value);
+      // remove from heap
+      this._heap.remove(nodei);
+      // set values
+      this._values[j] = nodej;
+      this._values[i] = undefined;
+    }
+    else if (nodei && nodej) {
+      // swap values
+      var v = nodei.value;
+      nodei.value = nodej.value;
+      nodej.value = v;
+    }
+  };
+  
+  return Spa;
+}
+
+exports.name = 'Spa';
+exports.path = 'type';
+exports.factory = factory;
diff --git a/lib/util/typed.js b/lib/util/typed.js
index e069db654..d767f6585 100644
--- a/lib/util/typed.js
+++ b/lib/util/typed.js
@@ -19,7 +19,7 @@ exports.create = function create(math) {
   typed.types['Help']       = function (x) { return x instanceof math.type.Help; };
   typed.types['ResultSet']  = function (x) { return x instanceof math.type.ResultSet; };
   typed.types['BigNumber']  = function (x) { return x instanceof math.type.BigNumber; };
-
+  
   typed.conversions = [
     {
       from: 'number',


Test patch:
diff --git a/test/function/algebra/decomposition/lup.test.js b/test/function/algebra/decomposition/lup.test.js
new file mode 100644
index 000000000..9e286c8c9
--- /dev/null
+++ b/test/function/algebra/decomposition/lup.test.js
@@ -0,0 +1,640 @@
+// test lup
+var assert = require('assert'),
+    approx = require('../../../../tools/approx'),
+    math = require('../../../../index');
+
+describe('lup', function () {
+
+  it('should decompose matrix, n x n, no permutations, array', function () {
+
+    var a = [[2, 1], [1, 4]];
+
+    var r = math.lup(a);
+    // L
+    assert.deepEqual(r.L.valueOf(), [[1, 0], [0.5, 1]]);
+    // U
+    assert.deepEqual(r.U.valueOf(), [[2, 1], [0, 3.5]]);
+    // P
+    assert.deepEqual(r.P.valueOf(), [[1, 0], [0, 1]]);
+    // verify
+    approx.deepEqual(math.multiply(r.P, a), math.multiply(r.L, r.U));
+  });
+
+  it('should decompose matrix, n x n, no permutations, ccs format', function () {
+
+    var m = math.matrix([[2, 1], [1, 4]], 'ccs');
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(r.L.valueOf(), [[1, 0], [0.5, 1]]);
+    // U
+    assert.deepEqual(r.U.valueOf(), [[2, 1], [0, 3.5]]);
+    // P
+    assert.deepEqual(r.P.valueOf(), [[1, 0], [0, 1]]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, n x n, no permutations, crs format', function () {
+
+    var m = math.matrix([[2, 1], [1, 4]], 'crs');
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(r.L.valueOf(), [[1, 0], [0.5, 1]]);
+    // U
+    assert.deepEqual(r.U.valueOf(), [[2, 1], [0, 3.5]]);
+    // P
+    assert.deepEqual(r.P.valueOf(), [[1, 0], [0, 1]]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, n x n, no permutations, dense format', function () {
+
+    var m = math.matrix([[2, 1], [1, 4]], 'dense');
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(r.L.valueOf(), [[1, 0], [0.5, 1]]);
+    // U
+    assert.deepEqual(r.U.valueOf(), [[2, 1], [0, 3.5]]);
+    // P
+    assert.deepEqual(r.P.valueOf(), [[1, 0], [0, 1]]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+  
+  it('should decompose matrix, m x n, m < n, no permutations, dense format', function () {
+    var m = math.matrix(
+      [
+        [2, 1, 1],
+        [1, 4, 5]
+      ]
+    );
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(
+      r.L,
+      math.matrix(
+        [
+          [1, 0],
+          [0.5, 1]
+        ]
+      ));
+    // U
+    assert.deepEqual(
+      r.U,
+      math.matrix(
+        [
+          [2, 1, 1],
+          [0, 3.5, 4.5]
+        ]
+      ));
+    // P
+    assert.deepEqual(
+      r.P,
+      math.matrix(
+        [
+          [1, 0],
+          [0, 1]
+        ]
+      ));
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, m x n, m > n, no permutations, dense format', function () {
+    var m = math.matrix(
+      [
+        [8, 2],
+        [6, 4],
+        [4, 1]
+      ]
+    );
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(
+      r.L,
+      math.matrix(
+        [
+          [1, 0],
+          [0.75, 1],
+          [0.5, 0]
+        ]
+      ));
+    // U
+    assert.deepEqual(
+      r.U,
+      math.matrix(
+        [
+          [8, 2],
+          [0, 2.5]
+        ]
+      ));
+    // P
+    assert.deepEqual(
+      r.P,
+      math.matrix(
+        [
+          [1, 0, 0],
+          [0, 1, 0],
+          [0, 0, 1]
+        ]
+      ));
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, n x n, dense format', function () {
+    var m = math.matrix(
+      [
+        [16, -120, 240, -140],
+        [-120, 1200, -2700, 1680],
+        [240, -2700, 6480, -4200],
+        [-140, 1680, -4200, 2800]
+      ]
+    );
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0, 0, 0],  
+        [-0.5, 1, 0, 0],
+        [-0.5833333333333334, -0.7, 1, 0],
+        [0.06666666666666667, -0.4, -0.5714285714285776, 1]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [240, -2700, 6480, -4200],
+        [0, -150, 540, -420], 
+        [0, 0, -42, 56],
+        [0, 0, 0, 4]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [0, 0, 1, 0],
+        [0, 1, 0, 0],      
+        [0, 0, 0, 1],
+        [1, 0, 0, 0]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, 3 x 3, zero pivote value, dense format', function () {
+    var m = math.matrix(
+      [
+        [1, 2, 3], 
+        [2, 4, 6], 
+        [4, 8, 9]
+      ]);
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0, 0],  
+        [0.5, 1, 0],
+        [0.25, 0, 1.0]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [4, 8, 9],
+        [0, 0, 1.5], 
+        [0, 0, 0.75]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [0, 0, 1],
+        [0, 1, 0],      
+        [1, 0, 0]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, 3 x 2, complex numbers, dense format', function () {
+    var m = math.matrix(
+      [
+        [math.complex(0, 3), 10],
+        [math.complex(0, 1), 1],
+        [math.complex(0, 1), 1]
+      ]);
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0],
+        [math.complex(0.3333333, 0), 1],
+        [math.complex(0.3333333, 0), 1]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [math.complex(0, 3), 10],
+        [0, math.complex(-2.3333333333, 0)]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+  
+  it('should decompose matrix, m x n, m < n, no permutations, ccs format', function () {
+    var m = math.matrix(
+      [
+        [2, 1, 1],
+        [1, 4, 5]
+      ], 
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0],
+        [0.5, 1]
+      ]);
+    // U
+    assert.deepEqual(
+      r.U.valueOf(),
+      [
+        [2, 1, 1],
+        [0, 3.5, 4.5]
+      ]);
+    // P
+    assert.deepEqual(
+      r.P.valueOf(),
+      [
+        [1, 0],
+        [0, 1]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, m x n, m > n, no permutations, ccs format', function () {
+    var m = math.matrix(
+      [
+        [8, 2],
+        [6, 4],
+        [4, 1]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0],
+        [0.75, 1],
+        [0.5, 0]
+      ]);
+    // U
+    assert.deepEqual(
+      r.U.valueOf(),
+      [
+        [8, 2],
+        [0, 2.5]
+      ]);
+    // P
+    assert.deepEqual(
+      r.P.valueOf(),
+      [
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, n x n, ccs format', function () {
+    var m = math.matrix(
+      [
+        [16, -120, 240, -140],
+        [-120, 1200, -2700, 1680],
+        [240, -2700, 6480, -4200],
+        [-140, 1680, -4200, 2800]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0, 0, 0],  
+        [-0.5, 1, 0, 0],
+        [-0.5833333333333334, -0.7, 1, 0],
+        [0.06666666666666667, -0.4, -0.5714285714285776, 1]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [240, -2700, 6480, -4200],
+        [0, -150, 540, -420], 
+        [0, 0, -42, 56],
+        [0, 0, 0, 4]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [0, 0, 1, 0],
+        [0, 1, 0, 0],      
+        [0, 0, 0, 1],
+        [1, 0, 0, 0]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, 3 x 3, zero pivote value, ccs format', function () {
+    var m = math.matrix(
+      [
+        [1, 2, 3], 
+        [2, 4, 6], 
+        [4, 8, 9]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0, 0],  
+        [0.5, 1, 0],
+        [0.25, 0, 1.0]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [4, 8, 9],
+        [0, 0, 1.5], 
+        [0, 0, 0.75]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [0, 0, 1],
+        [0, 1, 0],      
+        [1, 0, 0]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, 3 x 2, complex numbers, ccs format', function () {
+    var m = math.matrix(
+      [
+        [math.complex(0, 3), 10],
+        [math.complex(0, 1), 1],
+        [math.complex(0, 1), 1]
+      ], 'ccs');
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0],
+        [math.complex(0.3333333, 0), 1],
+        [math.complex(0.3333333, 0), 1]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [math.complex(0, 3), 10],
+        [0, math.complex(-2.3333333333, 0)]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+  
+  it('should decompose matrix, m x n, m < n, no permutations, crs format', function () {
+    var m = math.matrix(
+      [
+        [2, 1, 1],
+        [1, 4, 5]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0],
+        [0.5, 1]
+      ]);
+    // U
+    assert.deepEqual(
+      r.U.valueOf(),
+      [
+        [2, 1, 1],
+        [0, 3.5, 4.5]
+      ]);
+    // P
+    assert.deepEqual(
+      r.P.valueOf(),
+      [
+        [1, 0],
+        [0, 1]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, m x n, m > n, no permutations, crs format', function () {
+    var m = math.matrix(
+      [
+        [8, 2],
+        [6, 4],
+        [4, 1]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    assert.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0],
+        [0.75, 1],
+        [0.5, 0]
+      ]);
+    // U
+    assert.deepEqual(
+      r.U.valueOf(),
+      [
+        [8, 2],
+        [0, 2.5]
+      ]);
+    // P
+    assert.deepEqual(
+      r.P.valueOf(),
+      [
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, n x n, crs format', function () {
+    var m = math.matrix(
+      [
+        [16, -120, 240, -140],
+        [-120, 1200, -2700, 1680],
+        [240, -2700, 6480, -4200],
+        [-140, 1680, -4200, 2800]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0, 0, 0],  
+        [-0.5, 1, 0, 0],
+        [-0.5833333333333334, -0.7, 1, 0],
+        [0.06666666666666667, -0.4, -0.5714285714285776, 1]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [240, -2700, 6480, -4200],
+        [0, -150, 540, -420], 
+        [0, 0, -42, 56],
+        [0, 0, 0, 4]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [0, 0, 1, 0],
+        [0, 1, 0, 0],      
+        [0, 0, 0, 1],
+        [1, 0, 0, 0]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, 3 x 3, zero pivote value, crs format', function () {
+    var m = math.matrix(
+      [
+        [1, 2, 3], 
+        [2, 4, 6], 
+        [4, 8, 9]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0, 0],  
+        [0.5, 1, 0],
+        [0.25, 0, 1.0]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [4, 8, 9],
+        [0, 0, 1.5], 
+        [0, 0, 0.75]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [0, 0, 1],
+        [0, 1, 0],      
+        [1, 0, 0]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+
+  it('should decompose matrix, 3 x 2, complex numbers, crs format', function () {
+    var m = math.matrix(
+      [
+        [math.complex(0, 3), 10],
+        [math.complex(0, 1), 1],
+        [math.complex(0, 1), 1]
+      ],
+      'ccs');
+
+    var r = math.lup(m);
+    // L
+    approx.deepEqual(
+      r.L.valueOf(),
+      [
+        [1, 0],
+        [math.complex(0.3333333, 0), 1],
+        [math.complex(0.3333333, 0), 1]
+      ]);
+    // U
+    approx.deepEqual(
+      r.U.valueOf(),
+      [
+        [math.complex(0, 3), 10],
+        [0, math.complex(-2.3333333333, 0)]
+      ]);
+    // P
+    approx.deepEqual(
+      r.P.valueOf(),
+      [
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1]
+      ]);
+    // verify
+    approx.deepEqual(r.P.multiply(m), r.L.multiply(r.U));
+  });
+});
diff --git a/test/function/algebra/solver/lusolve.test.js b/test/function/algebra/solver/lusolve.test.js
new file mode 100644
index 000000000..f791ab08b
--- /dev/null
+++ b/test/function/algebra/solver/lusolve.test.js
@@ -0,0 +1,339 @@
+// test lusolve
+var assert = require('assert'),
+    approx = require('../../../../tools/approx'),
+    math = require('../../../../index');
+
+describe('lusolve', function () {
+  
+  it('should solve linear system 4 x 4, arrays', function () {
+    var m = 
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ];
+    var b = [-1, -1, -1, -1];
+    
+    var x = math.lusolve(m, b);
+    
+    approx.deepEqual(x, [-1, -0.5, -1/3, -0.25]);
+  });
+  
+  it('should solve linear system 4 x 4, array and column array', function () {
+    var m = 
+        [
+          [1, 0, 0, 0],
+          [0, 2, 0, 0],
+          [0, 0, 3, 0],
+          [0, 0, 0, 4]
+        ];
+    var b = [
+      [-1],
+      [-1], 
+      [-1],
+      [-1]
+    ];
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, [[-1], [-0.5], [-1/3], [-0.25]]);
+  });
+  
+  it('should solve linear system 4 x 4, matrices', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ]);
+    var b = math.matrix([-1, -1, -1, -1]);
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1/3], [-0.25]]));
+  });
+  
+  it('should solve linear system 4 x 4, crs matrices', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'crs');
+    var b = math.matrix([[-1], [-1], [-1], [-1]], 'crs');
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'crs'));
+  });
+  
+  it('should solve linear system 4 x 4, ccs matrices', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'ccs');
+    var b = math.matrix([[-1], [-1], [-1], [-1]], 'ccs');
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'ccs'));
+  });
+  
+  it('should solve linear system 4 x 4, matrix and column matrix', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ]);
+    var b = math.matrix([
+      [-1],
+      [-1], 
+      [-1],
+      [-1]
+    ]);
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1/3], [-0.25]]));
+  });
+  
+  it('should solve linear system 4 x 4, crs matrix and column matrix', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'crs');
+    var b = math.matrix([
+      [-1],
+      [-1], 
+      [-1],
+      [-1]
+    ], 'crs');
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'crs'));
+  });
+  
+  it('should solve linear system 4 x 4, ccs matrix and column matrix', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'ccs');
+    var b = math.matrix([
+      [-1],
+      [-1], 
+      [-1],
+      [-1]
+    ], 'ccs');
+
+    var x = math.lusolve(m, b);
+
+    assert(x instanceof math.type.Matrix);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'ccs'));
+  });
+  
+  it('should solve linear system 4 x 4, LUP decomposition (array)', function () {
+    var m = 
+        [
+          [1, 0, 0, 0],
+          [0, 2, 0, 0],
+          [0, 0, 3, 0],
+          [0, 0, 0, 4]
+        ];    
+    var lup = math.lup(m);
+    
+    var x = math.lusolve(lup, [-1, -1, -1, -1]);
+    approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1/3], [-0.25]]));
+    
+    var y = math.lusolve(lup, [1, 2, 1, -1]);
+    approx.deepEqual(y, math.matrix([[1], [1], [1/3], [-0.25]]));
+  });
+  
+  it('should solve linear system 4 x 4, LUP decomposition (matrix)', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ]);    
+    var lup = math.lup(m);
+
+    var x = math.lusolve(lup, [-1, -1, -1, -1]);
+    approx.deepEqual(x, math.matrix([[-1], [-0.5], [-1/3], [-0.25]]));
+
+    var y = math.lusolve(lup, [1, 2, 1, -1]);
+    approx.deepEqual(y, math.matrix([[1], [1], [1/3], [-0.25]]));
+  });
+  
+  it('should solve linear system 4 x 4, LUP decomposition (crs matrix)', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'crs');    
+    var lup = math.lup(m);
+
+    var x = math.lusolve(lup, [-1, -1, -1, -1]);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'crs'));
+
+    var y = math.lusolve(lup, [1, 2, 1, -1]);
+    approx.deepEqual(y, math.matrix([1, 1, 1/3, -0.25], 'crs'));
+  });
+  
+  it('should solve linear system 4 x 4, LUP decomposition (ccs matrix)', function () {
+    var m = math.matrix(
+      [
+        [1, 0, 0, 0],
+        [0, 2, 0, 0],
+        [0, 0, 3, 0],
+        [0, 0, 0, 4]
+      ], 'ccs');    
+    var lup = math.lup(m);
+
+    var x = math.lusolve(lup, [-1, -1, -1, -1]);
+    approx.deepEqual(x, math.matrix([-1, -0.5, -1/3, -0.25], 'ccs'));
+
+    var y = math.lusolve(lup, [1, 2, 1, -1]);
+    approx.deepEqual(y, math.matrix([1, 1, 1/3, -0.25], 'ccs'));
+  });
+  
+  it('should solve linear system 3 x 3, no permutations, arrays', function () {
+    var m = 
+        [
+          [2, 1, 1],
+          [1, 2, -1],
+          [1, 2, 1]
+        ];
+    var b = [-2, 4, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, [[-5/3], [7/3], [-1]]);
+  });
+  
+  it('should solve linear system 3 x 3, no permutations, matrix', function () {
+    var m = math.matrix(
+      [
+        [2, 1, 1],
+        [1, 2, -1],
+        [1, 2, 1]
+      ]);
+    var b = [-2, 4, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, math.matrix([[-5/3], [7/3], [-1]]));
+  });
+  
+  it('should solve linear system 3 x 3, no permutations, crs matrix', function () {
+    var m = math.matrix(
+      [
+        [2, 1, 1],
+        [1, 2, -1],
+        [1, 2, 1]
+      ], 'crs');
+    var b = [-2, 4, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, math.matrix([[-5/3], [7/3], [-1]], 'crs'));
+  });
+  
+  it('should solve linear system 3 x 3, no permutations, ccs matrix', function () {
+    var m = math.matrix(
+      [
+        [2, 1, 1],
+        [1, 2, -1],
+        [1, 2, 1]
+      ], 'ccs');
+    var b = [-2, 4, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, math.matrix([[-5/3], [7/3], [-1]], 'ccs'));
+  });
+  
+  it('should solve linear system 3 x 3, permutations, arrays', function () {
+    var m = 
+        [
+          [1, 2, -1],
+          [2, 1, 1],
+          [1, 2, 1]
+        ];
+    var b = [4, -2, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, [[-5/3], [7/3], [-1]]);
+  });
+
+  it('should solve linear system 3 x 3, permutations, matrix', function () {
+    var m = math.matrix(
+      [
+        [1, 2, -1],
+        [2, 1, 1],
+        [1, 2, 1]
+      ]);
+    var b = [4, -2, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, math.matrix([[-5/3], [7/3], [-1]]));
+  });
+
+  it('should solve linear system 3 x 3, permutations, crs matrix', function () {
+    var m = math.matrix(
+      [
+        [1, 2, -1],
+        [2, 1, 1],
+        [1, 2, 1]
+      ], 'crs');
+    var b = [4, -2, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, math.matrix([[-5/3], [7/3], [-1]], 'crs'));
+  });
+
+  it('should solve linear system 3 x 3, permutations, ccs matrix', function () {
+    var m = math.matrix(
+      [
+        [1, 2, -1],
+        [2, 1, 1],
+        [1, 2, 1]
+      ], 'ccs');
+    var b = [4, -2, 2];
+
+    var x = math.lusolve(m, b);
+
+    approx.deepEqual(x, math.matrix([[-5/3], [7/3], [-1]], 'ccs'));
+  });
+  
+  it('should throw exception when matrix is singular', function () {
+    assert.throws(function () { math.lusolve([[1, 1], [0, 0]], [1, 1]); }, /Error: Linear system cannot be solved since matrix is singular/);
+    assert.throws(function () { math.lusolve(math.matrix([[1, 1], [0, 0]], 'dense'), [1, 1]); }, /Error: Linear system cannot be solved since matrix is singular/);
+    assert.throws(function () { math.lusolve(math.matrix([[1, 1], [0, 0]], 'crs'), [1, 1]); }, /Error: Linear system cannot be solved since matrix is singular/);
+    assert.throws(function () { math.lusolve(math.matrix([[1, 1], [0, 0]], 'ccs'), [1, 1]); }, /Error: Linear system cannot be solved since matrix is singular/);
+  });
+});
\ No newline at end of file
diff --git a/test/function/construction/matrix.test.js b/test/function/construction/matrix.test.js
index 8a7692641..b47601bfe 100644
--- a/test/function/construction/matrix.test.js
+++ b/test/function/construction/matrix.test.js
@@ -1,6 +1,5 @@
 // test matrix construction
 var assert = require('assert'),
-    error = require('../../../lib/error/index'),
     math = require('../../../index'),
     matrix = math.matrix;
 
@@ -73,7 +72,7 @@ describe('matrix', function() {
   });
 
   it('should throw an error if called with too many arguments', function() {
-    assert.throws(function () {matrix([], 3, 3)}, /TypeError: Too many arguments/);
+    assert.throws(function () {matrix([], 3, 3);}, /TypeError: Too many arguments/);
   });
 
   it('should throw an error when called with an invalid storage format', function () {
diff --git a/test/type/FibonacciHeap.test.js b/test/type/FibonacciHeap.test.js
new file mode 100644
index 000000000..8e73b7aa6
--- /dev/null
+++ b/test/type/FibonacciHeap.test.js
@@ -0,0 +1,132 @@
+var assert = require('assert');
+var math = require('../../index');
+var FibonacciHeap = math.type.FibonacciHeap;
+
+describe('FibonacciHeap', function () {
+
+  describe('constructor', function () {
+    
+    it('should create heap', function () {
+      var h = new FibonacciHeap();
+      assert.equal(h._size, 0);
+      assert(h._minimum === null);
+    });
+    
+    it('should throw an error when called without new keyword', function () {
+      assert.throws(function () { FibonacciHeap(); }, /Constructor must be called with the new operator/);
+    });
+  });
+  
+  describe('insert', function () {
+    
+    it('should insert node when heap is empty', function () {
+      var h = new FibonacciHeap();
+      h.insert(1, 'v1');
+      assert.equal(h._size, 1);
+      assert(h._minimum !== null);
+      assert.equal(h._minimum.key, 1);
+      assert.equal(h._minimum.value, 'v1');
+    });
+    
+    it('should insert two nodes when heap is empty', function () {
+      var h = new FibonacciHeap();
+      h.insert(1, 'v1');
+      h.insert(10, 'v10');
+      assert.equal(h._size, 2);
+      assert(h._minimum !== null);
+      assert.equal(h._minimum.key, 1);
+      assert.equal(h._minimum.value, 'v1');
+    });
+    
+    it('should insert two nodes when heap is empty, reverse order', function () {
+      var h = new FibonacciHeap();
+      h.insert(10, 'v10');
+      h.insert(1, 'v1');
+      assert.equal(h._size, 2);
+      assert(h._minimum !== null);
+      assert.equal(h._minimum.key, 1);
+      assert.equal(h._minimum.value, 'v1');
+    });
+  });
+  
+  describe('extractMinimum', function () {
+    
+    it('should extract node from heap, one node', function () {
+      var h = new FibonacciHeap();
+      h.insert(1, 'v1');
+      var n = h.extractMinimum();
+      assert.equal(n.key, 1);
+      assert.equal(n.value, 'v1');
+      assert.equal(h._size, 0);
+      assert(h._minimum === null);
+    });
+    
+    it('should extract node from heap, two nodes', function () {
+      var h = new FibonacciHeap();
+      h.insert(1, 'v1');
+      h.insert(10, 'v10');
+      var n = h.extractMinimum();
+      assert.equal(n.key, 1);
+      assert.equal(n.value, 'v1');
+      assert.equal(h._size, 1);
+      assert.equal(h._minimum.key, 10);
+      assert.equal(h._minimum.value, 'v10');
+    });
+    
+    it('should extract nodes in ascending order', function () {
+      var h = new FibonacciHeap();
+      h.insert(5, 'v5');
+      h.insert(4, 'v4');
+      h.insert(1, 'v1');
+      h.insert(3, 'v3');
+      h.insert(2, 'v2');
+      // extract all nodes      
+      var n;
+      var l = h.extractMinimum();
+      var s = h._size;
+      while (true) {
+        n = h.extractMinimum();
+        if (!n)
+          break;
+        assert(n.key > l.key);
+        assert.equal(h._size, --s);
+        l = n;
+      }
+      assert.equal(h._size, 0);
+      assert(h._minimum === null);
+    });
+  });
+  
+  describe('remove', function () {
+    
+    it('should remove node, one node', function () {
+      var h = new FibonacciHeap();
+      var n = h.insert(1, 'v1');
+      h.remove(n);
+      assert.equal(h._size, 0);
+      assert(h._minimum === null);
+    });
+    
+    it('should remove node with smaller key', function () {
+      var h = new FibonacciHeap();
+      h.insert(20, 'v20');
+      var n = h.insert(1, 'v1');
+      h.insert(10, 'v10');
+      h.insert(5, 'v5');
+      h.insert(4, 'v4');
+      h.remove(n);
+      assert.equal(h._size, 4);
+    });
+    
+    it('should remove node with largest key', function () {
+      var h = new FibonacciHeap();      
+      h.insert(1, 'v1');
+      h.insert(10, 'v10');
+      var n = h.insert(20, 'v20');
+      h.insert(5, 'v5');
+      h.insert(4, 'v4');
+      h.remove(n);
+      assert.equal(h._size, 4);
+    });
+  });
+});
\ No newline at end of file
diff --git a/test/type/matrix/CcsMatrix.test.js b/test/type/matrix/CcsMatrix.test.js
index 1a25e2b8d..c92c3e363 100644
--- a/test/type/matrix/CcsMatrix.test.js
+++ b/test/type/matrix/CcsMatrix.test.js
@@ -3,7 +3,6 @@ var math = require('../../../index');
 var index = math.index;
 var Matrix = math.type.Matrix;
 var CcsMatrix = math.type.CcsMatrix;
-var DenseMatrix = math.type.DenseMatrix;
 var Complex = math.type.Complex;
 
 describe('CcsMatrix', function() {
@@ -91,7 +90,7 @@ describe('CcsMatrix', function() {
     });
     
     it('should create a CcsMatrix from a DenseMatrix', function () {
-      var m1 = new DenseMatrix(
+      var m1 = math.matrix(
         [
           [1, 2, 3],
           [4, 5, 6],
@@ -103,6 +102,19 @@ describe('CcsMatrix', function() {
       assert.deepEqual(m1.toArray(), m2.toArray());
     });
     
+    it('should create a CcsMatrix from a CrsMatrix', function () {
+      var m1 = math.matrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [7, 8, 9],
+          [10, 11, 12]
+        ], 'crs');
+      var m2 = new CcsMatrix(m1);
+      assert.deepEqual(m1.size(), m2.size());
+      assert.deepEqual(m1.toArray(), m2.toArray());
+    });
+    
     it('should throw an error when called without new keyword', function () {
       assert.throws(function () { CcsMatrix(); }, /Constructor must be called with the new operator/);
     });
@@ -1948,4 +1960,121 @@ describe('CcsMatrix', function() {
       assert.throws(function () {math.matrix([1,1,1], 'ccs').multiply([[1,1], [1,1]]);});
     });
   });
+  
+  describe('swapRows', function () {
+
+    it('should swap rows with values', function () {
+      var m = new CcsMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [7, 8, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(1, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [7, 8, 9],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap row with value and no values', function () {
+      var m = new CcsMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [0, 0, 0],
+          [10, 11, 12]
+        ]);
+      m.swapRows(1, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [0, 0, 0],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap row with no value and values', function () {
+      var m = new CcsMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [0, 0, 0],
+          [10, 11, 12]
+        ]);
+      m.swapRows(2, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [0, 0, 0],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap rows with missing values', function () {
+      var m = new CcsMatrix(
+        [
+          [1, 2, 3],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(2, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [7, 0, 9],
+          [0, 5, 0],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap last row with another row', function () {
+      var m = new CcsMatrix(
+        [
+          [1, 2, 3],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(3, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [10, 11, 12],
+          [7, 0, 9],
+          [0, 5, 0]
+        ]);
+    });
+
+    it('should swap first row with another row', function () {
+      var m = new CcsMatrix(
+        [
+          [0, 2, 0],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 0, 0]
+        ]);
+      m.swapRows(0, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [7, 0, 9],
+          [0, 5, 0],
+          [0, 2, 0],
+          [10, 0, 0]
+        ]);
+    });
+  });
 });
\ No newline at end of file
diff --git a/test/type/matrix/CrsMatrix.test.js b/test/type/matrix/CrsMatrix.test.js
index ffcb0459a..866a5b824 100644
--- a/test/type/matrix/CrsMatrix.test.js
+++ b/test/type/matrix/CrsMatrix.test.js
@@ -3,7 +3,6 @@ var math = require('../../../index');
 var index = math.index;
 var Matrix = math.type.Matrix;
 var CrsMatrix = math.type.CrsMatrix;
-var DenseMatrix = math.type.DenseMatrix;
 var Complex = math.type.Complex;
 
 describe('CrsMatrix', function() {
@@ -91,7 +90,7 @@ describe('CrsMatrix', function() {
     });
 
     it('should create a CrsMatrix from a DenseMatrix', function () {
-      var m1 = new DenseMatrix(
+      var m1 = math.matrix(
         [
           [1, 2, 3],
           [4, 5, 6],
@@ -102,6 +101,19 @@ describe('CrsMatrix', function() {
       assert.deepEqual(m1.size(), m2.size());
       assert.deepEqual(m1.toArray(), m2.toArray());
     });
+    
+    it('should create a CrsMatrix from a CcsMatrix', function () {
+      var m1 = math.matrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [7, 8, 9],
+          [10, 11, 12]
+        ], 'ccs');
+      var m2 = new CrsMatrix(m1);
+      assert.deepEqual(m1.size(), m2.size());
+      assert.deepEqual(m1.toArray(), m2.toArray());
+    });
 
     it('should throw an error when called without new keyword', function () {
       assert.throws(function () { CrsMatrix(); }, /Constructor must be called with the new operator/);
@@ -1930,4 +1942,121 @@ describe('CrsMatrix', function() {
       assert.throws(function () {math.matrix([1,1,1], 'crs').multiply([[1,1], [1,1]]);});
     });
   });
+  
+  describe('swapRows', function () {
+
+    it('should swap rows with values', function () {
+      var m = new CrsMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [7, 8, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(1, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [7, 8, 9],
+          [4, 5, 6],
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap row with value and no values', function () {
+      var m = new CrsMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [0, 0, 0],
+          [10, 11, 12]
+        ]);
+      m.swapRows(1, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [0, 0, 0],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap row with no value and values', function () {
+      var m = new CrsMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [0, 0, 0],
+          [10, 11, 12]
+        ]);
+      m.swapRows(2, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [0, 0, 0],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap rows with missing values', function () {
+      var m = new CrsMatrix(
+        [
+          [1, 2, 3],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(2, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [7, 0, 9],
+          [0, 5, 0],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap last row with another row', function () {
+      var m = new CrsMatrix(
+        [
+          [1, 2, 3],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(3, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [10, 11, 12],
+          [7, 0, 9],
+          [0, 5, 0]
+        ]);
+    });
+
+    it('should swap first row with another row', function () {
+      var m = new CrsMatrix(
+        [
+          [0, 2, 0],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 0, 0]
+        ]);
+      m.swapRows(0, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [7, 0, 9],
+          [0, 5, 0],
+          [0, 2, 0],
+          [10, 0, 0]
+        ]);
+    });
+  });
 });
\ No newline at end of file
diff --git a/test/type/matrix/DenseMatrix.test.js b/test/type/matrix/DenseMatrix.test.js
index bab8985b7..e10b7b8cd 100644
--- a/test/type/matrix/DenseMatrix.test.js
+++ b/test/type/matrix/DenseMatrix.test.js
@@ -1268,6 +1268,123 @@ describe('DenseMatrix', function() {
     });
   });
   
+  describe('swapRows', function () {
+
+    it('should swap rows with values', function () {
+      var m = new DenseMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [7, 8, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(1, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [7, 8, 9],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap row with value and no values', function () {
+      var m = new DenseMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [0, 0, 0],
+          [10, 11, 12]
+        ]);
+      m.swapRows(1, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [0, 0, 0],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap row with no value and values', function () {
+      var m = new DenseMatrix(
+        [
+          [1, 2, 3],
+          [4, 5, 6],
+          [0, 0, 0],
+          [10, 11, 12]
+        ]);
+      m.swapRows(2, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [0, 0, 0],
+          [4, 5, 6],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap rows with missing values', function () {
+      var m = new DenseMatrix(
+        [
+          [1, 2, 3],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(2, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [7, 0, 9],
+          [0, 5, 0],          
+          [10, 11, 12]
+        ]);
+    });
+
+    it('should swap last row with another row', function () {
+      var m = new DenseMatrix(
+        [
+          [1, 2, 3],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 11, 12]
+        ]);
+      m.swapRows(3, 1);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [1, 2, 3],
+          [10, 11, 12],
+          [7, 0, 9],
+          [0, 5, 0]
+        ]);
+    });
+
+    it('should swap first row with another row', function () {
+      var m = new DenseMatrix(
+        [
+          [0, 2, 0],
+          [0, 5, 0],
+          [7, 0, 9],
+          [10, 0, 0]
+        ]);
+      m.swapRows(0, 2);
+      assert.deepEqual(
+        m.valueOf(),
+        [
+          [7, 0, 9],
+          [0, 5, 0],
+          [0, 2, 0],
+          [10, 0, 0]
+        ]);
+    });
+  });
+  
   /**
    * Helper function to create an Array containing uninitialized values
    * Example: arr(uninit, uninit, 2);    // [ , , 2 ]
diff --git a/test/type/matrix/Spa.test.js b/test/type/matrix/Spa.test.js
new file mode 100644
index 000000000..19ba4fc39
--- /dev/null
+++ b/test/type/matrix/Spa.test.js
@@ -0,0 +1,175 @@
+var assert = require('assert');
+var math = require('../../../index');
+var Spa = math.type.Spa;
+
+describe('Spa', function() {
+
+  describe('constructor', function() {
+    
+    it('Should create Spa with initial length', function () {
+      var spa = new Spa(10);
+      assert(spa._values !== null);
+      assert(spa._heap !== null);
+      assert.equal(spa._values.length, 10);
+    });
+    
+    it('should throw an error when called without new keyword', function () {
+      assert.throws(function () { Spa(); }, /Constructor must be called with the new operator/);
+    });
+  });
+  
+  describe('test', function () {
+    
+    it('should add value when no value exists for row', function () {
+      var spa = new Spa(10);
+      spa.set(5, 0.5);
+      assert(spa._values[5]);
+      assert(spa._values[5].value);
+      assert.equal(spa._values[5].value, 0.5);
+    });
+    
+    it('should set value when value exists for row', function () {
+      var spa = new Spa(10);
+      spa.set(5, 0.5);
+      assert(spa._values[5]);
+      assert(spa._values[5].value);
+      assert.equal(spa._values[5].value, 0.5);
+      spa.set(5, 1.5);
+      assert(spa._values[5]);
+      assert(spa._values[5].value);
+      assert.equal(spa._values[5].value, 1.5);
+    });
+  });
+  
+  describe('get', function () {
+    
+    it('should get zero when no value exists for row', function () {
+      var spa = new Spa(10);
+      var v = spa.get(5);
+      assert.equal(v, 0);
+    });
+    
+    it('should get value when value exists for row', function () {
+      var spa = new Spa(10);
+      spa.set(5, 0.5);
+      var v = spa.get(5);
+      assert.equal(v, 0.5);
+    });
+  });
+  
+  describe('accumulate', function () {
+    
+    it('should add value when no value exists for row', function () {
+      var spa = new Spa(10);
+      spa.accumulate(5, 0.5);
+      assert(spa._values[5]);
+      assert(spa._values[5].value);
+      assert.equal(spa._values[5].value, 0.5);
+    });
+    
+    it('should accumulate value when value exists for row', function () {
+      var spa = new Spa(10);
+      spa.set(5, 0.5);
+      spa.accumulate(5, 1.5);
+      assert(spa._values[5]);
+      assert(spa._values[5].value);
+      assert.equal(spa._values[5].value, 2);
+    });
+  });
+  
+  describe('swap', function () {
+    
+    it('should swap two existing values', function () {
+      var spa = new Spa(10);
+      spa.set(5, 0.5);
+      spa.set(2, 0.2);
+      spa.swap(2, 5);
+      assert(spa._values[5]);
+      assert(spa._values[5].value);
+      assert.equal(spa._values[5].value, 0.2);
+      assert(spa._values[2]);
+      assert(spa._values[2].value);
+      assert.equal(spa._values[2].value, 0.5);
+    });
+    
+    it('should swap an existing and non existing values', function () {
+      var spa = new Spa(10);
+      spa.set(5, 0.5);
+      spa.swap(5, 2);
+      assert(!spa._values[5]);
+      assert(spa._values[2]);
+      assert(spa._values[2].value);
+      assert.equal(spa._values[2].value, 0.5);
+    });
+    
+    it('should swap a non existing and existing values', function () {
+      var spa = new Spa(10);
+      spa.set(5, 0.5);
+      spa.swap(2, 5);
+      assert(!spa._values[5]);
+      assert(spa._values[2]);
+      assert(spa._values[2].value);
+      assert.equal(spa._values[2].value, 0.5);
+    });
+    
+    it('should swap two non existing values', function () {
+      var spa = new Spa(10);
+      spa.swap(2, 5);
+      assert(!spa._values[5]);
+      assert(!spa._values[2]);
+    });
+  });
+  
+  describe('forEach', function () {
+    
+    it('should enumerate values in correct order', function () {
+      var spa = new Spa(10);
+      spa.set(2, 2);
+      spa.set(3, 3);
+      spa.set(4, 4);
+      spa.set(1, 2);
+      var x;
+      var c = 0;
+      spa.forEach(0, 9, function (i) {
+        if (!x) {
+          assert.equal(i, 1);
+          x = i;
+        }
+        else {
+          assert(i > x);
+          x = i;
+        }
+        c++;
+      });
+      assert.equal(c, 4);
+      assert.equal(x, 4);
+      assert(spa._heap._minimum !== null);
+      assert.equal(spa._heap._size, 4);
+    });
+    
+    it('should enumerate values in interval', function () {
+      var spa = new Spa(10);
+      spa.set(2, 2);
+      spa.set(3, 3);
+      spa.set(4, 4);
+      spa.set(1, 2);
+      var x;
+      var c = 0;
+      spa.forEach(2, 3, function (i) {
+        if (!x) {
+          assert.equal(i, 2);
+          x = i;
+        }
+        else {
+          assert(i > x);
+          x = i;
+        }
+        c++;
+      });
+      assert.equal(c, 2);
+      assert.equal(x, 3);
+      assert(spa._heap._minimum !== null);
+      assert.equal(spa._heap._size, 4);
+    });
+  });
+});