From d565a73a1543263d28557bfa0797700f26d39e1f Mon Sep 17 00:00:00 2001
From: Guillaume Bittoun <guillaume.bittoun@gmail.com>
Date: Mon, 19 Jan 2015 21:16:58 +0000
Subject: [PATCH] matrix: starting LU-decomposition

---
 src/lol/math/matrix.h | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/src/lol/math/matrix.h b/src/lol/math/matrix.h
index da316838..4df85e4f 100644
--- a/src/lol/math/matrix.h
+++ b/src/lol/math/matrix.h
@@ -528,6 +528,35 @@ T const & determinant(mat_t<T, 1, 1> const &m)
     return m[0][0];
 }
 
+/*
+ * Compute LU-decomposition
+ */
+
+template<typename T, int N>
+void lu_decomposition(mat_t<T, N, N> const &m, mat_t<T, N, N> & L, mat_t<T, N, N> & U)
+{
+    for (int i = 0; i < N; ++i)
+    {
+        for (int j = 0; j < N; ++j)
+        {
+            T sum = 0;
+            for (int k = 0 ; k < min(i, j) ; ++k)
+                sum += L[i][k] * U[k][j];
+
+            if (i < j)
+            {
+                L[i][j] = 0;
+                U[i][j] = (m[i][j] - sum) / L[i][i];
+            }
+            else /* if (i >= j) */
+            {
+                U[i][j] = i == j ? 1 : 0;
+                L[i][j] = (m[i][j] - sum) / U[j][j];
+            }
+        }
+    }
+}
+
 /*
  * Compute square matrix inverse
  */