LU complete, added some docs

README.md

@@ -10,3 +10,66 @@ Prerequisites
 This library require LAPACK to be built and installed as a shared library.
 In time the entire build process may be unified into this project, but that's
 some time away.
+
+Usage
+=====
+
+    var lapack = require('lapack');
+
+    /*
+    LAPACK functions
+    */
+
+    var result = lapack.sgeqrf([
+        [1, 2, 3],
+        [3, 4, 5],
+        [5, 6, 7]
+    ]);
+
+    console.log(result.R);
+    console.log(result.tau);
+
+    result = sgesvd('A', 'A', [
+        [1, 2, 3],
+        [3, 4, 5],
+        [5, 6, 7]
+    ]);
+
+    console.log(result.U);
+    console.log(result.S);
+    console.log(result.VT);
+
+    result = lapack.sgetrf([
+        [1, 2, 3],
+        [3, 4, 5],
+        [5, 6, 7]
+    ]);
+
+    console.log(result.LU);
+    console.log(result.IPIV);
+
+    /*
+    conveniently wrapped processes
+    */
+
+    // perform a complete LU factorization
+    var lu = lapack.lu([
+        [1, 2, 3],
+        [3, 4, 5],
+        [5, 6, 7]
+    ]);
+
+    console.log(lu.L);
+    console.log(lu.U);
+    console.log(lu.P);
+
+
+    // perform a complete QR factorization
+    var qr = lapack.qr([
+        [1, 2, 3],
+        [3, 4, 5],
+        [5, 6, 7]
+    ]);
+
+    console.log(qr.Q);
+    console.log(qr.R);

lib/node-lapack/index.js

@@ -3,4 +3,5 @@ var lapack = require('./lapack.js');
 exports.sgeqrf = lapack.sgeqrf;
 exports.sgesvd = lapack.sgesvd;
 exports.qr = lapack.qr;
+exports.lu = lapack.lu;
 exports.sgetrf = lapack.sgetrf;

lib/node-lapack/lapack.js

@@ -25,6 +25,24 @@ var FORTRAN_INT = 4;
 var FORTRAN_CHAR = 1;
 var FORTRAN_FLOAT = 4;
 
