diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 978cdecc..7d96fca2 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -7,8 +7,8 @@
 ]]
 
 # Build portable randomized GMRES example
-add_executable(rand_gmres.exe randGmres.cpp)
-target_link_libraries(rand_gmres.exe PRIVATE ReSolve)
+add_executable(randGmres.exe randGmres.cpp)
+target_link_libraries(randGmres.exe PRIVATE ReSolve)
 
 # Build an example with a system GMRES solver
 add_executable(sysGmres.exe sysGmres.cpp)
@@ -48,9 +48,7 @@ endif(RESOLVE_USE_KLU)
 set(installable_executables "")
 
 # Install all examples in bin directory
-list(APPEND installable_executables rand_gmres.exe)
-
-list(APPEND installable_executables sysGmres.exe)
+list(APPEND installable_executables randGmres.exe sysGmres.exe)
 
 if(RESOLVE_USE_KLU)
   list(APPEND installable_executables kluFactor.exe kluRefactor.exe
diff --git a/examples/ExampleHelper.hpp b/examples/ExampleHelper.hpp
index 3a4fb260..ce3d5163 100644
--- a/examples/ExampleHelper.hpp
+++ b/examples/ExampleHelper.hpp
@@ -57,7 +57,7 @@ namespace ReSolve
        * @pre Workspace handles are initialized
        *
        * @post Handlers are instantiated.
-       * allocated
+       *
        */
       ExampleHelper(workspace_type& workspace)
         : mh_(&workspace),
@@ -92,11 +92,6 @@ namespace ReSolve
           delete res_;
           res_ = nullptr;
         }
-        if (x_true_)
-        {
-          delete x_true_;
-          x_true_ = nullptr;
-        }
       }
 
       /// Returns the configured hardware backend
@@ -111,55 +106,6 @@ namespace ReSolve
         return memspace_;
       }
 
-      /**
-       * @brief Set the new linear system together with its computed solution
-       * and compute solution error and residual norms.
-       *
-       * This will set the new system A*x = r and compute related error norms.
-       *
-       * @param A[in] - Linear system matrix
-       * @param r[in] - Linear system right-hand side
-       * @param x[in] - Computed solution of the linear system
-       */
-      void setSystem(ReSolve::matrix::Sparse* A,
-                     ReSolve::vector::Vector* r,
-                     ReSolve::vector::Vector* x)
-      {
-        assert((res_ == nullptr) && (x_true_ == nullptr));
-        A_   = A;
-        r_   = r;
-        x_   = x;
-        res_ = new ReSolve::vector::Vector(A->getNumRows());
-        computeNorms();
-      }
-
-      /**
-       * @brief Set the new linear system together with its computed solution
-       * and compute solution error and residual norms.
-       *
-       * This is to be used after values in A and r are updated.
-       *
-       * @todo This method probably does not need any input parameters.
-       *
-       * @param A[in] - Linear system matrix
-       * @param r[in] - Linear system right-hand side
-       * @param x[in] - Computed solution of the linear system
-       */
-      void resetSystem(ReSolve::matrix::Sparse* A,
-                       ReSolve::vector::Vector* r,
-                       ReSolve::vector::Vector* x)
-      {
-        A_ = A;
-        r_ = r;
-        x_ = x;
-        if (res_ == nullptr)
-        {
-          res_ = new ReSolve::vector::Vector(A->getNumRows());
-        }
-
-        computeNorms();
-      }
-
       /// Return L2 norm of the linear system residual.
       ReSolve::real_type getNormResidual()
       {
@@ -173,16 +119,60 @@ namespace ReSolve
       }
 
       /// Minimalistic summary
-      void printShortSummary()
+      void printShortSummary(ReSolve::matrix::Sparse* A,
+                             ReSolve::vector::Vector* r,
+                             ReSolve::vector::Vector* x)
       {
+        A_ = A;
+        r_ = r;
+        x_ = x;
+
+        if (res_ == nullptr)
+        {
+          res_ = new ReSolve::vector::Vector(A->getNumRows());
+        }
+        else
+        {
+          if (res_->getSize() != A->getNumRows())
+          {
+            delete res_;
+            res_ = new ReSolve::vector::Vector(A->getNumRows());
+          }
+        }
+
+        res_->copyDataFrom(r_, memspace_, memspace_);
+        real_type norm  = computeResidualNorm(*A_, *x_, *res_, memspace_);
+        real_type rnorm = norm2(*r_, memspace_);
+
         std::cout << "\t2-Norm of the residual: "
                   << std::scientific << std::setprecision(16)
-                  << getNormRelativeResidual() << "\n";
+                  << norm / rnorm << "\n";
       }
 
       /// Summary of direct solve
