From 88cc0f2dd152b2163f81267beda4ee632e4d0a94 Mon Sep 17 00:00:00 2001
From: Waylon Flinn <waylonflinn@gmail.com>
Date: Thu, 20 Feb 2014 09:04:55 -0600
Subject: [PATCH 1/2] changed choice of pivot in guass-jordan method.

new code uses maximum instead of minimum. based on code and discussion in
NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
---
 lua/matrix.lua | 37 ++++++++++++++++++++++---------------
 1 file changed, 22 insertions(+), 15 deletions(-)

diff --git a/lua/matrix.lua b/lua/matrix.lua
index 83f7d57..83cf024 100644
--- a/lua/matrix.lua
+++ b/lua/matrix.lua
@@ -415,30 +415,37 @@ end
 -- returns on success: true,
 -- returns on failure: false,'rank of matrix'
 
+
 -- locals
--- checking here for the element nearest but not equal to zero (smallest pivot element).
--- This way the `factor` in `dogauss` will be >= 1, which
--- can give better results.
+-- perform an in-place pivot on the given matrix, if possible,
+-- returning whether or not a pivot was done
+
+-- we choose the pivot element by checking for the largest element in a column,
+-- with an index greater than a specified starting point
+-- mtx - the matrix to pivot
+-- i - the row to start scan for pivot on
+-- j - the column to scan for a pivot in
+-- norm2 - a function which gives the square of the norm for a matrix element
 local pivotOk = function( mtx,i,j,norm2 )
-	-- find min value
-	local iMin
-	local normMin = math.huge
+	-- find max value
+	local iMax
+	local normMax = 0
 	for _i = i,#mtx do
 		local e = mtx[_i][j]
 		local norm = math.abs(norm2(e))
-		if norm > 0 and norm < normMin then
-			iMin = _i
-			normMin = norm
+		if norm > normMax then
+			iMax = _i
+			normMax = norm
 			end
 		end
-	if iMin then
-		-- switch lines if not in position.
-		if iMin ~= i then
-			mtx[i],mtx[iMin] = mtx[iMin],mtx[i]
+	if iMax then
+		-- switch rows, if not in position.
+		if iMax ~= i then
+			mtx[i],mtx[iMax] = mtx[iMax],mtx[i]
 		end
 		return true
-		end
-	return false
+	end
+	return false -- column was all zeros below pivot
 end
 
 local function copy(x)

From d89d9df264654d31ea221051644330274de62e8c Mon Sep 17 00:00:00 2001
From: Waylon Flinn <waylonflinn@gmail.com>
Date: Thu, 20 Feb 2014 10:07:11 -0600
Subject: [PATCH 2/2] fixed floating point comparison in tests.

---
 test/test_matrix.lua | 23 ++++++++++++++++++++---
 1 file changed, 20 insertions(+), 3 deletions(-)

diff --git a/test/test_matrix.lua b/test/test_matrix.lua
index 38995fd..67fc7a9 100644
--- a/test/test_matrix.lua
+++ b/test/test_matrix.lua
@@ -4,6 +4,23 @@ local symbol = matrix.symbol
 
 local mtx, m1,m2,m3,m4,m5, ms,ms1,ms2,ms3,ms4
 
+-- compare matrices of floats using epsilon = 1e-9
+-- not perfect, but better than '=='
+-- improvements can be found here: http://stackoverflow.com/a/3423299/74291
+function matrix.fcompare(m1, m2)
+	local eps = 0.000000001
+	mtx = {}
+	for i = 1,#m1 do
+		mtx[i] = {}
+		for j = 1,#m1[1] do
+			if math.abs(m1[i][j] - m2[i][j]) > eps then
+				return false
+			end
+		end
+	end
+	return true
+end
+
 -- test matrix:new/matrix call function
 -- normal matrix
 mtx = matrix {{1,2},{3,4}}
@@ -72,7 +89,7 @@ assert( m1*m2^-1 == m1/m2 )
 m2 = matrix {{4,5},{6,7}}
 assert( m2/2 == matrix{{2,2.5},{3,3.5}} )
 mtx = matrix {{3,5,1},{2,4,5},{1,2,2}}
-assert( 2 / mtx == matrix{{4,16,-42},{-2,-10,26},{0,2,-4}} )
+assert( matrix.fcompare( 2 / mtx, matrix{{4,16,-42},{-2,-10,26},{0,2,-4}} ) )
 -- matrix.mulnum; symbol
 m1 = m1:replace(symbol)
 assert( m1*2 == matrix{{"(1)*(2)","(2)*(2)"},{"(3)*(2)","(4)*(2)"}}:replace(symbol) )
@@ -114,11 +131,11 @@ assert( mtx:det():round(10) == complex "5527+2687i" )
 --1
 mtx = matrix{{3,5,1},{2,4,5},{1,2,2}}
 local mtxinv = matrix{{2,8,-21},{-1,-5,13},{0,1,-2}}
-assert( mtx^-1 == mtxinv )
+assert( matrix.fcompare( mtx^-1, mtxinv) )
 --2
 mtx = matrix{{1,0,2},{4,1,1},{3,2,-7}} 
 local mtxinv = matrix{{-9,4,-2},{31,-13,7},{5,-2,1}}
-assert( mtx^-1 == mtxinv )
+assert( matrix.fcompare( mtx^-1, mtxinv) )
 -- matrix.invert; complex
 mtx = {
 { 3,"4-3i",1},