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) diff --git a/numerical_methods/lu_decompose.cpp b/numerical_methods/lu_decompose.cpp index a0a2d00ab..b27fed2ee 100644 --- a/numerical_methods/lu_decompose.cpp +++ b/numerical_methods/lu_decompose.cpp @@ -4,76 +4,18 @@ * square matrix * \author [Krishna Vedala](https://github.com/kvedala) */ +#include #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::ostream &operator<<(std::ostream &out, matrix 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> A(mat_size); - std::vector> L(mat_size); // output - std::vector> U(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 - 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); @@ -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 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; } diff --git a/numerical_methods/lu_decomposition.h b/numerical_methods/lu_decomposition.h new file mode 100644 index 000000000..4999fc40f --- /dev/null +++ b/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 +#include +#include +#ifdef _OPENMP +#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 + * \param[out] U output U matrix + * \returns 0 if no errors + * \returns negative if error occurred + */ +template +int lu_decomposition(const matrix &A, matrix *L, matrix *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 +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; + + double result = 1.f; + for (size_t i = 0; i < A.size(); i++) { + result *= L[i][i] * U[i][i]; + } + return result; +}