From 45d067d6f6c4418d38d36d93ad544e7f5ffe8a90 Mon Sep 17 00:00:00 2001 From: jdiestra Date: Mon, 16 Oct 2023 20:20:56 +0200 Subject: [PATCH 1/3] Updating forward substitution documentation --- documentation/main_page.md | 60 +++++++++++++++++++++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/documentation/main_page.md b/documentation/main_page.md index 9d57af6..11c39b4 100644 --- a/documentation/main_page.md +++ b/documentation/main_page.md @@ -85,7 +85,65 @@ A **Vector** object can be created as, ### Forward Substitution -TBD +Reference : Matrix Algorithms, G.W. Stewart. Page 88 + +For the following equation, + +\f$ Lx = b \f$ + +where, + +\f$ L \f$ : Nonsingular lower triangular matrix of order \f$ n \f$. + +\f$ b \f$ : Vector or \f$ n \f$ order. + +Then \f$ x \f$ can be calculated by, + + for k = 1 to n + x[k] = b[k] + for j = 1 to k-1 + x[k] = x[k] - L[k,j]x[j] + end for j + x[k] = x[k]/L[k,k] + end for k + +This algorihm is implemented by the function, + + marsvin::Vector forward_substitution(const ::marsvin::Matrix& L, const ::marsvin::Vector& b) + +For example, + + // L : Nonsingular lower triangular matrix of nxn size. + // b : Vector of n size. + marsvin::Vector x; + x = marsvin::forward_substitution(L, b); + +or, + + // L : Nonsingular lower triangular matrix of nxn size. + // b : Vector of n size. + marsvin::Vector x; + marsvin::forward_substitution(L, b, x); + +In order to avoid using memory to allocate \f$ x \f$, +we can overwrite the solution in vector \f$ b \f$, + + for k = 1 to n + for j = 1 to k-1 + b[k] = b[k] - L[k,j]b[j] + end for j + b[k] = b[k]/L[k,k] + end for k + +This algorihm is implemented by the function, + + marsvin::Vector forward_substitution_memory(const ::marsvin::Matrix& L, ::marsvin::Vector& b) + +For example, + + // L : Nonsingular lower triangular matrix of nxn size. + // b : Vector of n size. + marsvin::forward_substitution_memory(L, b); ### Backward Substitution From a2812724e6aacf7be7edb73b8769c9b1a3c3817e Mon Sep 17 00:00:00 2001 From: jdiestra Date: Tue, 26 Dec 2023 11:52:58 +0100 Subject: [PATCH 2/3] Adding backward substitution --- .../marsvin/error_handling/status_code.hpp | 1 + .../marsvin/error_handling/status_type.hpp | 5 +- .../backward_substitution.hpp | 49 +++++++++++++ .../forward_substitution.hpp | 1 + core/src/marsvin_exception.cpp | 8 +++ core/src/marsvin_status_code.cpp | 5 ++ documentation/main_page.md | 68 +++++++++++++++++-- test/CMakeLists.txt | 7 ++ test/test_marsvin_backward_substitution.cpp | 45 ++++++++++++ test/test_marsvin_base_matrix.cpp | 15 ++++ test/test_marsvin_status_code.cpp | 1 + 11 files changed, 196 insertions(+), 9 deletions(-) create mode 100644 core/include/marsvin/triangular_systems/backward_substitution.hpp create mode 100644 test/test_marsvin_backward_substitution.cpp diff --git a/core/include/marsvin/error_handling/status_code.hpp b/core/include/marsvin/error_handling/status_code.hpp index 2ba1f84..f4c27d8 100644 --- a/core/include/marsvin/error_handling/status_code.hpp +++ b/core/include/marsvin/error_handling/status_code.hpp @@ -36,6 +36,7 @@ class StatusCode { static StatusType TypeErrorSquareMatrix(); static StatusType TypeErrorEqualSize(); static StatusType TypeErrorSolveLinearEquantionDimensions(); + static StatusType TypeErrorMatrixDiagonalHasZero(); private: StatusType status_type_; diff --git a/core/include/marsvin/error_handling/status_type.hpp b/core/include/marsvin/error_handling/status_type.hpp index 8c7a2ed..439736f 100644 --- a/core/include/marsvin/error_handling/status_type.hpp +++ b/core/include/marsvin/error_handling/status_type.hpp @@ -31,8 +31,9 @@ enum class StatusType { kErrorSquareMatrix, // Error while a function is not applied to an square // matrix. kErrorEqualSize, // Error while same size of matrix/vector is needed. - kErrorSolveLinearEquantionDimensions // Error while A*x=b, A.columns() or - // A.rows() != b.size() + kErrorSolveLinearEquantionDimensions, // Error while A*x=b, A.columns() or + // A.rows() != b.size() + kErrorMatrixDiagonalHasZero // Error while diagonal of matrix has a zero. }; } // namespace marsvin diff --git a/core/include/marsvin/triangular_systems/backward_substitution.hpp b/core/include/marsvin/triangular_systems/backward_substitution.hpp new file mode 100644 index 0000000..9d98c93 --- /dev/null +++ b/core/include/marsvin/triangular_systems/backward_substitution.hpp @@ -0,0 +1,49 @@ +/** + * @file forward_substitution.hpp + * + */ + +#ifndef MARSVIN_TRIANGULAR_SYSTEMS_BACKWARD_SUBSTITUTION_HPP_ +#define MARSVIN_TRIANGULAR_SYSTEMS_BACKWARD_SUBSTITUTION_HPP_ + +#include "marsvin/containers/matrix.hpp" +#include "marsvin/containers/vector.hpp" + +namespace marsvin { + +template +void backward_substitution(const ::marsvin::Matrix& U, + const ::marsvin::Vector& b, + ::marsvin::Vector& x) { + // Check U is square + if (!U.is_square()) { + throw marsvin::Exception(marsvin::StatusCode::TypeErrorSquareMatrix()); + } + // Check if U rows or columns are equal to size of b + if (U.rows() != b.size()) { + throw marsvin::Exception( + marsvin::StatusCode::TypeErrorSolveLinearEquantionDimensions()); + } + // Check x and b sizes + if (x.empty()) { + x.resize(b.size()); + } else { + if (x.size() != b.size()) { + throw marsvin::Exception(marsvin::StatusCode::TypeErrorEqualSize()); + } + } + for (int k = U.rows() - 1; k >= 0; --k) { + // TBD : Add check if U.at(k,k) is close to Zero. + x.at(k) = b.at(k); + if (k >= 0) { + for (std::size_t j = k + 1; j <= U.rows() -1 ; ++j) { + x.at(k) = x.at(k) - U.at(k, j) * x.at(j); + } + } + x.at(k) = x.at(k) / U.at(k, k); + } +}; + +} // namespace marsvin + +#endif // MARSVIN_TRIANGULAR_SYSTEMS_BACKWARD_SUBSTITUTION_HPP_ diff --git a/core/include/marsvin/triangular_systems/forward_substitution.hpp b/core/include/marsvin/triangular_systems/forward_substitution.hpp index 0b62071..cd168b1 100644 --- a/core/include/marsvin/triangular_systems/forward_substitution.hpp +++ b/core/include/marsvin/triangular_systems/forward_substitution.hpp @@ -34,6 +34,7 @@ void forward_substitution(const ::marsvin::Matrix& L, } } for (std::size_t k = 0; k < L.rows(); ++k) { + // TBD : Add check if L.at(k,k) is close to Zero. x.at(k) = b.at(k); if (k > 0) { for (std::size_t j = 0; j <= (k - 1); ++j) { diff --git a/core/src/marsvin_exception.cpp b/core/src/marsvin_exception.cpp index 490d748..afafe34 100644 --- a/core/src/marsvin_exception.cpp +++ b/core/src/marsvin_exception.cpp @@ -56,10 +56,18 @@ std::string Exception::StatusCodeToString(const StatusCode& status_code) { break; case marsvin::StatusType::kErrorSquareMatrix: string_ = "Operation should be applied in an square matrix."; + break; case marsvin::StatusType::kErrorEqualSize: string_ = "Equal size of matrix or vector is required in operation."; + break; + case marsvin::StatusType::kErrorSolveLinearEquantionDimensions: + string_ = + "Error solving linear equation A*x=b. " + "Matrix number of rows or columns must be equal to vector size."; + break; default: + string_ = "Unknown error."; break; } return string_; diff --git a/core/src/marsvin_status_code.cpp b/core/src/marsvin_status_code.cpp index e81a5bb..e08543c 100644 --- a/core/src/marsvin_status_code.cpp +++ b/core/src/marsvin_status_code.cpp @@ -77,4 +77,9 @@ StatusType StatusCode::TypeErrorEqualSize() { StatusType StatusCode::TypeErrorSolveLinearEquantionDimensions() { return StatusType::kErrorSolveLinearEquantionDimensions; } + +StatusType StatusCode::TypeErrorMatrixDiagonalHasZero() { + return StatusType::kErrorMatrixDiagonalHasZero; +} + } // namespace marsvin diff --git a/documentation/main_page.md b/documentation/main_page.md index 11c39b4..3324db5 100644 --- a/documentation/main_page.md +++ b/documentation/main_page.md @@ -85,8 +85,6 @@ A **Vector** object can be created as, ### Forward Substitution -Reference : Matrix Algorithms, G.W. Stewart. Page 88 - For the following equation, \f$ Lx = b \f$ @@ -99,9 +97,9 @@ where, Then \f$ x \f$ can be calculated by, - for k = 1 to n + for k = 0 to n-1 x[k] = b[k] - for j = 1 to k-1 + for j = 0 to k-1 x[k] = x[k] - L[k,j]x[j] end for j x[k] = x[k]/L[k,k] @@ -128,8 +126,8 @@ or, In order to avoid using memory to allocate \f$ x \f$, we can overwrite the solution in vector \f$ b \f$, - for k = 1 to n - for j = 1 to k-1 + for k = 0 to n-1 + for j = 0 to k-1 b[k] = b[k] - L[k,j]b[j] end for j b[k] = b[k]/L[k,k] @@ -147,7 +145,63 @@ For example, ### Backward Substitution -TBD +For the following equation, + +\f$ Ux = b \f$ + +where, + +\f$ U \f$ : Nonsingular upper triangular matrix of order \f$ n \f$. + +\f$ b \f$ : Vector or \f$ n \f$ order. + +Then \f$ x \f$ can be calculated by, + + for k = n-1 to 0 + x[k] = b[k] + for j = k+1 to n-1 + x[k] = x[k] - U[k,j]x[j] + end for j + x[k] = x[k]/U[k,k] + end for k + +This algorihm is implemented by the function, + + marsvin::Vector backward_substitution(const ::marsvin::Matrix& L, const ::marsvin::Vector& b) + +For example, + + // U : Nonsingular upper triangular matrix of nxn size. + // b : Vector of n size. + marsvin::Vector x; + x = marsvin::backward_substitution(L, b); + +or, + + // U : Nonsingular upper triangular matrix of nxn size. + // b : Vector of n size. + marsvin::Vector x; + marsvin::backward_substitution(L, b, x); + +In order to avoid using memory to allocate \f$ x \f$, +we can overwrite the solution in vector \f$ b \f$, + + for k = n-1 to 0 + for j = k+1 to n-1 + b[k] = b[k] - U[k,j]b[j] + end for j + b[k] = b[k]/U[k,k] + end for k + +This algorihm is implemented by the function, + + marsvin::Vector backward_substitution_memory(const ::marsvin::Matrix& U, ::marsvin::Vector& b) + +For example, + + // U : Nonsingular upper triangular matrix of nxn size. + // b : Vector of n size. + marsvin::backward_substitution_memory(L, b); ### Gauss Elimination diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 47127fe..7a3616b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -65,3 +65,10 @@ MarsvinAddTest( INCLUDE_DIR ${PROJECT_SOURCE_DIR}/core/include LINK_LIBRARIES marsvin ) + +MarsvinAddTest( + NAME test_marsvin_backward_substitution + SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/test_marsvin_backward_substitution.cpp + INCLUDE_DIR ${PROJECT_SOURCE_DIR}/core/include + LINK_LIBRARIES marsvin +) diff --git a/test/test_marsvin_backward_substitution.cpp b/test/test_marsvin_backward_substitution.cpp new file mode 100644 index 0000000..19caeff --- /dev/null +++ b/test/test_marsvin_backward_substitution.cpp @@ -0,0 +1,45 @@ + +#include "gtest/gtest.h" +#include "marsvin/triangular_systems/backward_substitution.hpp" +#include "marsvin/tools/logger.hpp" +#include "marsvin/tools/compare.hpp" + +/* + * Description: + * Check algorithm by an example. + * Operation: U*x = b + * Choose value for U and b, and calculate "x". + * Then check if U*x == b + * x_expected = {1,2,3,4} + */ +TEST(BackwardSubstitution, Test_x_empty_4x4) { + marsvin::Logger logger_; + marsvin::Matrix U = {{3,2,1,1}, + {0,1,2,1}, + {0,0,1,1}, + {0,0,0,1}}; + marsvin::Vector b = {14,12,7,4}; + marsvin::Vector x_real = {1,2,3,4}; + marsvin::Vector x; + marsvin::backward_substitution(U, b, x); + EXPECT_EQ(x.size(), b.size()); + EXPECT_TRUE(marsvin::tools::compare(x,x_real)); + EXPECT_TRUE(marsvin::tools::compare(U*x,b)); +} + +TEST(BackwardSubstitution, Test_x_empty_5x5) { + marsvin::Logger logger_; + marsvin::Matrix U = {{1,2,3,4,5}, + {0,6,7,8,9}, + {0,0,10,11,12}, + {0,0,0,13,14}, + {0,0,0,0,15}}; + marsvin::Vector b = {55,110,134,122,75}; + marsvin::Vector x_real = {1,2,3,4,5}; + marsvin::Vector x; + marsvin::backward_substitution(U, b, x); + EXPECT_EQ(x.size(), b.size()); + EXPECT_TRUE(marsvin::tools::compare(x,x_real)); + EXPECT_TRUE(marsvin::tools::compare(U*x,b)); +} + diff --git a/test/test_marsvin_base_matrix.cpp b/test/test_marsvin_base_matrix.cpp index 695d4da..8e83a35 100644 --- a/test/test_marsvin_base_matrix.cpp +++ b/test/test_marsvin_base_matrix.cpp @@ -400,3 +400,18 @@ TEST(BaseMatrix, method_base_matrix_multiplication) { EXPECT_TRUE(marsvin::tools::compare(m_result,m_result_expected,tolerance)); } +TEST(BaseMatrix, multiplication_matrix_vector) { + constexpr std::size_t lhs_kRows = 5; + constexpr std::size_t lhs_kColumns = 5; + constexpr std::size_t rhs_kRows = 5; + constexpr std::size_t rhs_kColumns = 1; + marsvin::BaseMatrix m_lhs(lhs_kRows, lhs_kColumns, {1,2,3,4,5,0,6,7,8,9,0,0,10,11,12,0,0,0,13,14,0,0,0,0,15}); + marsvin::BaseMatrix m_rhs(rhs_kRows, rhs_kColumns, {1,2,3,4,5}); + marsvin::BaseMatrix m_result_expected(lhs_kRows, rhs_kColumns, {55,110,134,122,75}); + auto m_result = m_lhs*m_rhs; + EXPECT_EQ(m_result.rows(), lhs_kRows); + EXPECT_EQ(m_result.columns(), rhs_kColumns); + int tolerance = 0; + EXPECT_TRUE(marsvin::tools::compare(m_result, m_result_expected,tolerance)); +} + diff --git a/test/test_marsvin_status_code.cpp b/test/test_marsvin_status_code.cpp index 2f6688e..a5474c6 100644 --- a/test/test_marsvin_status_code.cpp +++ b/test/test_marsvin_status_code.cpp @@ -15,6 +15,7 @@ TEST(StatusCode, StaticMethods) { ASSERT_EQ(marsvin::StatusType::kErrorSquareMatrix, marsvin::StatusCode::TypeErrorSquareMatrix()); ASSERT_EQ(marsvin::StatusType::kErrorEqualSize, marsvin::StatusCode::TypeErrorEqualSize()); ASSERT_EQ(marsvin::StatusType::kErrorSolveLinearEquantionDimensions, marsvin::StatusCode::TypeErrorSolveLinearEquantionDimensions()); + ASSERT_EQ(marsvin::StatusType::kErrorMatrixDiagonalHasZero, marsvin::StatusCode::TypeErrorMatrixDiagonalHasZero()); } TEST(StatusCode, Constructor_01) { From 269dca4c062d2e7acbe09a5968fcf288fa81faa1 Mon Sep 17 00:00:00 2001 From: jdiestra Date: Tue, 26 Dec 2023 11:59:02 +0100 Subject: [PATCH 3/3] Fixing clang-format --- .../marsvin/triangular_systems/backward_substitution.hpp | 6 +++--- core/src/marsvin_exception.cpp | 3 ++- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/core/include/marsvin/triangular_systems/backward_substitution.hpp b/core/include/marsvin/triangular_systems/backward_substitution.hpp index 9d98c93..a1dd864 100644 --- a/core/include/marsvin/triangular_systems/backward_substitution.hpp +++ b/core/include/marsvin/triangular_systems/backward_substitution.hpp @@ -13,8 +13,8 @@ namespace marsvin { template void backward_substitution(const ::marsvin::Matrix& U, - const ::marsvin::Vector& b, - ::marsvin::Vector& x) { + const ::marsvin::Vector& b, + ::marsvin::Vector& x) { // Check U is square if (!U.is_square()) { throw marsvin::Exception(marsvin::StatusCode::TypeErrorSquareMatrix()); @@ -36,7 +36,7 @@ void backward_substitution(const ::marsvin::Matrix& U, // TBD : Add check if U.at(k,k) is close to Zero. x.at(k) = b.at(k); if (k >= 0) { - for (std::size_t j = k + 1; j <= U.rows() -1 ; ++j) { + for (std::size_t j = k + 1; j <= U.rows() - 1; ++j) { x.at(k) = x.at(k) - U.at(k, j) * x.at(j); } } diff --git a/core/src/marsvin_exception.cpp b/core/src/marsvin_exception.cpp index afafe34..3e7cfbe 100644 --- a/core/src/marsvin_exception.cpp +++ b/core/src/marsvin_exception.cpp @@ -64,7 +64,8 @@ std::string Exception::StatusCodeToString(const StatusCode& status_code) { case marsvin::StatusType::kErrorSolveLinearEquantionDimensions: string_ = "Error solving linear equation A*x=b. " - "Matrix number of rows or columns must be equal to vector size."; + "Matrix number of rows or columns must be equal to vector " + "size."; break; default: string_ = "Unknown error.";