Merge pull request #18 from solarcalf/update_accuracy

Update accuracy

Author: Amazingkivas

Poisson_Equation/headers/Solvers.hpp

@@ -60,6 +60,18 @@ namespace numcpp
         return norm;
     }
 
+    FP error(const std::vector<FP>& v1, const std::vector<FP>& v2)
+    {
+        FP error = 0.0;
+
+#pragma omp parallel for reduction(max: error)
+        for (size_t i = 0; i < v1.size(); ++i)
+        {
+            error = std::max(std::abs(v1[i] - v2[i]), error);
+        }
+        return error;
+    }
+
     std::vector<FP> vector_FMA(const std::vector<FP>& lhs, FP coef, const std::vector<FP>& rhs)
     {
         const size_t size = lhs.size();
@@ -210,15 +222,19 @@ namespace numcpp
 
             for (size_t i = 0; i < max_iterations; ++i)
             {
+                std::vector<FP> saved_approximation = approximation;
+
                 std::vector<FP> residual = (*system_matrix) * approximation - b;
 
-                FP residual_norm = norm(residual);
-                if (residual_norm <= required_precision) break;
+                //FP residual_norm = norm(residual);
 
                 std::vector<FP> Ar = (*system_matrix) * residual;
                 FP tau = scalar_product(Ar, residual) / scalar_product(Ar, Ar);
 
                 approximation = vector_FMA(residual, -tau, approximation);
+
+                FP approximation_error = error(saved_approximation, approximation);
+                if (approximation_error <= required_precision) break;
             }
 
             return approximation;
@@ -439,10 +455,11 @@ namespace numcpp
 
             for (size_t i = 1; i < max_iterations; ++i)
             {
+                std::vector<FP> saved_approximation = approximation;
+
                 residual = (*system_matrix) * approximation - b;
 
-                FP residual_norm = norm(residual);
-                if (residual_norm <= required_precision) break;
+                //FP residual_norm = norm(residual);
 
                 FP beta = scalar_product(Ah, residual) / scalar_product(Ah, h);
 
@@ -453,6 +470,9 @@ namespace numcpp
                 alpha = -scalar_product(residual, h) / scalar_product(Ah, h);
 
                 approximation = vector_FMA(h, alpha, approximation);
+
+                FP approximation_error = error(saved_approximation, approximation);
+                if (approximation_error <= required_precision) break;
             }
 
             return approximation;

Poisson_Equation/src/test.cpp

@@ -120,10 +120,15 @@ void test_task_custom_grid(size_t m, size_t n, std::unique_ptr<numcpp::ISolver>
 
 int main() 
 {
-    size_t m = 1024;
-    size_t n = 1024;
+    size_t m1 = 1024;
+    size_t n1 = 1024;
 
-    auto LS_solver = std::make_unique<numcpp::ConGrad>(std::vector<FP>(), 1000000, 0.00001, nullptr, std::vector<FP>()); 
+    auto LS_solver1 = std::make_unique<numcpp::ConGrad>(std::vector<FP>(), 1000000, 0.000000001, nullptr, std::vector<FP>());
+    test_task_custom_grid(m1, n1, std::move(LS_solver1));
 
-    test_task_custom_grid(m, n, std::move(LS_solver));
+    size_t m2 = 200;
+    size_t n2 = 200;
+
+    auto LS_solver2 = std::make_unique<numcpp::MinRes>(std::vector<FP>(), 1000000, 0.000001, nullptr, std::vector<FP>());
+    test_task(m2, n2, std::move(LS_solver2));
 }