From f29c14032a73160444f5be4df03a2c03ca4acbb5 Mon Sep 17 00:00:00 2001
From: Krishna Vedala <7001608+kvedala@users.noreply.github.com>
Date: Thu, 25 Jun 2020 14:51:54 -0400
Subject: [PATCH] added determinant computation using LU decomposition

---
 numerical_methods/lu_decomposition.h | 34 ++++++++++++++++++++++++++--
 1 file changed, 32 insertions(+), 2 deletions(-)

diff --git a/numerical_methods/lu_decomposition.h b/numerical_methods/lu_decomposition.h
index 9cb881ba9..42c29e4f2 100644
--- a/numerical_methods/lu_decomposition.h
+++ b/numerical_methods/lu_decomposition.h
@@ -36,7 +36,9 @@ int lu_decomposition(const std::vector<std::valarray<T>> &A,
         for (col = row; col < mat_size; col++) {
             // Summation of L[i,j] * U[j,k]
             double lu_sum = 0.;
-            for (j = 0; j < row; j++) lu_sum += L[0][row][j] * U[0][j][col];
+            for (j = 0; j < row; j++) {
+                lu_sum += L[0][row][j] * U[0][j][col];
+            }
 
             // Evaluate U[i,k]
             U[0][row][col] = A[row][col] - lu_sum;
@@ -54,7 +56,9 @@ int lu_decomposition(const std::vector<std::valarray<T>> &A,
 
             // Summation of L[i,j] * U[j,k]
             double lu_sum = 0.;
-            for (j = 0; j < row; j++) lu_sum += L[0][col][j] * U[0][j][row];
+            for (j = 0; j < row; j++) {
+                lu_sum += L[0][col][j] * U[0][j][row];
+            }
 
             // Evaluate U[i,k]
             L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
@@ -63,3 +67,29 @@ int lu_decomposition(const std::vector<std::valarray<T>> &A,
 
     return 0;
 }
+
+/**
+ * @brief Compute determinant of an NxN square matrix using LU decomposition.
+ * Using LU decomposition, the determinant is given by the product of diagonal
+ * elements of matrices L and U.
+ *
+ * @tparam T datatype of input matrix - int, unsigned int, double, etc
+ * @param A input square matrix
+ * @return determinant of matrix A
+ */
+template <typename T>
+double determinant_lu(const std::vector<std::valarray<T>> &A) {
+    std::vector<std::valarray<double>> L(A.size(),
+                                         std::valarray<double>(A.size()));
+    std::vector<std::valarray<double>> U(A.size(),
+                                         std::valarray<double>(A.size()));
+
+    if (lu_decomposition(A, &L, &U) < 0)
+        return 0;
+
+    double result = 1.f;
+    for (size_t i = 0; i < A.size(); i++) {
+        result *= L[i][i] * U[i][i];
+    }
+    return result;
+}