-      void printSummary()
+      void printSummary(ReSolve::matrix::Sparse* A,
+                        ReSolve::vector::Vector* r,
+                        ReSolve::vector::Vector* x)
       {
+        A_ = A;
+        r_ = r;
+        x_ = x;
+
+        if (res_ == nullptr)
+        {
+          res_ = new ReSolve::vector::Vector(A->getNumRows());
+        }
+        else
+        {
+          if (res_->getSize() != A->getNumRows())
+          {
+            delete res_;
+            res_ = new ReSolve::vector::Vector(A->getNumRows());
+          }
+        }
+
+        computeNorms();
+
         std::cout << "\t 2-Norm of the residual (before IR): "
                   << std::scientific << std::setprecision(16)
                   << getNormRelativeResidual() << "\n";
@@ -199,9 +189,9 @@ namespace ReSolve
       {
         std::cout << "FGMRES: init nrm: "
                   << std::scientific << std::setprecision(16)
-                  << ls->getInitResidualNorm() / norm_rhs_
+                  << ls->getInitResidualNorm()
                   << " final nrm: "
-                  << ls->getFinalResidualNorm() / norm_rhs_
+                  << ls->getFinalResidualNorm()
                   << " iter: " << ls->getNumIter() << "\n";
       }
 
@@ -209,33 +199,11 @@ namespace ReSolve
       void printIterativeSolverSummary(ReSolve::LinSolverIterative* ls)
       {
         std::cout << std::setprecision(16) << std::scientific;
-        std::cout << "\t Initial residual norm          ||b-A*x||       : " << ls->getInitResidualNorm() << "\n";
-        std::cout << "\t Initial relative residual norm ||b-A*x||/||b|| : " << ls->getInitResidualNorm() / norm_rhs_ << "\n";
-        std::cout << "\t Final residual norm            ||b-A*x||       : " << ls->getFinalResidualNorm() << "\n";
-        std::cout << "\t Final relative residual norm   ||b-A*x||/||b|| : " << ls->getFinalResidualNorm() / norm_rhs_ << "\n";
+        std::cout << "\t Initial relative residual norm ||b-A*x||/||b|| : " << ls->getInitResidualNorm() << "\n";
+        std::cout << "\t Final relative residual norm   ||b-A*x||/||b|| : " << ls->getFinalResidualNorm() << "\n";
         std::cout << "\t Number of iterations                           : " << ls->getNumIter() << "\n";
       }
 
-      /// Check the relative residual norm against `tolerance`.
-      int checkResult(ReSolve::real_type tolerance)
-      {
-        int                error_sum = 0;
-        ReSolve::real_type norm      = norm_res_ / norm_rhs_;
-
-        if (!std::isfinite(norm))
-        {
-          std::cout << "Result is not a finite number!\n";
-          error_sum++;
-        }
-        if (norm > tolerance)
-        {
-          std::cout << "Result inaccurate!\n";
-          error_sum++;
-        }
-
-        return error_sum;
-      }
-
       /**
        * @brief Verify the computation of the norm of scaled residuals.
        *
@@ -384,15 +352,14 @@ namespace ReSolve
       }
 
     private:
-      ReSolve::matrix::Sparse* A_; ///< pointer to system matrix
-      ReSolve::vector::Vector* r_; ///< pointer to system right-hand side
-      ReSolve::vector::Vector* x_; ///< pointer to the computed solution
+      ReSolve::matrix::Sparse* A_{nullptr}; ///< pointer to system matrix
+      ReSolve::vector::Vector* r_{nullptr}; ///< pointer to system right-hand side
+      ReSolve::vector::Vector* x_{nullptr}; ///< pointer to the computed solution
 
       ReSolve::MatrixHandler mh_; ///< matrix handler instance
       ReSolve::VectorHandler vh_; ///< vector handler instance
 
-      ReSolve::vector::Vector* res_{nullptr};    ///< pointer to residual vector
-      ReSolve::vector::Vector* x_true_{nullptr}; ///< pointer to solution error vector
+      ReSolve::vector::Vector* res_{nullptr}; ///< pointer to residual vector
 
       ReSolve::real_type norm_rhs_{0.0}; ///< right-hand side vector norm
       ReSolve::real_type norm_res_{0.0}; ///< residual vector norm
diff --git a/examples/experimental/CMakeLists.txt b/examples/experimental/CMakeLists.txt
index 5ad332a6..3fe1cb82 100644
--- a/examples/experimental/CMakeLists.txt
+++ b/examples/experimental/CMakeLists.txt
@@ -53,7 +53,6 @@ endif(RESOLVE_USE_KLU)
 set(installable_executables "")
 
 # Install all examples in bin directory
-list(APPEND installable_executables rand_gmres.exe)
 
 if(RESOLVE_USE_KLU)
 
diff --git a/examples/experimental/r_KLU_rocsolverrf_asym6x6.cpp b/examples/experimental/r_KLU_rocsolverrf_asym6x6.cpp
index 9a9d42e0..4ade639b 100644
--- a/examples/experimental/r_KLU_rocsolverrf_asym6x6.cpp
+++ b/examples/experimental/r_KLU_rocsolverrf_asym6x6.cpp
@@ -170,8 +170,7 @@ int main()
   }
   std::cout << "]" << std::endl;
 
-  helper.resetSystem(A, vec_rhs, vec_x);
-  helper.printShortSummary();
+  helper.printShortSummary(A, vec_rhs, vec_x);
   ReSolve::matrix::Csr* L = (ReSolve::matrix::Csr*) KLU.getLFactor();
   ReSolve::matrix::Csr* U = (ReSolve::matrix::Csr*) KLU.getUFactor();
   if (L == nullptr || U == nullptr)
@@ -214,7 +213,6 @@ int main()
       // Solve with refactorization
       status = Rf.solve(vec_rhs, vec_x);
       std::cout << "RocSolverRf solve status: " << status << std::endl;
-      helper.resetSystem(A, vec_rhs, vec_x);
 
       if (status == 0)
       {
diff --git a/examples/gluRefactor.cpp b/examples/gluRefactor.cpp
index 322b392b..6ba0f752 100644
--- a/examples/gluRefactor.cpp
+++ b/examples/gluRefactor.cpp
@@ -264,8 +264,7 @@ int gluRefactor(int argc, char* argv[])
     RESOLVE_RANGE_POP("Triangular solve");
 
     // Print summary of the results
-    helper.resetSystem(A, vec_rhs, vec_x);
-    helper.printSummary();
+    helper.printSummary(A, vec_rhs, vec_x);
 
   } // for (int i = 0; i < num_systems; ++i)
   RESOLVE_RANGE_POP(__FUNCTION__);
diff --git a/examples/gpuRefactor.cpp b/examples/gpuRefactor.cpp
index 90f5ee06..4b4b024b 100644
--- a/examples/gpuRefactor.cpp
+++ b/examples/gpuRefactor.cpp
@@ -260,8 +260,7 @@ int gpuRefactor(int argc, char* argv[])
       std::cout << "KLU solve status: " << status << std::endl;
 
       // Print summary of results
-      helper.resetSystem(A, vec_rhs, vec_x);
-      helper.printShortSummary();
+      helper.printShortSummary(A, vec_rhs, vec_x);
 
       if (i == 1)
       {
@@ -298,8 +297,7 @@ int gpuRefactor(int argc, char* argv[])
       RESOLVE_RANGE_POP("Refactorization");
 
       // Print summary of the results
-      helper.resetSystem(A, vec_rhs, vec_x);
-      helper.printSummary();
+      helper.printSummary(A, vec_rhs, vec_x);
 
       RESOLVE_RANGE_PUSH("Iterative refinement");
       if (is_iterative_refinement)
diff --git a/examples/kluFactor.cpp b/examples/kluFactor.cpp
index 78ad54d5..7703dbde 100644
--- a/examples/kluFactor.cpp
+++ b/examples/kluFactor.cpp
@@ -185,8 +185,7 @@ int main(int argc, char* argv[])
     status = KLU.solve(vec_rhs, vec_x);
     std::cout << "KLU solve status: " << status << std::endl;
 
-    helper.resetSystem(A, vec_rhs, vec_x);
-    helper.printShortSummary();
+    helper.printShortSummary(A, vec_rhs, vec_x);
     if (is_iterative_refinement)
     {
       // Setup iterative refinement
diff --git a/examples/kluRefactor.cpp b/examples/kluRefactor.cpp
index ddf188d4..559c1ad1 100644
--- a/examples/kluRefactor.cpp
+++ b/examples/kluRefactor.cpp
@@ -191,8 +191,7 @@ int main(int argc, char* argv[])
     status = KLU->solve(vec_rhs, vec_x);
     std::cout << "KLU solve status: " << status << std::endl;
 
-    helper.resetSystem(A, vec_rhs, vec_x);
-    helper.printShortSummary();
+    helper.printShortSummary(A, vec_rhs, vec_x);
     if (is_iterative_refinement)
     {
       // Setup iterative refinement
diff --git a/examples/randGmres.cpp b/examples/randGmres.cpp
index 629cc264..1bd0f62c 100644
--- a/examples/randGmres.cpp
+++ b/examples/randGmres.cpp
@@ -183,7 +183,6 @@ int runGmresExample(int argc, char* argv[])
   FGMRES.solve(vec_rhs, vec_x);
 
   // Print summary of results
-  helper.resetSystem(A, vec_rhs, vec_x);
   std::cout << "\nRandomized GMRES result on " << hwbackend << "\n";
   std::cout << "---------------------------------\n";
   helper.printIrSummary(&FGMRES);
diff --git a/examples/sysGmres.cpp b/examples/sysGmres.cpp
index 07085e09..faecbd01 100644
--- a/examples/sysGmres.cpp
+++ b/examples/sysGmres.cpp
@@ -268,8 +268,6 @@ int sysGmres(int argc, char* argv[])
 
   if (return_code == 0)
   {
-    helper.resetSystem(A, vec_rhs, vec_x);
-
     // Get reference to iterative solver and print results
     LinSolverIterative& iter_solver = solver.getIterativeSolver();
     helper.printIterativeSolverSummary(&iter_solver);
diff --git a/examples/sysRefactor.cpp b/examples/sysRefactor.cpp
index 3a095213..8ce61918 100644
--- a/examples/sysRefactor.cpp
+++ b/examples/sysRefactor.cpp
@@ -320,8 +320,7 @@ int sysRefactor(int argc, char* argv[])
     std::cout << "Triangular solve status: " << status << std::endl;
 
     // Print summary of results
-    helper.resetSystem(A, vec_rhs, vec_x);
-    helper.printShortSummary();
+    helper.printShortSummary(A, vec_rhs, vec_x);
     if ((i > 1) && is_iterative_refinement)
     {
       helper.printIrSummary(&(solver.getIterativeSolver()));
diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp
index 643f23a4..25d51c51 100644
--- a/resolve/LinSolverIterativeFGMRES.cpp
+++ b/resolve/LinSolverIterativeFGMRES.cpp
@@ -138,7 +138,7 @@ namespace ReSolve
     io::Logger::misc() << "it 0: norm of residual "
                        << std::scientific << std::setprecision(16)
                        << rnorm << " Norm of rhs: " << bnorm << "\n";
-    initial_residual_norm_ = rnorm;
+    initial_residual_norm_ = rnorm / bnorm; // relative residual norm
     while (outer_flag)
     {
       if (it == 0)
@@ -307,7 +307,7 @@ namespace ReSolve
 
       if (!outer_flag)
       {
-        final_residual_norm_ = rnorm;
+        final_residual_norm_ = rnorm / bnorm; // relative residual norm
         total_iters_         = it;
         io::Logger::misc() << "End of cycle, COMPUTED norm of residual "
                            << std::scientific << std::setprecision(16)
diff --git a/resolve/LinSolverIterativeRandFGMRES.cpp b/resolve/LinSolverIterativeRandFGMRES.cpp
index db26b007..8d315536 100644
--- a/resolve/LinSolverIterativeRandFGMRES.cpp
+++ b/resolve/LinSolverIterativeRandFGMRES.cpp
@@ -181,7 +181,7 @@ namespace ReSolve
                        << std::scientific << std::setprecision(16)
                        << rnorm << " Norm of rhs: " << bnorm << "\n";
 
-    initial_residual_norm_ = rnorm;
+    initial_residual_norm_ = rnorm / bnorm; // compute relative residual norm
     while (outer_flag)
     {
       if (it == 0)
@@ -399,7 +399,7 @@ namespace ReSolve
                            << std::scientific << std::setprecision(16)
                            << rnorm << "\n";
 
-        final_residual_norm_ = rnorm;
+        final_residual_norm_ = rnorm / bnorm; // relative residual norm
         total_iters_         = it;
       }
     } // outer while
diff --git a/tests/functionality/TestHelper.hpp b/tests/functionality/TestHelper.hpp
index 85e7a1dc..4db0bac0 100644
--- a/tests/functionality/TestHelper.hpp
+++ b/tests/functionality/TestHelper.hpp
@@ -362,7 +362,7 @@ class TestHelper
    * the residual norm for the system that has been set by the constructor
    * or (re)setSystem functions.
    *
-   * @param rrn_system - residual norm value to be verified
+   * @param rn_system - residual norm value to be verified
    * @return int - 0 if the result is correct, error code otherwise
    */
   int checkResidualNorm(ReSolve::real_type rn_system)
diff --git a/tests/functionality/testSysGmres.cpp b/tests/functionality/testSysGmres.cpp
index 1e45ee1b..277c1842 100644
--- a/tests/functionality/testSysGmres.cpp
+++ b/tests/functionality/testSysGmres.cpp
@@ -182,7 +182,7 @@ int test(int argc, char* argv[])
             << "\t Solver tolerance:               " << tol_out << "\n";
   helper.printIterativeSolverSummary(&(solver.getIterativeSolver()));
 
-  error_sum += helper.checkResidualNorm(solver.getIterativeSolver().getFinalResidualNorm());
+  error_sum += helper.checkRelativeResidualNorm(solver.getIterativeSolver().getFinalResidualNorm());
   error_sum += helper.checkResult(10.0 * tol_out);
   isTestPass(error_sum, "Test");