From 2c61414a838611dc19ce6296164133d642b185fa Mon Sep 17 00:00:00 2001 From: Krishna Vedala <7001608+kvedala@users.noreply.github.com> Date: Thu, 25 Jun 2020 13:38:11 -0400 Subject: [PATCH 1/6] split lu_decomposition to a header file and templated the function --- numerical_methods/lu_decompose.cpp | 74 +++------------------------- numerical_methods/lu_decomposition.h | 65 ++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 66 deletions(-) create mode 100644 numerical_methods/lu_decomposition.h diff --git a/numerical_methods/lu_decompose.cpp b/numerical_methods/lu_decompose.cpp index a0a2d00ab..2ce5fb653 100644 --- a/numerical_methods/lu_decompose.cpp +++ b/numerical_methods/lu_decompose.cpp @@ -7,73 +7,15 @@ #include #include #include -#include -#ifdef _OPENMP -#include -#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> &A, - std::vector> *L, - std::vector> *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 std::ostream &operator<<(std::ostream &out, - std::vector> const &v) { + std::vector> const &v) { const int width = 10; const char separator = ' '; @@ -99,14 +41,14 @@ int main(int argc, char **argv) { std::srand(std::time(NULL)); // random number initializer /* Create a square matrix with random values */ - std::vector> A(mat_size); - std::vector> L(mat_size); // output - std::vector> U(mat_size); // output + std::vector> A(mat_size, + std::valarray(mat_size)); + std::vector> L( + mat_size, std::valarray(mat_size)); // output + std::vector> U( + mat_size, std::valarray(mat_size)); // output for (int i = 0; i < mat_size; i++) { // calloc so that all valeus are '0' by default - A[i] = std::vector(mat_size); - L[i] = std::vector(mat_size); - U[i] = std::vector(mat_size); for (int j = 0; j < mat_size; j++) /* create random values in the limits [-range2, range-1] */ A[i][j] = static_cast(std::rand() % range - range2); diff --git a/numerical_methods/lu_decomposition.h b/numerical_methods/lu_decomposition.h new file mode 100644 index 000000000..9cb881ba9 --- /dev/null +++ b/numerical_methods/lu_decomposition.h @@ -0,0 +1,65 @@ +#pragma once + +#include +#include +#include +#ifdef _OPENMP +#include +#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 + */ +template +int lu_decomposition(const std::vector> &A, + std::vector> *L, + std::vector> *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; +} From e22f56c90628fa02017a43f581384a0c06551821 Mon Sep 17 00:00:00 2001 From: github-actions <${GITHUB_ACTOR}@users.noreply.github.com> Date: Thu, 25 Jun 2020 18:43:00 +0000 Subject: [PATCH 2/6] updating DIRECTORY.md --- DIRECTORY.md | 1 + 1 file changed, 1 insertion(+) diff --git a/DIRECTORY.md b/DIRECTORY.md index 257b4160d..679e6f7d8 100644 --- a/DIRECTORY.md +++ b/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) 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 3/6] 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> &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> &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> &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 +double determinant_lu(const std::vector> &A) { + std::vector> L(A.size(), + std::valarray(A.size())); + std::vector> U(A.size(), + std::valarray(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; +} From c1b0635f993c72fa32de11bc4df73b7dc77509b8 Mon Sep 17 00:00:00 2001 From: Krishna Vedala <7001608+kvedala@users.noreply.github.com> Date: Thu, 25 Jun 2020 15:03:12 -0400 Subject: [PATCH 4/6] create and added matrix type --- numerical_methods/lu_decomposition.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/numerical_methods/lu_decomposition.h b/numerical_methods/lu_decomposition.h index 42c29e4f2..22af249a5 100644 --- a/numerical_methods/lu_decomposition.h +++ b/numerical_methods/lu_decomposition.h @@ -7,6 +7,10 @@ #include #endif +/** Define matrix type as a `std::vector` of `std::valarray` */ +template +using matrix = std::vector>; + /** Perform LU decomposition on matrix * \param[in] A matrix to decompose * \param[out] L output L matrix @@ -15,9 +19,7 @@ * \returns negative if error occurred */ template -int lu_decomposition(const std::vector> &A, - std::vector> *L, - std::vector> *U) { +int lu_decomposition(const matrix &A, matrix *L, matrix *U) { int row, col, j; int mat_size = A.size(); @@ -78,11 +80,9 @@ int lu_decomposition(const std::vector> &A, * @return determinant of matrix A */ template -double determinant_lu(const std::vector> &A) { - std::vector> L(A.size(), - std::valarray(A.size())); - std::vector> U(A.size(), - std::valarray(A.size())); +double determinant_lu(const matrix &A) { + matrix L(A.size(), std::valarray(A.size())); + matrix U(A.size(), std::valarray(A.size())); if (lu_decomposition(A, &L, &U) < 0) return 0; From 84cf0da2bbf3dc7c1dd0dd8ff916317208068cef Mon Sep 17 00:00:00 2001 From: Krishna Vedala <7001608+kvedala@users.noreply.github.com> Date: Thu, 25 Jun 2020 15:03:57 -0400 Subject: [PATCH 5/6] automated self-test of LU decomposition using sample case and determinant checks --- numerical_methods/lu_decompose.cpp | 52 +++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 15 deletions(-) diff --git a/numerical_methods/lu_decompose.cpp b/numerical_methods/lu_decompose.cpp index 2ce5fb653..b27fed2ee 100644 --- a/numerical_methods/lu_decompose.cpp +++ b/numerical_methods/lu_decompose.cpp @@ -4,6 +4,7 @@ * square matrix * \author [Krishna Vedala](https://github.com/kvedala) */ +#include #include #include #include @@ -14,8 +15,7 @@ * operator to print a matrix */ template -std::ostream &operator<<(std::ostream &out, - std::vector> const &v) { +std::ostream &operator<<(std::ostream &out, matrix const &v) { const int width = 10; const char separator = ' '; @@ -29,24 +29,19 @@ 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> A(mat_size, - std::valarray(mat_size)); - std::vector> L( - mat_size, std::valarray(mat_size)); // output - std::vector> U( - mat_size, std::valarray(mat_size)); // output + matrix A(mat_size, std::valarray(mat_size)); + matrix L(mat_size, std::valarray(mat_size)); // output + matrix U(mat_size, std::valarray(mat_size)); // output for (int i = 0; i < mat_size; i++) { // calloc so that all valeus are '0' by default for (int j = 0; j < mat_size; j++) @@ -63,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 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 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 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; } From 68dd9b1235f1d2a33d0a873a33ff55846bc8a4f6 Mon Sep 17 00:00:00 2001 From: Krishna Vedala <7001608+kvedala@users.noreply.github.com> Date: Thu, 25 Jun 2020 15:22:02 -0400 Subject: [PATCH 6/6] added file documentation --- numerical_methods/lu_decomposition.h | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/numerical_methods/lu_decomposition.h b/numerical_methods/lu_decomposition.h index 22af249a5..4999fc40f 100644 --- a/numerical_methods/lu_decomposition.h +++ b/numerical_methods/lu_decomposition.h @@ -1,3 +1,10 @@ +/** + * @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