matrix: adding U-inverse function

src/lol/math/matrix.h

@@ -570,7 +570,7 @@ mat_t<T, N, N> l_inverse(mat_t<T, N, N> const & L)
         {
             T sum = 0;
 
-            for (int k = i ; k >= j ; --k)
+            for (int k = i ; k > j ; --k)
                 sum += ret[k][i] * L[j][k];
 
             ret[j][i] = ((i == j ? 1 : 0) - sum) / L[j][j];
@@ -580,6 +580,31 @@ mat_t<T, N, N> l_inverse(mat_t<T, N, N> const & L)
     return ret;
 }
 
+/*
+ * Compute U matrix inverse
+ */
+
+template<typename T, int N>
+mat_t<T, N, N> u_inverse(mat_t<T, N, N> const & U)
+{
+    mat_t<T, N, N> ret;
+
+    for (int i = 0 ; i < N ; ++i)
+    {
+        for (int j = i ; j < N ; ++j)
+        {
+            T sum = 0;
+
+            for (int k = i ; k < j ; ++k)
+                sum += ret[k][i] * U[j][k];
+
+            ret[j][i] = ((i == j ? 1 : 0) - sum) / U[j][j];
+        }
+    }
+
+    return ret;
+}
+
 /*
  * Compute square matrix inverse
  */

src/t/math/matrix.cpp

@@ -187,6 +187,34 @@ lolunit_declare_fixture(MatrixTest)
                 lolunit_assert_doubles_equal(identity[i][j], i == j ? 1 : 0, 1e-5);
     }
 
+    lolunit_declare_test(UInverse3x3)
+    {
+        mat3 m0 = inv3;
+        mat3 L, U;
+        lu_decomposition(inv3, L, U);
+        mat3 u_inv = u_inverse(U);
+
+        mat3 identity = u_inv * U;
+
+        for (int i = 0 ; i < 3 ; ++i)
+            for (int j = 0 ; j < 3 ; ++j)
+                lolunit_assert_doubles_equal(identity[i][j], i == j ? 1 : 0, 1e-5);
+    }
+
+    lolunit_declare_test(UInverse4x4)
+    {
+        mat4 m0 = inv4;
+        mat4 L, U;
+        lu_decomposition(inv4, L, U);
+        mat4 u_inv = u_inverse(U);
+
+        mat4 identity = u_inv * U;
+
+        for (int i = 0 ; i < 4 ; ++i)
+            for (int j = 0 ; j < 4 ; ++j)
+                lolunit_assert_doubles_equal(identity[i][j], i == j ? 1 : 0, 1e-5);
+    }
+
     lolunit_declare_test(Inverse3x3)
     {
         mat3 m0 = inv3;