Update forward and backward substitution

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_;

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

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<typename T>
+void backward_substitution(const ::marsvin::Matrix<T>& U,
+                           const ::marsvin::Vector<T>& b,
+                           ::marsvin::Vector<T>& 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_

core/include/marsvin/triangular_systems/forward_substitution.hpp

@@ -34,6 +34,7 @@ void forward_substitution(const ::marsvin::Matrix<T>& 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) {

core/src/marsvin_exception.cpp

@@ -56,10 +56,19 @@ 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_;

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

documentation/main_page.md

@@ -85,11 +85,123 @@ A **Vector** object can be created as,
 
 ### Forward Substitution
 
-TBD
+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 = 0 to n-1
+        x[k] = b[k]
+        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]
+    end for k
+
+This algorihm is implemented by the function,
+
+    marsvin::Vector<T> forward_substitution(const ::marsvin::Matrix<T>& L, const ::marsvin::Vector<T>& b)
+
+For example,
+
+    // L : Nonsingular lower triangular matrix of nxn size.
+    // b : Vector of n size.
+    marsvin::Vector<T> x;
+    x = marsvin::forward_substitution(L, b);
+
+or,
+
+    // L : Nonsingular lower triangular matrix of nxn size.
+    // b : Vector of n size.
+    marsvin::Vector<T> 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 = 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]
+    end for k
+
+This algorihm is implemented by the function,
+
+    marsvin::Vector<T> forward_substitution_memory(const ::marsvin::Matrix<T>& L, ::marsvin::Vector<T>& b)
+
+For example,
+
+    // L : Nonsingular lower triangular matrix of nxn size.
+    // b : Vector of n size.
+    marsvin::forward_substitution_memory(L, b);
 
 ### 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<T> backward_substitution(const ::marsvin::Matrix<T>& L, const ::marsvin::Vector<T>& b)
+
+For example,
+
+    // U : Nonsingular upper triangular matrix of nxn size.
+    // b : Vector of n size.
+    marsvin::Vector<T> x;
+    x = marsvin::backward_substitution(L, b);
+
+or,
+
+    // U : Nonsingular upper triangular matrix of nxn size.
+    // b : Vector of n size.
+    marsvin::Vector<T> 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<T> backward_substitution_memory(const ::marsvin::Matrix<T>& U, ::marsvin::Vector<T>& b)
+
+For example,
+
+    // U : Nonsingular upper triangular matrix of nxn size.
+    // b : Vector of n size.
+    marsvin::backward_substitution_memory(L, b);
 
 ### Gauss Elimination
 

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
+)

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<float> U = {{3,2,1,1},
+                                {0,1,2,1},
+                                {0,0,1,1},
+                                {0,0,0,1}};
+    marsvin::Vector<float> b = {14,12,7,4};
+    marsvin::Vector<float> x_real = {1,2,3,4};
+    marsvin::Vector<float> 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<float> 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<float> b = {55,110,134,122,75};
+    marsvin::Vector<float> x_real = {1,2,3,4,5};
+    marsvin::Vector<float> 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/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<int> 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<int> m_rhs(rhs_kRows, rhs_kColumns, {1,2,3,4,5});
+    marsvin::BaseMatrix<int> 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));
+}
+

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) {