portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Fix bugs in Flexible GMRES solvers

Browse Source 
Somehow I got the left and right-preconditioned versions of GMRES mixed up. As
of now this is right-preconditioned version, which is what we want.
This commit is contained in:
Daniel Richtmann 2018-02-08 16:00:31 +01:00
parent 4ded1ceeb0
commit 9f79a87102
No known key found for this signature in database
GPG Key ID: B33C490AF0772057
2 changed files with 26 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
lib/algorithms/iterative/FlexibleCommunicationAvoidingGeneralisedMinimalResidual.h
View File

@ -148,18 +148,14 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF
Field w(src._grid); Field w(src._grid);
Field r(src._grid); Field r(src._grid);
Field z(src._grid);
std::vector<Field> v(RestartLength + 1, src._grid); std::vector<Field> v(RestartLength + 1, src._grid);
std::vector<Field> z(RestartLength + 1, src._grid);
MatrixTimer.Start(); MatrixTimer.Start();
LinOp.Op(psi, z); LinOp.Op(psi, w);
MatrixTimer.Stop(); MatrixTimer.Stop();
PrecTimer.Start();
Preconditioner(z, w);
PrecTimer.Stop();
LinalgTimer.Start(); LinalgTimer.Start();
r = src - w; r = src - w;
@ -183,7 +179,7 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF
if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) { if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) {
computeSolution(v, psi, i); computeSolution(z, psi, i);
return cp; return cp;
} }
@ -193,16 +189,16 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF
return cp; return cp;
} }
void arnoldiStep(LinearOperatorBase<Field> &LinOp, std::vector<Field> &v, Field &z, Field &w, int iter) { void arnoldiStep(LinearOperatorBase<Field> &LinOp, std::vector<Field> &v, std::vector<Field> &z, Field &w, int iter) {
MatrixTimer.Start();
LinOp.Op(v[iter], z);
MatrixTimer.Stop();
PrecTimer.Start(); PrecTimer.Start();
Preconditioner(z, w); Preconditioner(v[iter], z[iter]);
PrecTimer.Stop(); PrecTimer.Stop();
MatrixTimer.Start();
LinOp.Op(z[iter], w);
MatrixTimer.Stop();
LinalgTimer.Start(); LinalgTimer.Start();
for (int i = 0; i <= iter; ++i) { for (int i = 0; i <= iter; ++i) {
H(iter, i) = innerProduct(v[i], w); H(iter, i) = innerProduct(v[i], w);
@ -237,7 +233,7 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF
QrTimer.Stop(); QrTimer.Stop();
} }
void computeSolution(std::vector<Field> const &v, Field &psi, int iter) { void computeSolution(std::vector<Field> const &z, Field &psi, int iter) {
CompSolutionTimer.Start(); CompSolutionTimer.Start();
for (int i = iter; i >= 0; i--) { for (int i = iter; i >= 0; i--) {
@ -249,12 +245,12 @@ class FlexibleCommunicationAvoidingGeneralisedMinimalResidual : public OperatorF
if (true) { if (true) {
for (int i = 0; i <= iter; i++) for (int i = 0; i <= iter; i++)
psi = psi + v[i] * y[i]; psi = psi + z[i] * y[i];
} }
else { else {
psi = y[0] * v[0]; psi = y[0] * z[0];
for (int i = 1; i <= iter; i++) for (int i = 1; i <= iter; i++)
psi = psi + v[i] * y[i]; psi = psi + z[i] * y[i];
} }
CompSolutionTimer.Stop(); CompSolutionTimer.Stop();
} }

30
lib/algorithms/iterative/FlexibleGeneralisedMinimalResidual.h
View File

@ -148,18 +148,14 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction<Field> {
Field w(src._grid); Field w(src._grid);
Field r(src._grid); Field r(src._grid);
Field z(src._grid);
std::vector<Field> v(RestartLength + 1, src._grid); std::vector<Field> v(RestartLength + 1, src._grid);
std::vector<Field> z(RestartLength + 1, src._grid);
MatrixTimer.Start(); MatrixTimer.Start();
LinOp.Op(psi, z); LinOp.Op(psi, w);
MatrixTimer.Stop(); MatrixTimer.Stop();
PrecTimer.Start();
Preconditioner(z, w);
PrecTimer.Stop();
LinalgTimer.Start(); LinalgTimer.Start();
r = src - w; r = src - w;
@ -183,7 +179,7 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction<Field> {
if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) { if ((i == RestartLength - 1) || (IterationCount == MaxIterations) || (cp <= rsq)) {
computeSolution(v, psi, i); computeSolution(z, psi, i);
return cp; return cp;
} }
@ -193,16 +189,16 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction<Field> {
return cp; return cp;
} }
void arnoldiStep(LinearOperatorBase<Field> &LinOp, std::vector<Field> &v, Field &z, Field &w, int iter) { void arnoldiStep(LinearOperatorBase<Field> &LinOp, std::vector<Field> &v, std::vector<Field> &z, Field &w, int iter) {
MatrixTimer.Start();
LinOp.Op(v[iter], z);
MatrixTimer.Stop();
PrecTimer.Start(); PrecTimer.Start();
Preconditioner(z, w); Preconditioner(v[iter], z[iter]);
PrecTimer.Stop(); PrecTimer.Stop();
MatrixTimer.Start();
LinOp.Op(z[iter], w);
MatrixTimer.Stop();
LinalgTimer.Start(); LinalgTimer.Start();
for (int i = 0; i <= iter; ++i) { for (int i = 0; i <= iter; ++i) {
H(iter, i) = innerProduct(v[i], w); H(iter, i) = innerProduct(v[i], w);
@ -237,7 +233,7 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction<Field> {
QrTimer.Stop(); QrTimer.Stop();
} }
void computeSolution(std::vector<Field> const &v, Field &psi, int iter) { void computeSolution(std::vector<Field> const &z, Field &psi, int iter) {
CompSolutionTimer.Start(); CompSolutionTimer.Start();
for (int i = iter; i >= 0; i--) { for (int i = iter; i >= 0; i--) {
@ -249,12 +245,12 @@ class FlexibleGeneralisedMinimalResidual : public OperatorFunction<Field> {
if (true) { if (true) {
for (int i = 0; i <= iter; i++) for (int i = 0; i <= iter; i++)
psi = psi + v[i] * y[i]; psi = psi + z[i] * y[i];
} }
else { else {
psi = y[0] * v[0]; psi = y[0] * z[0];
for (int i = 1; i <= iter; i++) for (int i = 1; i <= iter; i++)
psi = psi + v[i] * y[i]; psi = psi + z[i] * y[i];
} }
CompSolutionTimer.Stop(); CompSolutionTimer.Stop();
} }