+function eye(m) {
+    var matrix = [];
+
+    for(var i = 0; i < m; i++) {
+	var row = [];
+	matrix.push(row);
+
+	for(var j = 0; j < m; j++) {
+	    if(i == j)
+		row.push(1);
+	    else
+		row.push(0);
+	}
+    }
+
+    return matrix;
+}
+
 function matrixOp(matrix, callback) {
     var m = matrix.length;
     var n = matrix[0].length;    
@@ -40,21 +58,25 @@ function matrixOp(matrix, callback) {
     callback(m, n, f_m, f_n, f_a, f_lda);
 }
 
+function zeroBottomLeft(matrix) {
+    // zero out bottom left forming an upper right triangle matrix                                    
+    for(var i = 1; i < matrix.length; i++) {
+        for(var j = 0; j < i; j++)
+            matrix[i][j] = 0;
+    }
+
+    return matrix
+}
+
 function qr(matrix) {
     var result;
 
     sgeqrf(matrix, function(qr, m, n, f_m, f_n, f_a, f_lda, f_tau, f_work, f_lwork, f_info) {
 	var f_k = new FFI.Pointer(FORTRAN_INT);
 	f_k.putInt32(Math.min(m, n));
 	LAPACK.SORGQR(f_m, f_n, f_k, f_a, f_lda, f_tau, f_work, f_lwork, f_info);
-	qr['Q'] = fortranArray.fortranArrayToJSMatrix(f_a, m, n);
-
-	// zero out bottom left forming an upper right triangle matrix
-	for(var i = 1; i < m; i++) {
-	    for(var j = 0; j < i; j++)
-		qr.R[i][j] = 0;
-	}
-
+	qr.Q = fortranArray.fortranArrayToJSMatrix(f_a, m, n);
+	qr.R = zeroBottomLeft(qr.R);
 	result = qr;
     });
 
@@ -95,6 +117,55 @@ function sgeqrf(matrix, callback) {
     return qr;
 }
 
+function cloneMatrix(matrix) {
+    var clone = [];
+
+    for(var i = 0; i < matrix.length; i++) {
+	var row = [];
+	clone.push(row);
+
+	for(var j = 0; j < matrix.length; j++) {
+	    row.push(matrix[i][j]);
+	} 
+    } 
+
+    return clone;
+}
+
+function swapRows(matrix, i, j) {
+    var tmp = matrix[j];
+    matrix[j] = matrix[i];
+    matrix[i] = tmp;
+    return matrix;
+}
+
+function lu(matrix) {
+    var result = sgetrf(matrix);
+    var L = cloneMatrix(result.LU);
+    var P = eye(matrix.length);
+
+    for(var i = 0; i < L.length; i++) {
+	for(var j = i; j < L[i].length; j++) {
+	    if(i == j)
+		L[i][j] = 1;
+	    else
+		L[i][j] = 0;
+	}
+    }
+
+    for(var i = 0; i < result.IPIV.length; i++) {
+	if(i != result.IPIV[i] - 1)
+	    swapRows(P, i, result.IPIV[i] - 1);
+    }
+
+    return {
+	L: L,
+	U: zeroBottomLeft(cloneMatrix(result.LU)),
+	P: P,
+	IPIV: result.IPIV
+    };
+}
+
 function sgetrf(matrix) {
     var result = {};
 
@@ -157,3 +228,4 @@ exports.sgeqrf = sgeqrf;
 exports.sgesvd = sgesvd;
 exports.sgetrf = sgetrf;
 exports.qr = qr;
+exports.lu = lu;

package.json

@@ -1,7 +1,7 @@
 {
   "name": "lapack",
   "description": "",
-  "version": "0.0.3",
+  "version": "0.0.4",
   "homepage": "https://github.com/NaturalNode/node-lapack",
   "engines": {
     "node": ">=0.2.6"

spec/lapack_spec.js

@@ -9,14 +9,40 @@ describe('lapack', function() {
         [1, 1, 3]
     ];
 
+    var luIn = [
+	[11, 3, 11, 3],
+	[ 4, 2,  1, 4],
+	[-4, 5,  3, 1],
+	[-9, 4,  3, 9]
+    ];
+
     it('should lu', function() {
-	var luIn = [
-	    [ 4, 2, 1, 4],
-	    [-9, 4, 3, 9],
-	    [11, 3, 11, 3],
-	    [-4, 5, 3, 1]
-	];
-
+	var lu = lapack.lu(luIn);
+	expect(approxEql(lu.L, ([
+	    [ 1, 0, 0, 0 ],
+	    [ -0.8181818723678589, 1, 0, 0 ],
+	    [ 0.3636363744735718, 0.14084506034851074, 1, 0 ],
+	    [ -0.3636363744735718, 0.9436619877815247,
+	      0.9219222664833069, 1 ]	    
+	]))).toBeTruthy();
+
+	expect(approxEql(lu.U, [
+	    [ 11, 3, 11, 3 ],
+	    [ 0, 6.454545497894287,
+	      12.000000953674316, 11.454545974731445 ],
+	    [ 0, 0, -4.690140724182129, 1.2957748174667358 ],
+	    [ 0, 0, 0, -9.912914276123047 ]	    
+	])).toBeTruthy();
+
+	expect(approxEql(lu.P, [
+	    [ 1, 0, 0, 0 ],
+	    [ 0, 0, 0, 1 ],
+	    [ 0, 1, 0, 0 ],
+	    [ 0, 0, 1, 0 ]	    
+	])).toBeTruthy();
+    });
+
+    it('should sgetrf', function() {
 	var result = lapack.sgetrf(luIn);
 	expect(approxEql(result.LU, [ [ 11, 3, 11, 3 ],
 				      [ -0.8181818723678589,