From 1727357cc5b6e3030ae46c916b657c94874e4f84 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vladim=C3=ADr=20Vondru=C5=A1?= Date: Sun, 26 Aug 2012 22:43:33 +0200 Subject: [PATCH] Gauss-Jordan elimination. --- doc/namespaces.dox | 9 ++ src/Math/Algorithms/CMakeLists.txt | 8 ++ src/Math/Algorithms/GaussJordan.h | 113 +++++++++++++++++++ src/Math/Algorithms/Test/CMakeLists.txt | 1 + src/Math/Algorithms/Test/GaussJordanTest.cpp | 60 ++++++++++ src/Math/Algorithms/Test/GaussJordanTest.h | 32 ++++++ src/Math/CMakeLists.txt | 1 + 7 files changed, 224 insertions(+) create mode 100644 src/Math/Algorithms/CMakeLists.txt create mode 100644 src/Math/Algorithms/GaussJordan.h create mode 100644 src/Math/Algorithms/Test/CMakeLists.txt create mode 100644 src/Math/Algorithms/Test/GaussJordanTest.cpp create mode 100644 src/Math/Algorithms/Test/GaussJordanTest.h diff --git a/doc/namespaces.dox b/doc/namespaces.dox index 63587bca1..faf2cb178 100644 --- a/doc/namespaces.dox +++ b/doc/namespaces.dox @@ -28,6 +28,15 @@ Base classes for creating OpenGL contexts with various toolkits. Template classes for matrix and vector calculations. */ +/** @dir Math/Algorithms + * @brief Namespace Magnum::Math::Algorithms + */ +/** @namespace Magnum::Math::Algorithms +@brief %Algorithms + +Various matrix and vector algorithms. +*/ + /** @dir Math/Geometry * @brief Namespace Magnum::Math::Geometry */ diff --git a/src/Math/Algorithms/CMakeLists.txt b/src/Math/Algorithms/CMakeLists.txt new file mode 100644 index 000000000..cb5121b20 --- /dev/null +++ b/src/Math/Algorithms/CMakeLists.txt @@ -0,0 +1,8 @@ +set(MagnumMathAlgorithms_HEADERS + GaussJordan.h) +install(FILES ${MagnumMathAlgorithms_HEADERS} DESTINATION ${MAGNUM_INCLUDE_INSTALL_DIR}/Math/Algorithms) + +if(BUILD_TESTS) + enable_testing() + add_subdirectory(Test) +endif() diff --git a/src/Math/Algorithms/GaussJordan.h b/src/Math/Algorithms/GaussJordan.h new file mode 100644 index 000000000..d98bf89c8 --- /dev/null +++ b/src/Math/Algorithms/GaussJordan.h @@ -0,0 +1,113 @@ +#ifndef Magnum_Math_Algorithms_GaussJordan_h +#define Magnum_Math_Algorithms_GaussJordan_h +/* + Copyright © 2010, 2011, 2012 Vladimír Vondruš + + This file is part of Magnum. + + Magnum is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 3 + only, as published by the Free Software Foundation. + + Magnum is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License version 3 for more details. +*/ + +/** @file + * @brief Class GaussJordan + */ + +#include "Math/RectangularMatrix.h" + +namespace Magnum { namespace Math { namespace Algorithms { + +/** +@brief Gauss-Jordan elimination + +Based on ultra-compact Python code by Jarno Elonen, +http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html. +*/ +class GaussJordan { + public: + /** + * @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) { + for(size_t row = 0; row != size; ++row) { + /* Find max pivot */ + size_t rowMax = row; + for(size_t row2 = row+1; row2 != size; ++row2) + if(std::abs(a(row2, row)) > std::abs(a(rowMax, row))) + rowMax = row2; + + /* Swap the rows */ + std::swap(a[row], a[rowMax]); + std::swap(t[row], t[rowMax]); + + /* Singular */ + if(MathTypeTraits::equals(a(row, row), 0)) + return false; + + /* Eliminate column */ + for(size_t row2 = row+1; row2 != size; ++row2) { + T c = a(row2, row)/a(row, row); + + a[row2] -= a[row]*c; + t[row2] -= t[row]*c; + } + } + + /* Backsubstitute */ + for(size_t row = size; row != 0; --row) { + T c = T(1)/a(row-1, row-1); + + for(size_t row2 = 0; row2 != row-1; ++row2) + t[row2] -= t[row-1]*a(row2, row-1)*c; + + /* Normalize the row */ + t[row-1] *= c; + } + + return true; +} + +}}} + +#endif diff --git a/src/Math/Algorithms/Test/CMakeLists.txt b/src/Math/Algorithms/Test/CMakeLists.txt new file mode 100644 index 000000000..8ebf44bc1 --- /dev/null +++ b/src/Math/Algorithms/Test/CMakeLists.txt @@ -0,0 +1 @@ +corrade_add_test2(MathAlgorithmsGaussJordanTest GaussJordanTest.cpp) diff --git a/src/Math/Algorithms/Test/GaussJordanTest.cpp b/src/Math/Algorithms/Test/GaussJordanTest.cpp new file mode 100644 index 000000000..51e36d89d --- /dev/null +++ b/src/Math/Algorithms/Test/GaussJordanTest.cpp @@ -0,0 +1,60 @@ +/* + Copyright © 2010, 2011, 2012 Vladimír Vondruš + + This file is part of Magnum. + + Magnum is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 3 + only, as published by the Free Software Foundation. + + Magnum is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License version 3 for more details. +*/ + +#include "GaussJordanTest.h" + +#include "Math/Algorithms/GaussJordan.h" + +CORRADE_TEST_MAIN(Magnum::Math::Algorithms::Test::GaussJordanTest) + +namespace Magnum { namespace Math { namespace Algorithms { namespace Test { + +typedef RectangularMatrix<4, 4, float> Matrix4; + +GaussJordanTest::GaussJordanTest() { + addTests(&GaussJordanTest::singular, + &GaussJordanTest::invert); +} + +void GaussJordanTest::singular() { + Matrix4 a(1.0f, 2.0f, 3.0f, 4.0f, + 2.0f, 3.0f, -7.0f, 11.0f, + 2.0f, 4.0f, 6.0f, 8.0f, + 1.0f, 2.0f, 7.0f, 40.0f); + RectangularMatrix<4, 1, float> t; + + CORRADE_VERIFY(!GaussJordan::inPlaceTransposed(a, t)); +} + +void GaussJordanTest::invert() { + Matrix4 a(3.0f, 5.0f, 8.0f, 4.0f, + 4.0f, 4.0f, 7.0f, 3.0f, + 7.0f, -1.0f, 8.0f, 0.0f, + 9.0f, 4.0f, 5.0f, 9.0f); + + Matrix4 expectedInverse(-60/103.0f, 71/103.0f, -4/103.0f, 3/103.0f, + -66/103.0f, 109/103.0f, -25/103.0f, -7/103.0f, + 177/412.0f, -97/206.0f, 53/412.0f, -7/206.0f, + 259/412.0f, -185/206.0f, 31/412.0f, 27/206.0f); + + Matrix4 a2(a); + Matrix4 inverse; + CORRADE_VERIFY(GaussJordan::inPlace(a2, inverse)); + + CORRADE_COMPARE(inverse, expectedInverse); + CORRADE_COMPARE(a*inverse, Matrix4()); +} + +}}}} diff --git a/src/Math/Algorithms/Test/GaussJordanTest.h b/src/Math/Algorithms/Test/GaussJordanTest.h new file mode 100644 index 000000000..0c25d27e0 --- /dev/null +++ b/src/Math/Algorithms/Test/GaussJordanTest.h @@ -0,0 +1,32 @@ +#ifndef Magnum_Math_Algorithms_Test_GaussJordanTest_h +#define Magnum_Math_Algorithms_Test_GaussJordanTest_h +/* + Copyright © 2010, 2011, 2012 Vladimír Vondruš + + This file is part of Magnum. + + Magnum is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 3 + only, as published by the Free Software Foundation. + + Magnum is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License version 3 for more details. +*/ + +#include + +namespace Magnum { namespace Math { namespace Algorithms { namespace Test { + +class GaussJordanTest: public Corrade::TestSuite::Tester { + public: + GaussJordanTest(); + + void singular(); + void invert(); +}; + +}}}} + +#endif diff --git a/src/Math/CMakeLists.txt b/src/Math/CMakeLists.txt index 7053745a8..cf24e5d0b 100644 --- a/src/Math/CMakeLists.txt +++ b/src/Math/CMakeLists.txt @@ -13,6 +13,7 @@ set(MagnumMath_HEADERS Vector4.h) install(FILES ${MagnumMath_HEADERS} DESTINATION ${MAGNUM_INCLUDE_INSTALL_DIR}/Math) +add_subdirectory(Algorithms) add_subdirectory(Geometry) if(BUILD_TESTS)