From 4ded1ceeb0703bf208afe112327338bb7d68820f Mon Sep 17 00:00:00 2001
From: Daniel Richtmann <daniel.richtmann@gmail.com>
Date: Thu, 8 Feb 2018 15:29:44 +0100
Subject: [PATCH] Make GMRES solvers perform no more than MaxIterations steps

I noticed that it was possible to overrun this number.
---
 .../CommunicationAvoidingGeneralisedMinimalResidual.h       | 6 ++++--
 ...lexibleCommunicationAvoidingGeneralisedMinimalResidual.h | 6 ++++--
 .../iterative/FlexibleGeneralisedMinimalResidual.h          | 6 ++++--
 lib/algorithms/iterative/GeneralisedMinimalResidual.h       | 6 ++++--
 4 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h b/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h
index 3bf19131..b56c8afc 100644
--- a/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h
+++ b/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h
@@ -97,7 +97,9 @@ class CommunicationAvoidingGeneralisedMinimalResidual : public OperatorFunction<
     SolverTimer.Start();
 
     IterationCount = 0;
-    for (int k=0; k<MaxIterations; k++) {
+    auto outerLoopMax = MaxIterations/RestartLength + ((MaxIterations%RestartLength == 0) ? 0 : 1);
+
+    for (int k=0; k<outerLoopMax; k++) {
 
       cp = outerLoopBody(LinOp, src, psi, rsq);
 
@@ -167,7 +169,7 @@ class CommunicationAvoidingGeneralisedMinimalResidual : public OperatorFunction<
       std::cout << GridLogIterative << "CommunicationAvoidingGeneralisedMinimalResidual: Iteration " << IterationCount
                 << " residual " << cp << " target " << rsq << std::endl;
 
-      if ((i == RestartLength - 1) || (cp <= rsq)) {
+      if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) {
 
         computeSolution(v, psi, i);
 
diff --git a/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h b/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h
index 77936f29..057c1af3 100644
--- a/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h
+++ b/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h
@@ -103,7 +103,9 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF
     SolverTimer.Start();
 
     IterationCount = 0;
-    for (int k=0; k<MaxIterations; k++) {
+    auto outerLoopMax = MaxIterations/RestartLength + ((MaxIterations%RestartLength == 0) ? 0 : 1);
+
+    for (int k=0; k<outerLoopMax; k++) {
 
       cp = outerLoopBody(LinOp, src, psi, rsq);
 
@@ -179,7 +181,7 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF
       std::cout << GridLogIterative << "FlexibleCommunicationAvoidingGeneralisedMinimalResidual: Iteration " << IterationCount
                 << " residual " << cp << " target " << rsq << std::endl;
 
-      if ((i == RestartLength - 1) || (cp <= rsq)) {
+      if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) {
 
         computeSolution(v, psi, i);
 
diff --git a/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h b/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h
index 906dd35d..06aded09 100644
--- a/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h
+++ b/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h
@@ -103,7 +103,9 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction<Field> {
     SolverTimer.Start();
 
     IterationCount = 0;
-    for (int k=0; k<MaxIterations; k++) {
+    auto outerLoopMax = MaxIterations/RestartLength + ((MaxIterations%RestartLength == 0) ? 0 : 1);
+
+    for (int k=0; k<outerLoopMax; k++) {
 
       cp = outerLoopBody(LinOp, src, psi, rsq);
 
@@ -179,7 +181,7 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction<Field> {
       std::cout << GridLogIterative << "FlexibleGeneralisedMinimalResidual: Iteration " << IterationCount
                 << " residual " << cp << " target " << rsq << std::endl;
 
-      if ((i == RestartLength - 1) || (cp <= rsq)) {
+      if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) {
 
         computeSolution(v, psi, i);
 
diff --git a/lib/algorithms/iterative/GeneralisedMinimalResidual.h b/lib/algorithms/iterative/GeneralisedMinimalResidual.h
index 78dbf5ca..e9d54a53 100644
--- a/lib/algorithms/iterative/GeneralisedMinimalResidual.h
+++ b/lib/algorithms/iterative/GeneralisedMinimalResidual.h
@@ -97,7 +97,9 @@ class GeneralisedMinimalResidual : public OperatorFunction<Field> {
     SolverTimer.Start();
 
     IterationCount = 0;
-    for (int k=0; k<MaxIterations; k++) {
+    auto outerLoopMax = MaxIterations/RestartLength + ((MaxIterations%RestartLength == 0) ? 0 : 1);
+
+    for (int k=0; k<outerLoopMax; k++) {
 
       cp = outerLoopBody(LinOp, src, psi, rsq);
 
@@ -167,7 +169,7 @@ class GeneralisedMinimalResidual : public OperatorFunction<Field> {
       std::cout << GridLogIterative << "GeneralisedMinimalResidual: Iteration " << IterationCount
                 << " residual " << cp << " target " << rsq << std::endl;
 
-      if ((i == RestartLength - 1) || (cp <= rsq)) {
+      if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) {
 
         computeSolution(v, psi, i);