Merge pull request #906 from kvedala/lu_decompose

enhancement: LU decompose to header file and self-tests

DIRECTORY.md

@@ -140,6 +140,7 @@
   * [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/gaussian_elimination.cpp)
   * [Golden Search Extrema](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/golden_search_extrema.cpp)
   * [Lu Decompose](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/lu_decompose.cpp)
+  * [Lu Decomposition](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/lu_decomposition.h)
   * [Newton Raphson Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/newton_raphson_method.cpp)
   * [Ode Forward Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_forward_euler.cpp)
   * [Ode Midpoint Euler](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ode_midpoint_euler.cpp)

numerical_methods/lu_decompose.cpp

@@ -4,76 +4,18 @@
  * square matrix
  * \author [Krishna Vedala](https://github.com/kvedala)
  */
+#include <cassert>
 #include <ctime>
 #include <iomanip>
 #include <iostream>
-#include <vector>
-#ifdef _OPENMP
-#include <omp.h>
-#endif
 
-/** Perform LU decomposition on matrix
- * \param[in] A matrix to decompose
- * \param[out] L output L matrix
- * \param[out] U output U matrix
- * \returns 0 if no errors
- * \returns negative if error occurred
- */
-int lu_decomposition(const std::vector<std::vector<double>> &A,
-                     std::vector<std::vector<double>> *L,
-                     std::vector<std::vector<double>> *U) {
-    int row, col, j;
-    int mat_size = A.size();
-
-    if (mat_size != A[0].size()) {
-        // check matrix is a square matrix
-        std::cerr << "Not a square matrix!\n";
-        return -1;
-    }
-
-    // regularize each row
-    for (row = 0; row < mat_size; row++) {
-        // Upper triangular matrix
-#ifdef _OPENMP
-#pragma omp for
-#endif
-        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];
-
-            // Evaluate U[i,k]
-            U[0][row][col] = A[row][col] - lu_sum;
-        }
-
-        // Lower triangular matrix
-#ifdef _OPENMP
-#pragma omp for
-#endif
-        for (col = row; col < mat_size; col++) {
-            if (row == col) {
-                L[0][row][col] = 1.;
-                continue;
-            }
-
-            // 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];
-
-            // Evaluate U[i,k]
-            L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
-        }
-    }
-
-    return 0;
-}
+#include "./lu_decomposition.h"
 
 /**
  * operator to print a matrix
  */
 template <typename T>
-std::ostream &operator<<(std::ostream &out,
-                         std::vector<std::vector<T>> const &v) {
+std::ostream &operator<<(std::ostream &out, matrix<T> const &v) {
     const int width = 10;
     const char separator = ' ';
 
@@ -87,26 +29,21 @@ std::ostream &operator<<(std::ostream &out,
     return out;
 }
 
-/** Main function */
-int main(int argc, char **argv) {
+/**
+ * Test LU decomposition
+ * \todo better ways to self-check a matrix output?
+ */
+void test1() {
     int mat_size = 3;  // default matrix size
     const int range = 50;
     const int range2 = range >> 1;
 
-    if (argc == 2)
-        mat_size = atoi(argv[1]);
-
-    std::srand(std::time(NULL));  // random number initializer
-
     /* Create a square matrix with random values */
-    std::vector<std::vector<double>> A(mat_size);
-    std::vector<std::vector<double>> L(mat_size);  // output
-    std::vector<std::vector<double>> U(mat_size);  // output
+    matrix<double> A(mat_size, std::valarray<double>(mat_size));
+    matrix<double> L(mat_size, std::valarray<double>(mat_size));  // output
+    matrix<double> U(mat_size, std::valarray<double>(mat_size));  // output
     for (int i = 0; i < mat_size; i++) {
         // calloc so that all valeus are '0' by default
-        A[i] = std::vector<double>(mat_size);
-        L[i] = std::vector<double>(mat_size);
-        U[i] = std::vector<double>(mat_size);
         for (int j = 0; j < mat_size; j++)
             /* create random values in the limits [-range2, range-1] */
             A[i][j] = static_cast<double>(std::rand() % range - range2);
@@ -121,6 +58,33 @@ int main(int argc, char **argv) {
     std::cout << "A = \n" << A << "\n";
     std::cout << "L = \n" << L << "\n";
     std::cout << "U = \n" << U << "\n";
+}
 
+/**
+ * @brief Test determinant computation using LU decomposition
+ */
+void test2() {
+    std::cout << "Determinant test 1...";
+    matrix<int> A1({{1, 2, 3}, {4, 9, 6}, {7, 8, 9}});
+    assert(determinant_lu(A1) == -48);
+    std::cout << "passed\n";
+
+    std::cout << "Determinant test 2...";
+    matrix<int> A2({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}});
+    assert(determinant_lu(A2) == 0);
+    std::cout << "passed\n";
+
+    std::cout << "Determinant test 3...";
+    matrix<float> A3({{1.2, 2.3, 3.4}, {4.5, 5.6, 6.7}, {7.8, 8.9, 9.0}});
+    assert(determinant_lu(A3) == 3.63);
+    std::cout << "passed\n";
+}
+
+/** Main function */
+int main(int argc, char **argv) {
+    std::srand(std::time(NULL));  // random number initializer
+
+    test1();
+    test2();
     return 0;
 }

numerical_methods/lu_decomposition.h

@@ -0,0 +1,102 @@
+/**
+ * @file lu_decomposition.h
+ * @author [Krishna Vedala](https://github.com/kvedala)
+ * @brief Functions associated with [LU
+ * Decomposition](https://en.wikipedia.org/wiki/LU_decomposition)
+ * of a square matrix.
+ */
+#pragma once
+
+#include <iostream>
+#include <valarray>
+#include <vector>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+/** Define matrix type as a `std::vector` of `std::valarray` */
+template <typename T>
+using matrix = std::vector<std::valarray<T>>;
+
+/** Perform LU decomposition on matrix
+ * \param[in] A matrix to decompose
+ * \param[out] L output L matrix
+ * \param[out] U output U matrix
+ * \returns 0 if no errors
+ * \returns negative if error occurred
+ */
+template <typename T>
+int lu_decomposition(const matrix<T> &A, matrix<double> *L, matrix<double> *U) {
+    int row, col, j;
+    int mat_size = A.size();
+
+    if (mat_size != A[0].size()) {
+        // check matrix is a square matrix
+        std::cerr << "Not a square matrix!\n";
+        return -1;
+    }
+
+    // regularize each row
+    for (row = 0; row < mat_size; row++) {
+        // Upper triangular matrix
+#ifdef _OPENMP
+#pragma omp for
+#endif
+        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];
+            }
+
+            // Evaluate U[i,k]
+            U[0][row][col] = A[row][col] - lu_sum;
+        }
+
+        // Lower triangular matrix
+#ifdef _OPENMP
+#pragma omp for
+#endif
+        for (col = row; col < mat_size; col++) {
+            if (row == col) {
+                L[0][row][col] = 1.;
+                continue;
+            }
+
+            // 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];
+            }
+
+            // Evaluate U[i,k]
+            L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
+        }
+    }
+
+    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 matrix<T> &A) {
+    matrix<double> L(A.size(), std::valarray<double>(A.size()));
+    matrix<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;
+}