diff --git a/11_12_2019/solve/solve.cpp b/11_12_2019/solve/solve.cpp
index 42f750e..3b3891e 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<iostream>
#include <cassert>
+#include<typeinfo>
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,22 @@ void testSolve( )
A( 2, 1 ) = -3.0;
A( 2, 2 ) = 1.0;
- std::vector<double> rhs = {3., 15., 14.};
+ std::cout << A;
+ std::vector<double> rhs = {3.0, 15.0, 14.0};
std::vector<double> answer = matrixUtilities::solveByGaussElimination( A, rhs );
-
- assert( answer.size() == 3 );
- assert( std::abs( answer[0] - 3.0 ) < std::numeric_limits<double>::epsilon( ) );
- assert( std::abs( answer[1] - 1.0 ) < std::numeric_limits<double>::epsilon( ) );
- assert( std::abs( answer[2] - 2.0 ) < std::numeric_limits<double>::epsilon( ) );
+ //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/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<double>&);
std::vector<double> solveByGaussElimination( const Matrix& matrix, const std::vector<double>& rhs);
//additional functions
diff --git a/Matrix/src/MatrixUtils.cpp b/Matrix/src/MatrixUtils.cpp
index bd60a80..ce80a6e 100644
--- a/Matrix/src/MatrixUtils.cpp
+++ b/Matrix/src/MatrixUtils.cpp
@@ -21,7 +21,57 @@ namespace matrixUtilities
std::vector<double> solveByGaussElimination( const Matrix& matrix, const std::vector<double>& rhs)
{
// Fill Gauss Elimination
- return std::vector<double>();
+ std::vector<double> b = rhs;
+ std::vector<double> result;
+ result.reserve(rhs.size());
+ //Assign input matrix
+ Matrix A = matrix;
+
+ int n = A.numberOfRows();
+ for(int i=0; i<n; i++)
+ for(int j=i+1; j<n;j++)
+ if (abs(A(i, i) < abs(A(j, i))))
+ {
+ for (int k = 0; k < n; k++)
+ {
+ A(i, k) = A(i, k) + A(j, k);
+ A(j, k) = A(i, k) - A(j, k);
+ A(i, k) = A(i, k) - A(j, k);
+
+ }
+
+ b[i] = b[i] + b[j];
+ b[j] = b[i] - b[j];
+ b[i] = b[i] - b[j];
+ }
+
+ std::cout << A;
+
+ //performing Gaussian elimination
+ for(int i=0; i<n;i++)
+ for (int j = i + 1; j < n; j++)
+ {
+ double f = A(j, i) / A(i, i);
+ for (int k = 0; k < n; k++)
+ {
+ A(j, k) = A(j, k) - f * A(i, k);
+ }
+ b[j] = b[j] - f * b[i];
+
+ }
+ //backward substitution
+ for (int i = n - 1; 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;
}
Matrix matrixSum(const Matrix &m1, const Matrix &m2)
{
@@ -143,6 +193,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<double>&b)
+{
+ std::vector<double> 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)
{