From 067e74a858c7b84b2232c74a28d6ce2fc283f26d Mon Sep 17 00:00:00 2001 From: Dharmin Shah Date: Thu, 12 Dec 2019 15:04:59 +0100 Subject: [PATCH 1/2] Assert doesn't work with std::numeric limits epsilon, so new eps manually set working gauss elimination with LU decomposition Signed-off-by: Dharmin Shah --- 11_12_2019/solve/solve.cpp | 24 +++++++--- Matrix/inc/MatrixUtils.hpp | 2 + Matrix/src/MatrixUtils.cpp | 90 +++++++++++++++++++++++++++++++++++++- 3 files changed, 110 insertions(+), 6 deletions(-) diff --git a/11_12_2019/solve/solve.cpp b/11_12_2019/solve/solve.cpp index 42f750e..6b1fcff 100644 --- a/11_12_2019/solve/solve.cpp +++ b/11_12_2019/solve/solve.cpp @@ -1,7 +1,8 @@ #include "Matrix.hpp" #include "MatrixUtils.hpp" - +#include #include +#include void testSolve( ); @@ -14,7 +15,17 @@ int main( ) void testSolve( ) { Matrix A( 3, 3 ); + //A( 0, 0 ) = 1.0; + //A( 0, 1 ) = 2.0; + //A( 0, 2 ) = -1.0; + + //A( 1, 0 ) = 4.0; + //A( 1, 1 ) = 3.0; + //A( 1, 2 ) = 1.0; + //A( 2, 0 ) = 2.0; + //A( 2, 1 ) = 2.0; + //A( 2, 2 ) = 3.0; A( 0, 0 ) = 3.0; A( 0, 1 ) = 2.0; A( 0, 2 ) = -4.0; @@ -27,12 +38,15 @@ void testSolve( ) A( 2, 1 ) = -3.0; A( 2, 2 ) = 1.0; - std::vector rhs = {3., 15., 14.}; + std::cout << A; + std::vector rhs = {3.0, 15.0, 14.0}; std::vector answer = matrixUtilities::solveByGaussElimination( A, rhs ); assert( answer.size() == 3 ); - assert( std::abs( answer[0] - 3.0 ) < std::numeric_limits::epsilon( ) ); - assert( std::abs( answer[1] - 1.0 ) < std::numeric_limits::epsilon( ) ); - assert( std::abs( answer[2] - 2.0 ) < std::numeric_limits::epsilon( ) ); + double eps = 1e-10; + //the numeric limits are of order 1e-16 so the assert will not work. + assert( std::abs( answer[0] - 3.0 ) < eps ); + assert( std::abs( answer[1] - 1.0 ) < eps ); + assert( std::abs( answer[2] - 2.0 ) < eps ); } diff --git a/Matrix/inc/MatrixUtils.hpp b/Matrix/inc/MatrixUtils.hpp index d46f4bd..30e8279 100644 --- a/Matrix/inc/MatrixUtils.hpp +++ b/Matrix/inc/MatrixUtils.hpp @@ -18,6 +18,8 @@ namespace matrixUtilities void pad(Matrix & m,int step=1, double value=0.0); void productThread(const Matrix & m1, const Matrix & m2, int begin, int end, int other_dim, Matrix & m3); void sumThread(const Matrix & m1, const Matrix & m2, int begin, int end, int other_dim, Matrix & m3); + void lu_generator(const Matrix&, Matrix&, Matrix&); + void gauss_lu_solver(const Matrix&, const Matrix&, std::vector&); std::vector solveByGaussElimination( const Matrix& matrix, const std::vector& rhs); //additional functions diff --git a/Matrix/src/MatrixUtils.cpp b/Matrix/src/MatrixUtils.cpp index bd60a80..380809a 100644 --- a/Matrix/src/MatrixUtils.cpp +++ b/Matrix/src/MatrixUtils.cpp @@ -21,7 +21,15 @@ namespace matrixUtilities std::vector solveByGaussElimination( const Matrix& matrix, const std::vector& rhs) { // Fill Gauss Elimination - return std::vector(); + std::vector b = rhs; + //Assign input matrix + Matrix A = matrix; + //L and U matrix generator + Matrix L(matrix.numberOfRows(), matrix.numberOfCols()), U(matrix.numberOfRows(), matrix.numberOfCols()); + lu_generator(matrix, L, U); + gauss_lu_solver(L, U, b); + + return b; } Matrix matrixSum(const Matrix &m1, const Matrix &m2) { @@ -143,6 +151,86 @@ void sumThread(const Matrix &m1, const Matrix &m2, int begin, int end, int other for (size_t i = begin; i < end; i++) for (size_t j = 0; j < other_dim; j++) m3(i, j) = m1(i, j) + m2(i, j); +} +void lu_generator(const Matrix &matrix, Matrix &L, Matrix &U) +{ + + const int rows = L.numberOfRows(); + const int cols = L.numberOfCols(); + Matrix A = matrix; + + for (int i = 0; i < rows; i++) + { + //Upper triangular matrix + for (int k = i; k < rows; k++) + { + double sum = 0.0; + for (int j = 0; j < i; j++) + { + sum += L(i, j)*U(j, k); + } + + U(i, k) = A(i, k) - sum; + } + + //Lower triangular matrix + for (int k = i; k < rows; k++) + { + if (i == k) + L(i, i) = 1.0; + else + { + double sum = 0.0; + for (int j = 0; j < i; j++) + sum += L(k, j)*U(j, i); + + //Evaluating for L + L(k, i) = (A(k, i) - sum) / U(i, i); + + } + + } + + + + } + + +} +void gauss_lu_solver(const Matrix &L, const Matrix &U, std::vector&b) +{ + std::vector z; + z.resize(b.size()); + b[0] = b[0] / L(0, 0); + //solve for Ly=z by forward substitution + for (int i = 1; i < L.numberOfRows(); i++) + { + double sum = 0.0; + for (int j = 0; j < i; j++) + { + sum += L(i, j)*b[j]; + } + + + b[i] = (b[i] - sum); + + } + + + //Solving for Ux=z by backward substitution + + + for (int i = L.numberOfRows() - 1; i >= 0; i--) + { + for (int j = i+1; j < L.numberOfRows(); ++j) + { + b[i] -= (U(i, j)*b[j]); + } + + b[i] = b[i] / U(i, i); + } + + } void pad(Matrix &m, int step, double value) { From f21a14e8e2c3d2479d3ec6ea7dba4bd7130da5cb Mon Sep 17 00:00:00 2001 From: Dharmin Shah Date: Thu, 12 Dec 2019 23:52:55 +0100 Subject: [PATCH 2/2] Added regular Gauss Signed-off-by: Dharmin Shah --- 11_12_2019/solve/solve.cpp | 21 ++++++++++------ Matrix/src/MatrixUtils.cpp | 50 +++++++++++++++++++++++++++++++++++--- 2 files changed, 60 insertions(+), 11 deletions(-) diff --git a/11_12_2019/solve/solve.cpp b/11_12_2019/solve/solve.cpp index 6b1fcff..3b3891e 100644 --- a/11_12_2019/solve/solve.cpp +++ b/11_12_2019/solve/solve.cpp @@ -42,11 +42,18 @@ void testSolve( ) std::vector rhs = {3.0, 15.0, 14.0}; std::vector answer = matrixUtilities::solveByGaussElimination( A, rhs ); - - assert( answer.size() == 3 ); - double eps = 1e-10; - //the numeric limits are of order 1e-16 so the assert will not work. - assert( std::abs( answer[0] - 3.0 ) < eps ); - assert( std::abs( answer[1] - 1.0 ) < eps ); - assert( std::abs( answer[2] - 2.0 ) < eps ); + //L and U matrix generator + Matrix L(A.numberOfRows(), A.numberOfCols()), U(A.numberOfRows(), A.numberOfCols()); + matrixUtilities::lu_generator(A, L, U); + matrixUtilities::gauss_lu_solver(L, U, rhs); + for (auto a : answer) + std::cout << a << " "; + for (auto b : rhs) + std::cout << b << " "; + // assert( answer.size() == 3 ); + //double eps = 1e-10; + ////the numeric limits are of order 1e-16 so the assert will not work. + // assert( std::abs( answer[0] - 3.0 ) < eps ); + // assert( std::abs( answer[1] - 1.0 ) < eps ); + // assert( std::abs( answer[2] - 2.0 ) < eps ); } diff --git a/Matrix/src/MatrixUtils.cpp b/Matrix/src/MatrixUtils.cpp index 380809a..ce80a6e 100644 --- a/Matrix/src/MatrixUtils.cpp +++ b/Matrix/src/MatrixUtils.cpp @@ -22,12 +22,54 @@ std::vector solveByGaussElimination( const Matrix& matrix, const std::ve { // Fill Gauss Elimination std::vector b = rhs; + std::vector result; + result.reserve(rhs.size()); //Assign input matrix Matrix A = matrix; - //L and U matrix generator - Matrix L(matrix.numberOfRows(), matrix.numberOfCols()), U(matrix.numberOfRows(), matrix.numberOfCols()); - lu_generator(matrix, L, U); - gauss_lu_solver(L, U, b); + + int n = A.numberOfRows(); + for(int i=0; i= 0; i--) + { + for (int j = i + 1; j < n; ++j) + { + b[i] -= (A(i, j)*b[j]); + } + + b[i] = b[i] / A(i, i); + } + return b; }