From af9dfe6f8ca866a7f29ba88cc90fc46740ce0224 Mon Sep 17 00:00:00 2001 From: Kacper Donat Date: Wed, 2 May 2018 23:29:41 +0200 Subject: [PATCH] solvers --- Makefile | 18 +++++++ main.cpp | 22 ++++++-- matrix.h | 40 +++++++++++---- solver.h | 149 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 216 insertions(+), 13 deletions(-) create mode 100644 Makefile create mode 100644 solver.h diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..371c068 --- /dev/null +++ b/Makefile @@ -0,0 +1,18 @@ +cpp=cl +link=link + +cppflags=/EHsc /std:c++latest /nologo /O2 +linkflags=/nologo + +all: solver.exe + +.cpp.obj: + $(cpp) /c $(cppdebug) $(cppflags) $(cppvars) $*.cpp + +solver.exe: main.obj + $(link) $(ldebug) $(conflags) /out:$@ $** $(conlibs) + +clean: + del *.obj + del *.lst + del *.exe diff --git a/main.cpp b/main.cpp index a7b6893..ee6f783 100644 --- a/main.cpp +++ b/main.cpp @@ -1,10 +1,12 @@ #include "matrix.h" +#include "solver.h" #include #include const size_t N = 10; -std::pair, Matrix> prepare(size_t index, size_t n) { +std::pair, Matrix> prepare(size_t index, size_t n) +{ double a1, a2, a3; a2 = a3 = -1; a1 = 5 + (index / 100) % 10; @@ -21,9 +23,23 @@ std::pair, Matrix> prepare(size_t index, size_t n) { return std::make_pair(M, b); } -int main(int argc, const char argv[]) +int main(int argc, const char* argv[]) { auto tuple = prepare(165581, N); auto M = std::get<0>(tuple); auto b = std::get<1>(tuple); -} \ No newline at end of file + + JacobiSolver jacobi; + std::cout << jacobi.solve(M, b); + + GaussSolver gauss; + std::cout << gauss.solve(M, b); + + LUSolver lu; + + auto pair = LUSolver::decompose(M); + auto L = std::get<0>(pair); + auto U = std::get<1>(pair); + + std::cout << lu.solve(L, U, b); +} diff --git a/matrix.h b/matrix.h index 36936f4..6352dfc 100644 --- a/matrix.h +++ b/matrix.h @@ -29,11 +29,14 @@ public: T& operator()(size_t n, size_t m); T operator()(size_t n, size_t m) const; - self_t operator*(const self_t& rhs); - self_t operator+(const self_t& rhs); + self_t operator*(const self_t& rhs) const; + self_t operator+(const self_t& rhs) const; - self_t operator*(T rhs); - self_t operator+(T rhs); + self_t operator*(T rhs) const; + self_t operator+(T rhs) const; + + self_t operator-() const; + self_t operator-(const self_t& rhs) const; self_t& operator=(self_t&& rhs) noexcept; @@ -75,7 +78,7 @@ T Matrix::operator()(size_t n, size_t m) const { } template -Matrix Matrix::operator+(const Matrix &rhs) { +Matrix Matrix::operator+(const Matrix &rhs) const { assert(this->_cols == rhs._cols && this->_rows == rhs._rows); self_t result(this->_rows, this->_cols); @@ -87,7 +90,7 @@ Matrix Matrix::operator+(const Matrix &rhs) { } template -Matrix Matrix::operator*(const Matrix &rhs) { +Matrix Matrix::operator*(const Matrix &rhs) const { assert(_cols == rhs._rows); self_t result(_rows, rhs._cols); @@ -105,23 +108,40 @@ Matrix Matrix::operator*(const Matrix &rhs) { } template -Matrix Matrix::operator*(const T rhs) { +Matrix Matrix::operator-() const { + self_t result(*this); + return result * -1; +} + +template +Matrix Matrix::operator-(const Matrix& rhs) const { self_t result(*this); for (size_t i = 0; i < _rows; i++) for (size_t j = 0; j < _cols; j++) - result(i, j) *= rhs; + result(i, j) -= rhs(i, j); return result; } template -Matrix Matrix::operator+(const T rhs) { +Matrix Matrix::operator*(const T rhs) const { self_t result(*this); for (size_t i = 0; i < _rows; i++) for (size_t j = 0; j < _cols; j++) - result(i, j) += rhs; + result(i, j) *= rhs(i, j); + + return result; +} + +template +Matrix Matrix::operator+(const T rhs) const { + self_t result(*this); + + for (size_t i = 0; i < _rows; i++) + for (size_t j = 0; j < _cols; j++) + result(i, j) += rhs(i, j); return result; } diff --git a/solver.h b/solver.h new file mode 100644 index 0000000..8b0df6d --- /dev/null +++ b/solver.h @@ -0,0 +1,149 @@ +#ifndef P02_SOLVER_H +#define P02_SOLVER_H + +#include "matrix.h" + +template +double residuum(const Matrix& M, const Matrix& b, const Matrix& x) +{ + return norm(M*x - b); +} + +class Solver { +public: + virtual Matrix solve(const Matrix& M, const Matrix& b) = 0; +}; + +class LUSolver : public Solver { +public: + static std::pair, Matrix> decompose(const Matrix& M); + Matrix solve(const Matrix& L, const Matrix& U, const Matrix& b); + Matrix solve(const Matrix& M, const Matrix& b) override; +}; + +class IterativeSolver : public Solver { +protected: + double error; +public: + IterativeSolver(); + IterativeSolver(double error); + + virtual Matrix iteration(const Matrix& M, const Matrix& b, const Matrix& old) = 0; + virtual Matrix solve(const Matrix& M, const Matrix& b) override; +}; + +class JacobiSolver : public IterativeSolver { +public: + using IterativeSolver::IterativeSolver; + virtual Matrix iteration(const Matrix& M, const Matrix& b, const Matrix& old) override; +}; + +class GaussSolver : public IterativeSolver { +public: + using IterativeSolver::IterativeSolver; + virtual Matrix iteration(const Matrix& M, const Matrix& b, const Matrix& old) override; +}; + +Matrix IterativeSolver::solve(const Matrix& M, const Matrix& b) +{ + assert(M.rows() == M.cols()); + assert(b.rows() == M.cols()); + + size_t n = M.cols(); + + Matrix x(n, 1); + for (size_t i = 0; i < n; ++i) x(i, 0) = 1.; + + size_t iterations = 0; + while(residuum(M, b, x) > error) { + x = iteration(M, b, x); + iterations ++; + } + + return x; +} + +IterativeSolver::IterativeSolver() : error(1e-9) {} +IterativeSolver::IterativeSolver(double error) : error(error) {} + +Matrix JacobiSolver::iteration(const Matrix& M, const Matrix& b, const Matrix& old) +{ + Matrix x(b.rows(), 1); + + for (size_t i = 0; i < b.rows(); ++i) { + double t = b(i, 0); + for (size_t j = 0; j < b.rows(); ++j) { + if (j == i) continue; + t -= M(i, j)*old(j, 0); + } + + x(i, 0) = t / M(i,i); + } + + return x; +} + +Matrix GaussSolver::iteration(const Matrix& M, const Matrix& b, const Matrix& old) +{ + Matrix x(b.rows(), 1); + + for (size_t i = 0; i < b.rows(); ++i) { + double t = b(i, 0); + + for (size_t j = 0; j < i; ++j) + t -= M(i, j)*x(j, 0); + for (size_t j = i+1; j < b.rows(); ++j) + t -= M(i, j)*old(j, 0); + + x(i, 0) = t / M(i,i); + } + + return x; +} + +Matrix LUSolver::solve(const Matrix& M, const Matrix& b) +{ + auto LU = LUSolver::decompose(M); + return this->solve(std::get<0>(LU), std::get<1>(LU), b); +} + +Matrix LUSolver::solve(const Matrix& L, const Matrix& U, const Matrix& b) +{ + Matrix x(b); + const size_t n = x.rows(); + + for (size_t i = 0; i < n; i++) { + for (size_t j = 0; j < i; j++) { + x(i, 0) -= L(i,j)*x(j, 0); + } + } + + for (size_t i = n; i != 0; i--) { + for (size_t j = n-1; j > i-1; j--) { + x(i-1, 0) -= U(i-1,j)*x(j, 0); + } + x(i-1, 0) /= U(i-1,i-1); + } + + return x; +} + +std::pair, Matrix> LUSolver::decompose(const Matrix& M) +{ + const size_t n = M.rows(); + + Matrix U(M); + Matrix L = Matrix::diag(n, 1); + + for (size_t i = 0; i < n - 1; ++i) { + for (size_t j = i+1; j < n; ++j) { + L(j, i) = U(j, i) / U(i, i); + for (size_t k = i; k < n; k++) { + U(j, k) -= L(j, i)*U(i, k); + } + } + } + + return std::make_pair(L, U); +} +#endif