Gauss-Jordan elimination.

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
  */

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()

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š <mosra@centrum.cz>
+
+    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<size_t size, size_t rows, class T> static bool inPlaceTransposed(RectangularMatrix<size, size, T>& a, RectangularMatrix<size, rows, T>& t);
+
+        /**
+         * @brief Eliminate in place
+         *
+         * Transposes the matrices, calls inPlaceTransposed() on them and then
+         * transposes them back.
+         */
+        template<size_t size, size_t cols, class T> static bool inPlace(RectangularMatrix<size, size, T>& a, RectangularMatrix<cols, size, T>& t) {
+            a = a.transposed();
+            RectangularMatrix<size, cols, T> tTransposed = t.transposed();
+
+            bool ret = inPlaceTransposed(a, tTransposed);
+
+            a = a.transposed();
+            t = tTransposed.transposed();
+
+            return ret;
+        }
+};
+
+template<size_t size, size_t cols, class T> bool GaussJordan::inPlaceTransposed(RectangularMatrix<size, size, T>& a, RectangularMatrix<size, cols, T>& 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<T>::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

src/Math/Algorithms/Test/CMakeLists.txt

@@ -0,0 +1 @@
+corrade_add_test2(MathAlgorithmsGaussJordanTest GaussJordanTest.cpp)

src/Math/Algorithms/Test/GaussJordanTest.cpp

@@ -0,0 +1,60 @@
+/*
+    Copyright © 2010, 2011, 2012 Vladimír Vondruš <mosra@centrum.cz>
+
+    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());
+}
+
+}}}}

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š <mosra@centrum.cz>
+
+    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 <TestSuite/Tester.h>
+
+namespace Magnum { namespace Math { namespace Algorithms { namespace Test {
+
+class GaussJordanTest: public Corrade::TestSuite::Tester<GaussJordanTest> {
+    public:
+        GaussJordanTest();
+
+        void singular();
+        void invert();
+};
+
+}}}}
+
+#endif

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)