From d8a43ffd3f4968ec16c51d5be6ee6438c5498e99 Mon Sep 17 00:00:00 2001 From: David Manura Date: Mon, 16 Apr 2012 23:32:45 -0400 Subject: [PATCH] Fix gauss-jordan bug. Patch submission from Martin Ossmann. --- README.txt | 6 ++++++ dist.info | 2 +- doc/fit_changelog.txt | 7 ------ doc/matrix_changelog.txt | 5 ++++- lua/matrix.lua | 46 +++++++++++++++++----------------------- test/test_matrix.lua | 15 +++++++++++++ 6 files changed, 46 insertions(+), 35 deletions(-) delete mode 100644 doc/fit_changelog.txt diff --git a/README.txt b/README.txt index a0e8b2b..a0c86b5 100644 --- a/README.txt +++ b/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. diff --git a/dist.info b/dist.info index f2426f1..37b5730 100644 --- a/dist.info +++ b/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" diff --git a/doc/fit_changelog.txt b/doc/fit_changelog.txt deleted file mode 100644 index 69b62af..0000000 --- a/doc/fit_changelog.txt +++ /dev/null @@ -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 \ No newline at end of file diff --git a/doc/matrix_changelog.txt b/doc/matrix_changelog.txt index abae275..43940cc 100644 --- a/doc/matrix_changelog.txt +++ b/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 diff --git a/lua/matrix.lua b/lua/matrix.lua index 94eeca1..5e587a3 100644 --- a/lua/matrix.lua +++ b/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 diff --git a/test/test_matrix.lua b/test/test_matrix.lua index 021faf7..38995fd 100644 --- a/test/test_matrix.lua +++ b/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}}