Fix gauss-jordan bug.

Patch submission from Martin Ossmann.

README.txt

@@ -17,6 +17,12 @@ luamatrix").
 
 http://lua-users.org/wiki/LuaMatrix
 
+== Credits ==
+
+Authors: Michael Lutz (original author), David Manura (maintainer).
+Bug fixes/comments: Geoff Richards, SatheeshJM, Martin Ossmann,
+and others.
+
 == License ==
 
 See the LICENSE.txt file for licensing details.

dist.info

@@ -1,7 +1,7 @@
 --- This file is part of LuaDist project
 
 name = "lua-matrix"
-version = "0.2.10"
+version = "0.2.11.20120416"
 
 desc = "'LuaMatrix' provides a good selection of matrix functions."
 author = "Michael Lutz, David Manura"

doc/fit_changelog.txt

@@ -1,7 +0,0 @@
-fit changelog:
-
-v 0.2: 2007-08-26
-	- some optimising, for the input
-
-v 0.1: 2007-08-11
-	- built from former fit

doc/matrix_changelog.txt

@@ -1,6 +1,9 @@
 matrix changelog
 
-v 0.2.10
+v 0.2.11.20120416
+	- Gauss-Jordan bug fix.  Also affects matrix inverse.
+
+v 0.2.10.20111203
 	- Add _VERSION to modules.
 
 v 0.2.9: 2008-08-26

lua/matrix.lua

@@ -114,7 +114,7 @@ LICENSE
 --// matrix //
 --////////////
 
-local matrix = {_TYPE='module', _NAME='matrix', _VERSION='0.2.10.20111203'}
+local matrix = {_TYPE='module', _NAME='matrix', _VERSION='0.2.11.20120416'}
 
 -- access to the metatable we set at the end of the file
 local matrix_meta = {}
@@ -416,35 +416,29 @@ end
 -- returns on failure: false,'rank of matrix'
 
 -- locals
--- checking here for the nearest element to 1 or -1; (smallest pivot element)
--- this way the factor of the evolving number division should be > 1 or the
--- divided number itself, what gives better results
-local setelementtosmallest = function( mtx,i,j,zero,one,norm2 )
-	-- check if element is one
-	if mtx[i][j] == one then return true end
-	-- check for lowest value
-	local _ilow
+-- 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.
+local pivotOk = function( mtx,i,j,norm2 )
+	-- find min value
+	local iMin
+	local normMin = math.huge
 	for _i = i,#mtx do
 		local e = mtx[_i][j]
-		if e == one then
-			break
-		end
-		if not _ilow then
-			if e ~= zero then
-				_ilow = _i
+		local norm = math.abs(norm2(e))
+		if norm > 0 and norm < normMin then
+			iMin = _i
+			normMin = norm
 			end
-		elseif (e ~= zero) and math.abs(norm2(e)-1) < math.abs(norm2(mtx[_ilow][j])-1) then
-			_ilow = _i
 		end
-	end
-	if _ilow then
-		-- switch lines if not input line
-		-- legal operation
-		if _ilow ~= i then
-			mtx[i],mtx[_ilow] = mtx[_ilow],mtx[i]
+	if iMin then
+		-- switch lines if not in position.
+		if iMin ~= i then
+			mtx[i],mtx[iMin] = mtx[iMin],mtx[i]
 		end
 		return true
-	end
+		end
+	return false
 end
 
 local function copy(x)
@@ -463,7 +457,7 @@ function matrix.dogauss( mtx )
 	-- stairs left -> right
 	for j = 1,rows do
 		-- check if element can be setted to one
-		if setelementtosmallest( mtx,j,j,zero,one,norm2 ) then
+		if pivotOk( mtx,j,j,norm2 ) then
 			-- start parsing rows
 			for i = j+1,rows do
 				-- check if element is not already zero
@@ -540,7 +534,7 @@ function matrix.invert( m1 )
 end
 
 --// matrix.sqrt ( m1 [,iters] )
--- calculate the square root of a matrix using "Denman–Beavers square root iteration"
+-- calculate the square root of a matrix using "Denman Beavers square root iteration"
 -- condition: matrix rows == matrix columns; must have a invers matrix and a square root
 -- if called without additional arguments, the function finds the first nearest square root to
 -- input matrix, there are others but the error between them is very small

test/test_matrix.lua

@@ -135,6 +135,21 @@ mtxinvcomp = matrix( mtxinvcomp )
 mtxinv:round( 5 )
 mtxinv = mtxinv:elementstostrings()
 assert( mtxinvcomp == mtxinv  )
+-- Fixed in v0.2.11 failed (Gauss-Jordan)
+local mtx = matrix {{ 1, -1,  1 }, {-1,  1,  0 }, { 1,  0,  0 } }
+assert((mtx^-1)^-1 == mtx)
+-- Fixed in v0.2.11 failed (Gauss-Jordan)
+local mtx = matrix {{ 0,  0,  1 }, { 0,  1,  0 }, { 1,  0,  0 } }
+assert((mtx^-1)^-1 == mtx)
+-- similar to above but with complex
+local mtx = matrix.replace({{ 0,  0,  "1i" }, { 0,  "1i",  0 }, { "1i",  0,  0 } }, complex)
+assert((mtx^-1)^-1 == mtx)
+-- random
+for i=1,10 do
+  local mtx = matrix(4, 4):random(-20, 20, 5)
+  assert(matrix.normmax((mtx^-1)^-1 - mtx) < 1E-13)
+end
+
 
 -- matrix.sqrt; number
 local m1 = matrix{{4,2,1},{1,5,4},{1,5,2}}