From 019c9c0c5fc1f61dee93442b2c82926e7676ba21 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vladim=C3=ADr=20Vondru=C5=A1?= Date: Mon, 4 Feb 2013 18:46:37 +0100 Subject: [PATCH] Math: Gauss-Jordan algorithm as pure functions. No need to wrap this in a class. --- src/Math/Algorithms/GaussJordan.h | 79 +++++++++----------- src/Math/Algorithms/Test/GaussJordanTest.cpp | 4 +- 2 files changed, 36 insertions(+), 47 deletions(-) diff --git a/src/Math/Algorithms/GaussJordan.h b/src/Math/Algorithms/GaussJordan.h index dc1734346..130403216 100644 --- a/src/Math/Algorithms/GaussJordan.h +++ b/src/Math/Algorithms/GaussJordan.h @@ -16,7 +16,7 @@ */ /** @file - * @brief Class Magnum::Math::Algorithms::GaussJordan + * @brief Function Magnum::Math::Algorithms::gaussJordanInPlaceTransposed(), Magnum::Math::Algorithms::gaussJordanInPlace() */ #include "Math/RectangularMatrix.h" @@ -24,54 +24,25 @@ namespace Magnum { namespace Math { namespace Algorithms { /** -@brief Gauss-Jordan elimination +@brief In-place Gauss-Jordan elimination on transposed matrices +@param a Transposed left side of augmented matrix +@param t Transposed right side of augmented matrix +@return True if @p a is regular, false if @p a is singular (and thus the + system cannot be solved). + +As Gauss-Jordan elimination works on rows and matrices in OpenGL are +column-major, it is more efficient to operate on transposed matrices and treat +columns as rows. See also gaussJordanInPlace() which works with non-transposed matrices. + +The function eliminates matrix @p a and solves @p t in place. For efficiency +reasons, only pure Gaussian elimination is done on @p a and the final +backsubstitution is done only on @p t, as @p a would always end with identity +matrix anyway. Based on ultra-compact Python code by Jarno Elonen, http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html. */ -class GaussJordan { - public: - GaussJordan() = delete; - - /** - * @brief Eliminate transposed matrices in place - * @param a Transposed left side of augmented matrix - * @param t Transposed right side of augmented matrix - * @return True if @p a is regular, false if @p a is singular (and - * thus the system cannot be solved). - * - * As Gauss-Jordan elimination works on rows and matrices in OpenGL - * are column-major, it is more efficient to operate on transposed - * matrices and treat columns as rows. See also inPlace() which works - * with non-transposed matrices. - * - * The function eliminates matrix @p a and solves @p t in place. For - * efficiency reasons, only pure Gaussian elimination is done on @p a - * and the final backsubstitution is done only on @p t, as @p a would - * always end with identity matrix anyway. - */ - template static bool inPlaceTransposed(RectangularMatrix& a, RectangularMatrix& t); - - /** - * @brief Eliminate in place - * - * Transposes the matrices, calls inPlaceTransposed() on them and then - * transposes them back. - */ - template static bool inPlace(RectangularMatrix& a, RectangularMatrix& t) { - a = a.transposed(); - RectangularMatrix tTransposed = t.transposed(); - - bool ret = inPlaceTransposed(a, tTransposed); - - a = a.transposed(); - t = tTransposed.transposed(); - - return ret; - } -}; - -template bool GaussJordan::inPlaceTransposed(RectangularMatrix& a, RectangularMatrix& t) { +template bool gaussJordanInPlaceTransposed(RectangularMatrix& a, RectangularMatrix& t) { for(std::size_t row = 0; row != size; ++row) { /* Find max pivot */ std::size_t rowMax = row; @@ -110,6 +81,24 @@ template bool GaussJordan::inPlaceT return true; } +/** +@brief In-place Gauss-Jordan elimination + +Transposes the matrices, calls gaussJordanInPlaceTransposed() on them and then +transposes them back. +*/ +template bool gaussJordanInPlace(RectangularMatrix& a, RectangularMatrix& t) { + a = a.transposed(); + RectangularMatrix tTransposed = t.transposed(); + + bool ret = gaussJordanInPlaceTransposed(a, tTransposed); + + a = a.transposed(); + t = tTransposed.transposed(); + + return ret; +} + }}} #endif diff --git a/src/Math/Algorithms/Test/GaussJordanTest.cpp b/src/Math/Algorithms/Test/GaussJordanTest.cpp index 1825769df..155144f84 100644 --- a/src/Math/Algorithms/Test/GaussJordanTest.cpp +++ b/src/Math/Algorithms/Test/GaussJordanTest.cpp @@ -42,7 +42,7 @@ void GaussJordanTest::singular() { Vector4(1.0f, 2.0f, 7.0f, 40.0f)); RectangularMatrix<4, 1, float> t; - CORRADE_VERIFY(!GaussJordan::inPlaceTransposed(a, t)); + CORRADE_VERIFY(!gaussJordanInPlaceTransposed(a, t)); } void GaussJordanTest::invert() { @@ -58,7 +58,7 @@ void GaussJordanTest::invert() { Matrix4 a2(a); Matrix4 inverse = Matrix4::fromDiagonal(Vector4(1.0f)); - CORRADE_VERIFY(GaussJordan::inPlace(a2, inverse)); + CORRADE_VERIFY(gaussJordanInPlace(a2, inverse)); CORRADE_COMPARE(inverse, expectedInverse); CORRADE_COMPARE(a*inverse, Matrix4::fromDiagonal(Vector4(1.0f)));