From 2ab9d4bc56a666b0a66b0de440106ec1bcb40fb8 Mon Sep 17 00:00:00 2001 From: Daniel Richtmann Date: Mon, 11 Jun 2018 15:54:32 +0200 Subject: [PATCH] WilsonMG: Fix random behavior in GMRES From time to time I saw random since the basis vectors were not initialized properly. --- .../CommunicationAvoidingGeneralisedMinimalResidual.h | 3 ++- ...FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h | 5 +++-- .../iterative/FlexibleGeneralisedMinimalResidual.h | 5 +++-- lib/algorithms/iterative/GeneralisedMinimalResidual.h | 3 ++- .../MixedPrecisionFlexibleGeneralisedMinimalResidual.h | 5 +++-- 5 files changed, 13 insertions(+), 8 deletions(-) diff --git a/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h b/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h index 1f5d293a..f0289683 100644 --- a/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/CommunicationAvoidingGeneralisedMinimalResidual.h @@ -145,7 +145,8 @@ class CommunicationAvoidingGeneralisedMinimalResidual : public OperatorFunction< Field w(src._grid); Field r(src._grid); - std::vector v(RestartLength + 1, src._grid); + // this should probably be made a class member so that it is only allocated once, not in every restart + std::vector v(RestartLength + 1, src._grid); for (auto &elem : v) elem = zero; MatrixTimer.Start(); LinOp.Op(psi, w); diff --git a/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h b/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h index b992f760..db857248 100644 --- a/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h @@ -152,8 +152,9 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF Field w(src._grid); Field r(src._grid); - std::vector v(RestartLength + 1, src._grid); - std::vector z(RestartLength + 1, src._grid); + // these should probably be made class members so that they are only allocated once, not in every restart + std::vector v(RestartLength + 1, src._grid); for (auto &elem : v) elem = zero; + std::vector z(RestartLength + 1, src._grid); for (auto &elem : z) elem = zero; MatrixTimer.Start(); LinOp.Op(psi, w); diff --git a/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h b/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h index bc5184d4..efc8c787 100644 --- a/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h @@ -150,8 +150,9 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction { Field w(src._grid); Field r(src._grid); - std::vector v(RestartLength + 1, src._grid); - std::vector z(RestartLength + 1, src._grid); + // these should probably be made class members so that they are only allocated once, not in every restart + std::vector v(RestartLength + 1, src._grid); for (auto &elem : v) elem = zero; + std::vector z(RestartLength + 1, src._grid); for (auto &elem : z) elem = zero; MatrixTimer.Start(); LinOp.Op(psi, w); diff --git a/lib/algorithms/iterative/GeneralisedMinimalResidual.h b/lib/algorithms/iterative/GeneralisedMinimalResidual.h index eaa43563..10636234 100644 --- a/lib/algorithms/iterative/GeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/GeneralisedMinimalResidual.h @@ -143,7 +143,8 @@ class GeneralisedMinimalResidual : public OperatorFunction { Field w(src._grid); Field r(src._grid); - std::vector v(RestartLength + 1, src._grid); + // this should probably be made a class member so that it is only allocated once, not in every restart + std::vector v(RestartLength + 1, src._grid); for (auto &elem : v) elem = zero; MatrixTimer.Start(); LinOp.Op(psi, w); diff --git a/lib/algorithms/iterative/MixedPrecisionFlexibleGeneralisedMinimalResidual.h b/lib/algorithms/iterative/MixedPrecisionFlexibleGeneralisedMinimalResidual.h index d38d41a9..04113684 100644 --- a/lib/algorithms/iterative/MixedPrecisionFlexibleGeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/MixedPrecisionFlexibleGeneralisedMinimalResidual.h @@ -157,8 +157,9 @@ class MixedPrecisionFlexibleGeneralisedMinimalResidual : public OperatorFunction FieldD w(src._grid); FieldD r(src._grid); - std::vector v(RestartLength + 1, src._grid); // these should probably be made class members - std::vector z(RestartLength + 1, src._grid); // so that they are only allocated once, not in every restart + // these should probably be made class members so that they are only allocated once, not in every restart + std::vector v(RestartLength + 1, src._grid); for (auto &elem : v) elem = zero; + std::vector z(RestartLength + 1, src._grid); for (auto &elem : z) elem = zero; MatrixTimer.Start(); LinOp.Op(psi, w);