ENH: Add a flag to compute the residual

Add a flag to computer the residual for iterative methods
changes to SOR and Gauss-Seidel methods

Solvers/BackwardGaussSeidel.cpp

@@ -24,7 +24,9 @@
 #include "Solvers/BackwardGaussSeidel.h"
 #include <memory>
 
-BackwardGaussSeidel::BackwardGaussSeidel(const core::SparseMatrixd &A, const core::Vectord &rhs):
+BackwardGaussSeidel::BackwardGaussSeidel(
+    const core::SparseMatrixd &A,
+    const core::Vectord &rhs):
     U(A.triangularView<Eigen::Upper>()),
     L(A.triangularView<Eigen::StrictlyLower>())
 {
@@ -32,9 +34,13 @@ BackwardGaussSeidel::BackwardGaussSeidel(const core::SparseMatrixd &A, const cor
 }
 
 //---------------------------------------------------------------------------
-void BackwardGaussSeidel::iterate(core::Vectord &x)
+void BackwardGaussSeidel::iterate(core::Vectord &x, bool updateResidual)
 {
     x = this->linearSystem->getRHSVector() - L * x;
     U.solveInPlace(x);
-    this->linearSystem->computeResidual(x, this->residual);
+
+    if (updateResidual)
+    {
+        this->linearSystem->computeResidual(x, this->residual);
+    }
 }

Solvers/BackwardGaussSeidel.h

@@ -49,7 +49,7 @@ public:
     ///
     /// \brief Do one iteration of the method
     ///
-    void iterate(core::Vectord &x) override;
+    void iterate(core::Vectord &x, bool updateResidual = true) override;
 
 private:
     Eigen::SparseTriangularView<core::SparseMatrixd, Eigen::Upper> U;

Solvers/BackwardSOR.cpp

@@ -29,10 +29,10 @@ BackwardSOR::BackwardSOR(const core::SparseMatrixd &A,
 {}
 
 //---------------------------------------------------------------------------
-void BackwardSOR::iterate(core::Vectord &x)
+void BackwardSOR::iterate(core::Vectord &x, bool updateResidual)
 {
     auto old = x; // necessary copy
-    this->gaussSeidel.iterate(x);
+    this->gaussSeidel.iterate(x, updateResidual);
     x = this->weight * x + (1 - this->weight) * old;
 }
 

Solvers/BackwardSOR.h

@@ -43,7 +43,7 @@ public:
     ///
     /// \brief Do one iteration of the method
     ///
-    void iterate(core::Vectord &x) override;
+    void iterate(core::Vectord &x, bool updateResidual = true) override;
 
     ///
     /// \brief Set Weight

Solvers/ForwardGaussSeidel.cpp

@@ -23,7 +23,9 @@
 
 #include "Solvers/ForwardGaussSeidel.h"
 
-ForwardGaussSeidel::ForwardGaussSeidel(const core::SparseMatrixd &A, const core::Vectord &rhs):
+ForwardGaussSeidel::ForwardGaussSeidel(
+    const core::SparseMatrixd &A,
+    const core::Vectord &rhs):
     U(A.triangularView<Eigen::StrictlyUpper>()),
     L(A.triangularView<Eigen::Lower>())
 {
@@ -31,9 +33,13 @@ ForwardGaussSeidel::ForwardGaussSeidel(const core::SparseMatrixd &A, const core:
 }
 
 //---------------------------------------------------------------------------
-void ForwardGaussSeidel::iterate(core::Vectord &x)
+void ForwardGaussSeidel::iterate(core::Vectord &x, bool updateResidual)
 {
     x = this->linearSystem->getRHSVector() - U * x;
     L.solveInPlace(x);
-    this->linearSystem->computeResidual(x, this->residual);
+
+    if (updateResidual)
+    {
+        this->linearSystem->computeResidual(x, this->residual);
+    }
 }

Solvers/ForwardGaussSeidel.h

@@ -48,7 +48,7 @@ public:
     ///
     /// \brief Do one iteration of the method
     ///
-    void iterate(core::Vectord &x) override;
+    void iterate(core::Vectord &x, bool updateResidual = true) override;
 
 private:
     Eigen::SparseTriangularView<core::SparseMatrixd, Eigen::StrictlyUpper> U;

Solvers/ForwardSOR.cpp

@@ -29,10 +29,10 @@ ForwardSOR::ForwardSOR(const core::SparseMatrixd &A,
 {}
 
 //---------------------------------------------------------------------------
-void ForwardSOR::iterate(core::Vectord &x)
+void ForwardSOR::iterate(core::Vectord &x, bool updateResidual)
 {
     auto old = x; // necessary copy
-    this->gaussSeidel.iterate(x);
+    this->gaussSeidel.iterate(x, updateResidual);
     x = this->weight * x + (1 - this->weight) * old;
 }
 

Solvers/ForwardSOR.h

@@ -49,7 +49,7 @@ public:
     ///
     /// \brief Do one iteration of the method
     ///
-    void iterate(core::Vectord &x) override;
+    void iterate(core::Vectord &x, bool updateResidual = true) override;
 
     ///
     /// \brief Set Weight

Solvers/IterativeLinearSolver.h

@@ -42,7 +42,7 @@ public:
     ///
     /// \brief Do one iteration of the method
     ///
-    virtual void iterate(core::Vectord &x) = 0;
+    virtual void iterate(core::Vectord &x, bool updateResidual = true) = 0;
 
     ///
     /// \brief Solve the linear system using Gauss-Seidel iterations