solvers

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

main.cpp

@@ -1,10 +1,12 @@
 #include "matrix.h"
+#include "solver.h"
 #include <iostream>
 #include <utility>
 
 const size_t N = 10;
 
-std::pair<Matrix<double>, Matrix<double>> prepare(size_t index, size_t n) {
+std::pair<Matrix<double>, Matrix<double>> 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<double>, Matrix<double>> 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);
-}
+
+    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);
+}

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) const;
-    self_t operator+(const self_t& rhs);
+    self_t operator+(const self_t& rhs) const;
 
-    self_t operator*(T rhs);
+    self_t operator*(T rhs) const;
-    self_t operator+(T rhs);
+    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<T>::operator()(size_t n, size_t m) const {
 }
 
 template<typename T>
-Matrix<T> Matrix<T>::operator+(const Matrix<T> &rhs) {
+Matrix<T> Matrix<T>::operator+(const Matrix<T> &rhs) const {
     assert(this->_cols == rhs._cols && this->_rows == rhs._rows);
     self_t result(this->_rows, this->_cols);
 
@@ -87,7 +90,7 @@ Matrix<T> Matrix<T>::operator+(const Matrix<T> &rhs) {
 }
 
 template<typename T>
-Matrix<T> Matrix<T>::operator*(const Matrix<T> &rhs) {
+Matrix<T> Matrix<T>::operator*(const Matrix<T> &rhs) const {
     assert(_cols == rhs._rows);
     self_t result(_rows, rhs._cols);
 
@@ -105,23 +108,40 @@ Matrix<T> Matrix<T>::operator*(const Matrix<T> &rhs) {
 }
 
 template<typename T>
-Matrix<T> Matrix<T>::operator*(const T rhs) {
+Matrix<T> Matrix<T>::operator-() const {
+    self_t result(*this);
+    return result * -1;
+}
+
+template<typename T>
+Matrix<T> Matrix<T>::operator-(const Matrix<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<typename T>
-Matrix<T> Matrix<T>::operator+(const T rhs) {
+Matrix<T> Matrix<T>::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<typename T>
+Matrix<T> Matrix<T>::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;
 }

solver.h

@@ -0,0 +1,149 @@
+#ifndef P02_SOLVER_H
+#define P02_SOLVER_H
+
+#include "matrix.h"
+
+template <typename T>
+double residuum(const Matrix<T>& M, const Matrix<T>& b, const Matrix<T>& x)
+{
+    return norm(M*x - b);
+}
+
+class Solver {
+public:
+    virtual Matrix<double> solve(const Matrix<double>& M, const Matrix<double>& b) = 0;
+};
+
+class LUSolver : public Solver {
+public: 
+    static std::pair<Matrix<double>, Matrix<double>> decompose(const Matrix<double>& M);
+    Matrix<double> solve(const Matrix<double>& L, const Matrix<double>& U, const Matrix<double>& b);
+    Matrix<double> solve(const Matrix<double>& M, const Matrix<double>& b) override;
+};
+
+class IterativeSolver : public Solver {
+protected:
+    double error;
+public:
+    IterativeSolver();
+    IterativeSolver(double error);   
+
+    virtual Matrix<double> iteration(const Matrix<double>& M, const Matrix<double>& b, const Matrix<double>& old) = 0;
+    virtual Matrix<double> solve(const Matrix<double>& M, const Matrix<double>& b) override;
+};
+
+class JacobiSolver : public IterativeSolver {
+public:  
+    using IterativeSolver::IterativeSolver;
+    virtual Matrix<double> iteration(const Matrix<double>& M, const Matrix<double>& b, const Matrix<double>& old) override;
+};
+
+class GaussSolver : public IterativeSolver {
+public:  
+    using IterativeSolver::IterativeSolver;
+    virtual Matrix<double> iteration(const Matrix<double>& M, const Matrix<double>& b, const Matrix<double>& old) override;
+};
+
+Matrix<double> IterativeSolver::solve(const Matrix<double>& M, const Matrix<double>& b)
+{
+    assert(M.rows() == M.cols());
+    assert(b.rows() == M.cols());
+
+    size_t n = M.cols();
+
+    Matrix<double> 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<double> JacobiSolver::iteration(const Matrix<double>& M, const Matrix<double>& b, const Matrix<double>& old)
+{
+    Matrix<double> 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<double> GaussSolver::iteration(const Matrix<double>& M, const Matrix<double>& b, const Matrix<double>& old)
+{
+    Matrix<double> 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<double> LUSolver::solve(const Matrix<double>& M, const Matrix<double>& b)
+{
+    auto LU = LUSolver::decompose(M);
+    return this->solve(std::get<0>(LU), std::get<1>(LU), b);
+}
+
+Matrix<double> LUSolver::solve(const Matrix<double>& L, const Matrix<double>& U, const Matrix<double>& b)
+{
+    Matrix<double> 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<double>, Matrix<double>> LUSolver::decompose(const Matrix<double>& M)
+{
+    const size_t n = M.rows();
+
+    Matrix<double> U(M);
+    Matrix<double> L = Matrix<double>::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