Permalink
Browse files

Fix gauss-jordan bug.

Patch submission from Martin Ossmann.
  • Loading branch information...
1 parent 8dd7aa2 commit a3ff63798ce61c1d2a02c473a2376a7c083e536d @davidm committed Apr 17, 2012
Showing with 46 additions and 26 deletions.
  1. +4 −0 CHANGES.txt
  2. +1 −0 MANIFEST
  3. +6 −0 README.txt
  4. +20 −26 lua/matrix.lua
  5. +15 −0 test/test_matrix.lua
View
@@ -0,0 +1,4 @@
+0.2.11.20120416
+ Gauss-Jordan bug fix. Also affects matrix inverse.
+
+0.2.10.20111203
View
@@ -1,3 +1,4 @@
+CHANGES.txt
README.txt
LICENSE.txt
dist.info
View
@@ -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.
View
@@ -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
View
@@ -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}}

0 comments on commit a3ff637

Please sign in to comment.